src/c/qrdecomp.c

   1 /* QR Decomposition from LINPACK */
   2
   3 #include "linalg.h"
   4
   5 VOID linpack_dqrdc P8C(double *, x,
   6                        int, ldx,
   7                        int, n,
   8                        int, p,
   9                        double *, qraux,
  10                        int *, jpvt,
  11                        double *, work,
  12                        int, job)
  13 {
  14   /* System generated locals */
  15   int x_dim1, x_offset, i__1, i__2, i__3;
  16   double d__1, d__2;
  17
  18   /* Local variables */
  19   int negj;
  20   int maxj;
  21   int j, l;
  22   double t;
  23   int swapj;
  24   double nrmxl;
  25   int jj, jp, pl, pu;
  26   double tt, maxnrm;
  27   int lp1, lup;
  28
  29   /* dqrdc uses householder transformations to compute the qr */
  30   /* factorization of an n by p matrix x.  column pivoting */
  31   /* based on the 2-norms of the reduced columns may be */
  32   /* performed at the users option. */
  33
  34   /* on entry */
  35
  36   /*    x       double precision(ldx,p), where ldx .ge. n. */
  37   /*            x contains the matrix whose decomposition is to be */
  38   /*            computed. */
  39
  40   /*    ldx     integer. */
  41   /*            ldx is the leading dimension of the array x. */
  42
  43   /*    n       integer. */
  44   /*            n is the number of rows of the matrix x. */
  45
  46   /*    p       integer. */
  47   /*            p is the number of columns of the matrix x. */
  48
  49   /*    jpvt    integer(p). */
  50   /*            jpvt contains integers that control the selection */
  51   /*            of the pivot columns.  the k-th column x(k) of x */
  52   /*            is placed in one of three classes according to the */
  53   /*            value of jpvt(k). */
  54
  55   /*               if jpvt(k) .gt. 0, then x(k) is an initial */
  56   /*                                  column. */
  57
  58   /*               if jpvt(k) .eq. 0, then x(k) is a free column. */
  59
  60   /*               if jpvt(k) .lt. 0, then x(k) is a final column. */
  61
  62   /*            before the decomposition is computed, initial columns */
  63   /*            are moved to the beginning of the array x and final */
  64   /*            columns to the end.  both initial and final columns */
  65   /*            are frozen in place during the computation and only */
  66   /*            free columns are moved.  at the k-th stage of the */
  67   /*            reduction, if x(k) is occupied by a free column */
  68   /*            it is interchanged with the free column of largest */
  69   /*            reduced norm.  jpvt is not referenced if */
  70   /*            job .eq. 0. */
  71
  72   /*    work    double precision(p). */
  73   /*            work is a work array.  work is not referenced if */
  74   /*            job .eq. 0. */
  75
  76   /*    job     integer. */
  77   /*            job is an integer that initiates column pivoting. */
  78   /*            if job .eq. 0, no pivoting is done. */
  79   /*            if job .ne. 0, pivoting is done. */
  80
  81   /* on return */
  82
  83   /*    x       x contains in its upper triangle the upper */
  84   /*            triangular matrix r of the qr factorization. */
  85   /*            below its diagonal x contains information from */
  86   /*            which the orthogonal part of the decomposition */
  87   /*            can be recovered.  note that if pivoting has */
  88   /*            been requested, the decomposition is not that */
  89   /*            of the original matrix x but that of x */
  90   /*            with its columns permuted as described by jpvt. */
  91
  92   /*    qraux   double precision(p). */
  93   /*            qraux contains further information required to recover */
  94   /*            the orthogonal part of the decomposition. */
  95
  96   /*    jpvt    jpvt(k) contains the index of the column of the */
  97   /*            original matrix that has been interchanged into */
  98   /*            the k-th column, if pivoting was requested. */
  99
 100   /* linpack. this version dated 08/14/78 . */
 101   /* g.w. stewart, university of maryland, argonne national lab. */
 102
 103   /* dqrdc uses the following functions and subprograms. */
 104
 105   /* blas daxpy,ddot,dscal,dswap,dnrm2 */
 106   /* fortran dabs,dmax1,min0,dsqrt */
 107
 108   /* Parameter adjustments */
 109   x_dim1 = ldx;
 110   x_offset = x_dim1 + 1;
 111   x -= x_offset;
 112   --qraux;
 113   --jpvt;
 114   --work;
 115
 116   /* Function Body */
 117   pl = 1;
 118   pu = 0;
 119   if (job == 0) {
 120     goto L60;
 121   }
 122
 123   /*    pivoting has been requested.  rearrange the columns */
 124   /*    according to jpvt. */
 125
 126   i__1 = p;
 127   for (j = 1; j <= i__1; ++j) {
 128     swapj = jpvt[j] > 0;
 129     negj = jpvt[j] < 0;
 130     jpvt[j] = j;
 131     if (negj) {
 132       jpvt[j] = -j;
 133     }
 134     if (! swapj) {
 135       goto L10;
 136     }
 137     if (j != pl) {
 138       blas_dswap(n, &x[pl * x_dim1 + 1], 1, &x[j * x_dim1 + 1], 1);
 139     }
 140     jpvt[j] = jpvt[pl];
 141     jpvt[pl] = j;
 142     ++pl;
 143   L10:
 144     /* L20: */
 145     ;
 146   }
 147   pu = p;
 148   i__1 = p;
 149   for (jj = 1; jj <= i__1; ++jj) {
 150     j = p - jj + 1;
 151     if (jpvt[j] >= 0) {
 152       goto L40;
 153     }
 154     jpvt[j] = -jpvt[j];
 155     if (j == pu) {
 156       goto L30;
 157     }
 158     blas_dswap(n, &x[pu * x_dim1 + 1], 1, &x[j * x_dim1 + 1], 1);
 159     jp = jpvt[pu];
 160     jpvt[pu] = jpvt[j];
 161     jpvt[j] = jp;
 162   L30:
 163     --pu;
 164   L40:
 165     /* L50: */
 166     ;
 167   }
 168  L60:
 169
 170   /* compute the norms of the free columns. */
 171
 172   if (pu < pl) {
 173     goto L80;
 174   }
 175   i__1 = pu;
 176   for (j = pl; j <= i__1; ++j) {
 177     qraux[j] = blas_dnrm2(n, &x[j * x_dim1 + 1], 1);
 178     work[j] = qraux[j];
 179     /* L70: */
 180   }
 181  L80:
 182
 183   /* perform the householder reduction of x. */
 184
 185   lup = min(n,p);
 186   i__1 = lup;
 187   for (l = 1; l <= i__1; ++l) {
 188     if (l < pl || l >= pu) {
 189       goto L120;
 190     }
 191
 192     /*       locate the column of largest norm and bring it */
 193     /*       into the pivot position. */
 194
 195     maxnrm = 0.;
 196     maxj = l;
 197     i__2 = pu;
 198     for (j = l; j <= i__2; ++j) {
 199       if (qraux[j] <= maxnrm) {
 200         goto L90;
 201       }
 202       maxnrm = qraux[j];
 203       maxj = j;
 204     L90:
 205       /* L100: */
 206       ;
 207     }
 208     if (maxj == l) {
 209       goto L110;
 210     }
 211     blas_dswap(n, &x[l * x_dim1 + 1], 1, &x[maxj * x_dim1 + 1], 1);
 212     qraux[maxj] = qraux[l];
 213     work[maxj] = work[l];
 214     jp = jpvt[maxj];
 215     jpvt[maxj] = jpvt[l];
 216     jpvt[l] = jp;
 217   L110:
 218   L120:
 219     qraux[l] = 0.;
 220     if (l == n) {
 221       goto L190;
 222     }
 223
 224     /*       compute the householder transformation for column l. */
 225
 226     i__2 = n - l + 1;
 227     nrmxl = blas_dnrm2(i__2, &x[l + l * x_dim1], 1);
 228     if (nrmxl == 0.) {
 229       goto L180;
 230     }
 231     if (x[l + l * x_dim1] != 0.) {
 232       nrmxl = d_sign(&nrmxl, &x[l + l * x_dim1]);
 233     }
 234     i__2 = n - l + 1;
 235     d__1 = 1. / nrmxl;
 236     blas_dscal(i__2, d__1, &x[l + l * x_dim1], 1);
 237     x[l + l * x_dim1] += 1.;
 238
 239     /*          apply the transformation to the remaining columns, */
 240     /*          updating the norms. */
 241
 242     lp1 = l + 1;
 243     if (p < lp1) {
 244       goto L170;
 245     }
 246     i__2 = p;
 247     for (j = lp1; j <= i__2; ++j) {
 248       i__3 = n - l + 1;
 249       t = -blas_ddot(i__3, &x[l + l * x_dim1], 1, &x[l + j * x_dim1],
 250                      1) / x[l + l * x_dim1];
 251       i__3 = n - l + 1;
 252       blas_daxpy(i__3, t, &x[l + l * x_dim1], 1, &x[l + j * x_dim1], 1);
 253       if (j < pl || j > pu) {
 254         goto L150;
 255       }
 256       if (qraux[j] == 0.) {
 257         goto L150;
 258       }
 259       /* Computing 2nd power */
 260       d__2 = (d__1 = x[l + j * x_dim1], abs(d__1)) / qraux[j];
 261       tt = 1. - d__2 * d__2;
 262       tt = max(tt,0.);
 263       t = tt;
 264       /* Computing 2nd power */
 265       d__1 = qraux[j] / work[j];
 266       tt = tt * .05 * (d__1 * d__1) + 1.;
 267       if (tt == 1.) {
 268         goto L130;
 269       }
 270       qraux[j] *= sqrt(t);
 271       goto L140;
 272     L130:
 273       i__3 = n - l;
 274       qraux[j] = blas_dnrm2(i__3, &x[l + 1 + j * x_dim1], 1);
 275       work[j] = qraux[j];
 276     L140:
 277     L150:
 278       /* L160: */
 279       ;
 280     }
 281   L170:
 282
 283     /*          save the transformation. */
 284
 285     qraux[l] = x[l + l * x_dim1];
 286     x[l + l * x_dim1] = -nrmxl;
 287   L180:
 288   L190:
 289     /* L200: */
 290     ;
 291   }
 292   return;
 293 }
 294
 295 VOID linpack_dqrsl P13C(double *, x,
 296                         int, ldx,
 297                         int, n,
 298                         int, k,
 299                         double *, qraux,
 300                         double *, y,
 301                         double *, qy,
 302                         double *, qty,
 303                         double *, b,
 304                         double *, rsd,
 305                         double *, xb,
 306                         int, job,
 307                         int *, info)
 308 {
 309   /* System generated locals */
 310   integer x_dim1, x_offset, i__1, i__2;
 311
 312   /* Local variables */
 313   double temp;
 314   int cqty;
 315   int i, j;
 316   double t;
 317   int cb;
 318   int jj;
 319   int cr;
 320   int ju, kp1;
 321   int cxb, cqy;
 322
 323   /* dqrsl applies the output of dqrdc to compute coordinate */
 324   /* transformations, projections, and least squares solutions. */
 325   /* for k .le. min(n,p), let xk be the matrix */
 326
 327   /*        xk = (x(jpvt(1)),x(jpvt(2)), ... ,x(jpvt(k))) */
 328
 329   /* formed from columnns jpvt(1), ... ,jpvt(k) of the original */
 330   /* n x p matrix x that was input to dqrdc (if no pivoting was */
 331   /* done, xk consists of the first k columns of x in their */
 332   /* original order).  dqrdc produces a factored orthogonal matrix q */
 333   /* and an upper triangular matrix r such that */
 334
 335   /*          xk = q * (r) */
 336   /*                   (0) */
 337
 338   /* this information is contained in coded form in the arrays */
 339   /* x and qraux. */
 340
 341   /* on entry */
 342
 343   /*    x      double precision(ldx,p). */
 344   /*           x contains the output of dqrdc. */
 345
 346   /*    ldx    integer. */
 347   /*           ldx is the leading dimension of the array x. */
 348
 349   /*    n      integer. */
 350   /*           n is the number of rows of the matrix xk.  it must */
 351   /*           have the same value as n in dqrdc. */
 352
 353   /*    k      integer. */
 354   /*           k is the number of columns of the matrix xk.  k */
 355   /*           must nnot be greater than min(n,p), where p is the */
 356   /*           same as in the calling sequence to dqrdc. */
 357
 358   /*    qraux  double precision(p). */
 359   /*           qraux contains the auxiliary output from dqrdc. */
 360
 361   /*    y      double precision(n) */
 362   /*           y contains an n-vector that is to be manipulated */
 363   /*           by dqrsl. */
 364
 365   /*    job    integer. */
 366   /*           job specifies what is to be computed.  job has */
 367   /*           the decimal expansion abcde, with the following */
 368   /*           meaning. */
 369
 370   /*                if a.ne.0, compute qy. */
 371   /*                if b,c,d, or e .ne. 0, compute qty. */
 372   /*                if c.ne.0, compute b. */
 373   /*                if d.ne.0, compute rsd. */
 374   /*                if e.ne.0, compute xb. */
 375
 376   /*           note that a request to compute b, rsd, or xb */
 377   /*           automatically triggers the computation of qty, for */
 378   /*           which an array must be provided in the calling */
 379   /*           sequence. */
 380
 381   /* on return */
 382
 383   /*    qy     double precision(n). */
 384   /*           qy conntains q*y, if its computation has been */
 385   /*           requested. */
 386
 387   /*    qty    double precision(n). */
 388   /*           qty contains trans(q)*y, if its computation has */
 389   /*           been requested.  here trans(q) is the */
 390   /*           transpose of the matrix q. */
 391
 392   /*    b      double precision(k) */
 393   /*           b contains the solution of the least squares problem */
 394
 395   /*                minimize norm2(y - xk*b), */
 396
 397   /*           if its computation has been requested.  (note that */
 398   /*           if pivoting was requested in dqrdc, the j-th */
 399   /*           component of b will be associated with column jpvt(j) */
 400   /*           of the original matrix x that was input into dqrdc.) */
 401
 402   /*    rsd    double precision(n). */
 403   /*           rsd contains the least squares residual y - xk*b, */
 404   /*           if its computation has been requested.  rsd is */
 405   /*           also the orthogonal projection of y onto the */
 406   /*           orthogonal complement of the column space of xk. */
 407
 408   /*    xb     double precision(n). */
 409   /*           xb contains the least squares approximation xk*b, */
 410   /*           if its computation has been requested.  xb is also */
 411   /*           the orthogonal projection of y onto the column space */
 412   /*           of x. */
 413
 414   /*    info   integer. */
 415   /*           info is zero unless the computation of b has */
 416   /*           been requested and r is exactly singular.  in */
 417   /*           this case, info is the index of the first zero */
 418   /*           diagonal element of r and b is left unaltered. */
 419
 420   /* the parameters qy, qty, b, rsd, and xb are not referenced */
 421   /* if their computation is not requested and in this case */
 422   /* can be replaced by dummy variables in the calling program. */
 423   /* to save storage, the user may in some cases use the same */
 424   /* array for different parameters in the calling sequence.  a */
 425   /* frequently occuring example is when one wishes to compute */
 426   /* any of b, rsd, or xb and does not need y or qty.  in this */
 427   /* case one may identify y, qty, and one of b, rsd, or xb, while */
 428   /* providing separate arrays for anything else that is to be */
 429   /* computed.  thus the calling sequence */
 430
 431   /*      call dqrsl(x,ldx,n,k,qraux,y,dum,y,b,y,dum,110,info) */
 432
 433   /* will result in the computation of b and rsd, with rsd */
 434   /* overwriting y.  more generally, each item in the following */
 435   /* list contains groups of permissible identifications for */
 436   /* a single callinng sequence. */
 437
 438   /*      1. (y,qty,b) (rsd) (xb) (qy) */
 439
 440   /*      2. (y,qty,rsd) (b) (xb) (qy) */
 441
 442   /*      3. (y,qty,xb) (b) (rsd) (qy) */
 443
 444   /*      4. (y,qy) (qty,b) (rsd) (xb) */
 445
 446   /*      5. (y,qy) (qty,rsd) (b) (xb) */
 447
 448   /*      6. (y,qy) (qty,xb) (b) (rsd) */
 449
 450   /* in any group the value returned in the array allocated to */
 451   /* the group corresponds to the last member of the group. */
 452
 453   /* linpack. this version dated 08/14/78 . */
 454   /* g.w. stewart, university of maryland, argonne national lab. */
 455
 456   /* dqrsl uses the following functions and subprograms. */
 457
 458   /* blas daxpy,dcopy,ddot */
 459   /* fortran dabs,min0,mod */
 460
 461   /* Parameter adjustments */
 462   x_dim1 = ldx;
 463   x_offset = x_dim1 + 1;
 464   x -= x_offset;
 465   --qraux;
 466   --y;
 467   --qy;
 468   --qty;
 469   --b;
 470   --rsd;
 471   --xb;
 472
 473   /* set info flag. */
 474
 475   *info = 0;
 476
 477   /* determine what is to be computed. */
 478
 479   cqy = job / 10000 != 0;
 480   cqty = job % 10000 != 0;
 481   cb = job % 1000 / 100 != 0;
 482   cr = job % 100 / 10 != 0;
 483   cxb = job % 10 != 0;
 484   /* Computing MIN */
 485   i__1 = k, i__2 = n - 1;
 486   ju = min(i__1,i__2);
 487
 488   /* special action when n=1. */
 489
 490   if (ju != 0) {
 491     goto L40;
 492   }
 493   if (cqy) {
 494     qy[1] = y[1];
 495   }
 496   if (cqty) {
 497     qty[1] = y[1];
 498   }
 499   if (cxb) {
 500     xb[1] = y[1];
 501   }
 502   if (! cb) {
 503     goto L30;
 504   }
 505   if (x[x_dim1 + 1] != 0.) {
 506     goto L10;
 507   }
 508   *info = 1;
 509   goto L20;
 510  L10:
 511   b[1] = y[1] / x[x_dim1 + 1];
 512  L20:
 513  L30:
 514   if (cr) {
 515     rsd[1] = 0.;
 516   }
 517   goto L250;
 518  L40:
 519
 520   /*    set up to compute qy or qty. */
 521
 522   if (cqy) {
 523     blas_dcopy(n, &y[1], 1, &qy[1], 1);
 524   }
 525   if (cqty) {
 526     blas_dcopy(n, &y[1], 1, &qty[1], 1);
 527   }
 528   if (! cqy) {
 529     goto L70;
 530   }
 531
 532   /*       compute qy. */
 533
 534   i__1 = ju;
 535   for (jj = 1; jj <= i__1; ++jj) {
 536     j = ju - jj + 1;
 537     if (qraux[j] == 0.) {
 538       goto L50;
 539     }
 540     temp = x[j + j * x_dim1];
 541     x[j + j * x_dim1] = qraux[j];
 542     i__2 = n - j + 1;
 543     t = -blas_ddot(i__2, &x[j + j * x_dim1], 1, &qy[j], 1) / x[j + j
 544                                                                * x_dim1];
 545     i__2 = n - j + 1;
 546     blas_daxpy(i__2, t, &x[j + j * x_dim1], 1, &qy[j], 1);
 547     x[j + j * x_dim1] = temp;
 548   L50:
 549     /* L60: */
 550     ;
 551   }
 552  L70:
 553   if (! cqty) {
 554     goto L100;
 555   }
 556
 557   /*       compute trans(q)*y. */
 558
 559   i__1 = ju;
 560   for (j = 1; j <= i__1; ++j) {
 561     if (qraux[j] == 0.) {
 562       goto L80;
 563     }
 564     temp = x[j + j * x_dim1];
 565     x[j + j * x_dim1] = qraux[j];
 566     i__2 = n - j + 1;
 567     t = -blas_ddot(i__2, &x[j + j * x_dim1], 1, &qty[j], 1) / x[j +
 568                                                                 j * x_dim1];
 569     i__2 = n - j + 1;
 570     blas_daxpy(i__2, t, &x[j + j * x_dim1], 1, &qty[j], 1);
 571     x[j + j * x_dim1] = temp;
 572   L80:
 573     /* L90: */
 574     ;
 575   }
 576  L100:
 577
 578   /*    set up to compute b, rsd, or xb. */
 579
 580   if (cb) {
 581     blas_dcopy(k, &qty[1], 1, &b[1], 1);
 582   }
 583   kp1 = k + 1;
 584   if (cxb) {
 585     blas_dcopy(k, &qty[1], 1, &xb[1], 1);
 586   }
 587   if (cr && k < n) {
 588     i__1 = n - k;
 589     blas_dcopy(i__1, &qty[kp1], 1, &rsd[kp1], 1);
 590   }
 591   if (! cxb || kp1 > n) {
 592     goto L120;
 593   }
 594   i__1 = n;
 595   for (i = kp1; i <= i__1; ++i) {
 596     xb[i] = 0.;
 597     /* L110: */
 598   }
 599  L120:
 600   if (! cr) {
 601     goto L140;
 602   }
 603   i__1 = k;
 604   for (i = 1; i <= i__1; ++i) {
 605     rsd[i] = 0.;
 606     /* L130: */
 607   }
 608  L140:
 609   if (! cb) {
 610     goto L190;
 611   }
 612
 613   /*       compute b. */
 614
 615   i__1 = k;
 616   for (jj = 1; jj <= i__1; ++jj) {
 617     j = k - jj + 1;
 618     if (x[j + j * x_dim1] != 0.) {
 619       goto L150;
 620     }
 621     *info = j;
 622     /*       ......exit */
 623     goto L180;
 624   L150:
 625     b[j] /= x[j + j * x_dim1];
 626     if (j == 1) {
 627       goto L160;
 628     }
 629     t = -b[j];
 630     i__2 = j - 1;
 631     blas_daxpy(i__2, t, &x[j * x_dim1 + 1], 1, &b[1], 1);
 632   L160:
 633     /* L170: */
 634     ;
 635   }
 636  L180:
 637  L190:
 638   if (! cr && ! cxb) {
 639     goto L240;
 640   }
 641
 642   /*       compute rsd or xb as required. */
 643
 644   i__1 = ju;
 645   for (jj = 1; jj <= i__1; ++jj) {
 646     j = ju - jj + 1;
 647     if (qraux[j] == 0.) {
 648       goto L220;
 649     }
 650     temp = x[j + j * x_dim1];
 651     x[j + j * x_dim1] = qraux[j];
 652     if (! cr) {
 653       goto L200;
 654     }
 655     i__2 = n - j + 1;
 656     t = -blas_ddot(i__2, &x[j + j * x_dim1], 1, &rsd[j], 1) / x[j +
 657                                                                 j * x_dim1];
 658     i__2 = n - j + 1;
 659     blas_daxpy(i__2, t, &x[j + j * x_dim1], 1, &rsd[j], 1);
 660   L200:
 661     if (! cxb) {
 662       goto L210;
 663     }
 664     i__2 = n - j + 1;
 665     t = -blas_ddot(i__2, &x[j + j * x_dim1], 1, &xb[j], 1) / x[j + j
 666                                                                * x_dim1];
 667     i__2 = n - j + 1;
 668     blas_daxpy(i__2, t, &x[j + j * x_dim1], 1, &xb[j], 1);
 669   L210:
 670     x[j + j * x_dim1] = temp;
 671   L220:
 672     /* L230: */
 673     ;
 674   }
 675  L240:
 676  L250:
 677   return;
 678 }
 679
 680 VOID linpack_zqrdc P8C(dcomplex *, x,
 681                        int, ldx,
 682                        int, n,
 683                        int, p,
 684                        dcomplex *, qraux,
 685                        int *, jpvt,
 686                        dcomplex *, work,
 687                        int, job)
 688 {
 689   /* System generated locals */
 690   int x_dim1, x_offset, i__1, i__2, i__3, i__4;
 691   double d__1, d__2, d__3, d__4;
 692   dcomplex z__1, z__2, z__3;
 693   static dcomplex c_b28 = {1.,0.};
 694
 695   /* Local variables */
 696   int negj;
 697   int maxj, j, l;
 698   dcomplex t;
 699   int swapj;
 700   dcomplex nrmxl;
 701   int jj, jp, pl, pu;
 702   double tt;
 703   double maxnrm;
 704   int lp1, lup;
 705
 706
 707   /*     zqrdc uses householder transformations to compute the qr */
 708   /*     factorization of an n by p matrix x.  column pivoting */
 709   /*     based on the 2-norms of the reduced columns may be */
 710   /*     performed at the users option. */
 711
 712   /*     on entry */
 713
 714   /*        x       complex*16(ldx,p), where ldx .ge. n. */
 715   /*                x contains the matrix whose decomposition is to be */
 716   /*                computed. */
 717
 718   /*        ldx     integer. */
 719   /*                ldx is the leading dimension of the array x. */
 720
 721   /*        n       integer. */
 722   /*                n is the number of rows of the matrix x. */
 723
 724   /*        p       integer. */
 725   /*                p is the number of columns of the matrix x. */
 726
 727   /*        jpvt    integer(p). */
 728   /*                jpvt contains integers that control the selection */
 729   /*                of the pivot columns.  the k-th column x(k) of x */
 730   /*                is placed in one of three classes according to the */
 731   /*                value of jpvt(k). */
 732
 733   /*                   if jpvt(k) .gt. 0, then x(k) is an initial */
 734   /*                                      column. */
 735
 736   /*                   if jpvt(k) .eq. 0, then x(k) is a free column. */
 737
 738   /*                   if jpvt(k) .lt. 0, then x(k) is a final column. */
 739
 740   /*                before the decomposition is computed, initial columns */
 741   /*                are moved to the beginning of the array x and final */
 742   /*                columns to the end.  both initial and final columns */
 743   /*                are frozen in place during the computation and only */
 744   /*                free columns are moved.  at the k-th stage of the */
 745   /*                reduction, if x(k) is occupied by a free column */
 746   /*                it is interchanged with the free column of largest */
 747   /*                reduced norm.  jpvt is not referenced if */
 748   /*                job .eq. 0. */
 749
 750   /*        work    complex*16(p). */
 751   /*                work is a work array.  work is not referenced if */
 752   /*                job .eq. 0. */
 753
 754   /*        job     integer. */
 755   /*                job is an integer that initiates column pivoting. */
 756   /*                if job .eq. 0, no pivoting is done. */
 757   /*                if job .ne. 0, pivoting is done. */
 758
 759   /*     on return */
 760
 761   /*        x       x contains in its upper triangle the upper */
 762   /*                triangular matrix r of the qr factorization. */
 763   /*                below its diagonal x contains information from */
 764   /*                which the unitary part of the decomposition */
 765   /*                can be recovered.  note that if pivoting has */
 766   /*                been requested, the decomposition is not that */
 767   /*                of the original matrix x but that of x */
 768   /*                with its columns permuted as described by jpvt. */
 769
 770   /*        qraux   complex*16(p). */
 771   /*                qraux contains further information required to recover
 772    */
 773   /*                the unitary part of the decomposition. */
 774
 775   /*        jpvt    jpvt(k) contains the index of the column of the */
 776   /*                original matrix that has been interchanged into */
 777   /*                the k-th column, if pivoting was requested. */
 778
 779   /*     linpack. this version dated 08/14/78 . */
 780   /*     g.w. stewart, university of maryland, argonne national lab. */
 781
 782   /*     zqrdc uses the following functions and subprograms. */
 783
 784   /*     blas zaxpy,zdotc,zscal,zswap,dznrm2 */
 785   /*     fortran dabs,dmax1,cdabs,dcmplx,cdsqrt,min0 */
 786
 787   /*     internal variables */
 788
 789
 790
 791   /* Parameter adjustments */
 792   x_dim1 = ldx;
 793   x_offset = x_dim1 + 1;
 794   x -= x_offset;
 795   --qraux;
 796   --jpvt;
 797   --work;
 798
 799   /* Function Body */
 800   pl = 1;
 801   pu = 0;
 802   if (job == 0) {
 803     goto L60;
 804   }
 805
 806   /*        pivoting has been requested.  rearrange the columns */
 807   /*        according to jpvt. */
 808
 809   i__1 = p;
 810   for (j = 1; j <= i__1; ++j) {
 811     swapj = jpvt[j] > 0;
 812     negj = jpvt[j] < 0;
 813     jpvt[j] = j;
 814     if (negj) {
 815       jpvt[j] = -j;
 816     }
 817     if (! swapj) {
 818       goto L10;
 819     }
 820     if (j != pl) {
 821       blas_zswap(n, &x[pl * x_dim1 + 1], 1, &x[j * x_dim1 + 1], 1);
 822     }
 823     jpvt[j] = jpvt[pl];
 824     jpvt[pl] = j;
 825     ++pl;
 826   L10:
 827     /* L20: */
 828     ;
 829   }
 830   pu = p;
 831   i__1 = p;
 832   for (jj = 1; jj <= i__1; ++jj) {
 833     j = p - jj + 1;
 834     if (jpvt[j] >= 0) {
 835       goto L40;
 836     }
 837     jpvt[j] = -jpvt[j];
 838     if (j == pu) {
 839       goto L30;
 840     }
 841     blas_zswap(n, &x[pu * x_dim1 + 1], 1, &x[j * x_dim1 + 1], 1);
 842     jp = jpvt[pu];
 843     jpvt[pu] = jpvt[j];
 844     jpvt[j] = jp;
 845   L30:
 846     --pu;
 847   L40:
 848     /* L50: */
 849     ;
 850   }
 851  L60:
 852
 853   /*     compute the norms of the free columns. */
 854
 855   if (pu < pl) {
 856     goto L80;
 857   }
 858   i__1 = pu;
 859   for (j = pl; j <= i__1; ++j) {
 860     i__2 = j;
 861     d__1 = blas_dznrm2(n, &x[j * x_dim1 + 1], 1);
 862     z__1.r = d__1, z__1.i = 0.;
 863     qraux[i__2].r = z__1.r, qraux[i__2].i = z__1.i;
 864     i__2 = j;
 865     i__3 = j;
 866     work[i__2].r = qraux[i__3].r, work[i__2].i = qraux[i__3].i;
 867     /* L70: */
 868   }
 869  L80:
 870
 871   /*     perform the householder reduction of x. */
 872
 873   lup = min(n,p);
 874   i__1 = lup;
 875   for (l = 1; l <= i__1; ++l) {
 876     if (l < pl || l >= pu) {
 877       goto L120;
 878     }
 879
 880     /*           locate the column of largest norm and bring it */
 881     /*           into the pivot position. */
 882
 883     maxnrm = 0.;
 884     maxj = l;
 885     i__2 = pu;
 886     for (j = l; j <= i__2; ++j) {
 887       i__3 = j;
 888       if (qraux[i__3].r <= maxnrm) {
 889         goto L90;
 890       }
 891       i__3 = j;
 892       maxnrm = qraux[i__3].r;
 893       maxj = j;
 894     L90:
 895       /* L100: */
 896       ;
 897     }
 898     if (maxj == l) {
 899       goto L110;
 900     }
 901     blas_zswap(n, &x[l * x_dim1 + 1], 1, &x[maxj * x_dim1 + 1], 1);
 902     i__2 = maxj;
 903     i__3 = l;
 904     qraux[i__2].r = qraux[i__3].r, qraux[i__2].i = qraux[i__3].i;
 905     i__2 = maxj;
 906     i__3 = l;
 907     work[i__2].r = work[i__3].r, work[i__2].i = work[i__3].i;
 908     jp = jpvt[maxj];
 909     jpvt[maxj] = jpvt[l];
 910     jpvt[l] = jp;
 911   L110:
 912   L120:
 913     i__2 = l;
 914     qraux[i__2].r = 0., qraux[i__2].i = 0.;
 915     if (l == n) {
 916       goto L190;
 917     }
 918
 919     /*           compute the householder transformation for column l. */
 920
 921     i__2 = n - l + 1;
 922     d__1 = blas_dznrm2(i__2, &x[l + l * x_dim1], 1);
 923     z__1.r = d__1, z__1.i = 0.;
 924     nrmxl.r = z__1.r, nrmxl.i = z__1.i;
 925     z__1.r = nrmxl.r * 0. - nrmxl.i * -1., z__1.i = nrmxl.r * -1. +
 926       nrmxl.i * 0.;
 927     if ((d__1 = nrmxl.r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) == 0.) {
 928       goto L180;
 929     }
 930     i__2 = l + l * x_dim1;
 931     i__3 = l + l * x_dim1;
 932     z__1.r = x[i__3].r * 0. - x[i__3].i * -1., z__1.i = x[i__3].r * -1. +
 933       x[i__3].i * 0.;
 934     if ((d__1 = x[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) != 0.)
 935       {
 936         d__3 = z_abs(&nrmxl);
 937         i__4 = l + l * x_dim1;
 938         d__4 = z_abs(&x[l + l * x_dim1]);
 939         z__3.r = x[i__4].r / d__4, z__3.i = x[i__4].i / d__4;
 940         z__2.r = d__3 * z__3.r, z__2.i = d__3 * z__3.i;
 941         nrmxl.r = z__2.r, nrmxl.i = z__2.i;
 942       }
 943     i__2 = n - l + 1;
 944     z_div(&z__1, &c_b28, &nrmxl);
 945     blas_zscal(i__2, &z__1, &x[l + l * x_dim1], 1);
 946     i__2 = l + l * x_dim1;
 947     i__3 = l + l * x_dim1;
 948     z__1.r = x[i__3].r + 1., z__1.i = x[i__3].i + 0.;
 949     x[i__2].r = z__1.r, x[i__2].i = z__1.i;
 950
 951     /*              apply the transformation to the remaining columns, */
 952     /*              updating the norms. */
 953
 954     lp1 = l + 1;
 955     if (p < lp1) {
 956       goto L170;
 957     }
 958     i__2 = p;
 959     for (j = lp1; j <= i__2; ++j) {
 960       i__3 = n - l + 1;
 961       blas_zdotc(&z__3, i__3, &x[l + l * x_dim1], 1, &x[l + j * x_dim1], 1);
 962       z__2.r = -z__3.r, z__2.i = -z__3.i;
 963       z_div(&z__1, &z__2, &x[l + l * x_dim1]);
 964       t.r = z__1.r, t.i = z__1.i;
 965       i__3 = n - l + 1;
 966       blas_zaxpy(i__3, &t, &x[l + l * x_dim1], 1, &x[l + j * x_dim1], 1);
 967       if (j < pl || j > pu) {
 968         goto L150;
 969       }
 970       i__3 = j;
 971       i__4 = j;
 972       z__1.r = qraux[i__4].r * 0. - qraux[i__4].i * -1., z__1.i = qraux[
 973                                                                         i__4].r * -1. + qraux[i__4].i * 0.;
 974       if ((d__1 = qraux[i__3].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2))
 975           == 0.) {
 976         goto L150;
 977       }
 978       i__3 = j;
 979       /* Computing 2nd power */
 980       d__1 = z_abs(&x[l + j * x_dim1]) / qraux[i__3].r;
 981       tt = 1. - d__1 * d__1;
 982       tt = max(tt,0.);
 983       z__1.r = tt, z__1.i = 0.;
 984       t.r = z__1.r, t.i = z__1.i;
 985       i__3 = j;
 986       i__4 = j;
 987       /* Computing 2nd power */
 988       d__1 = qraux[i__3].r / work[i__4].r;
 989       tt = tt * .05 * (d__1 * d__1) + 1.;
 990       if (tt == 1.) {
 991         goto L130;
 992       }
 993       i__3 = j;
 994       i__4 = j;
 995       z_sqrt(&z__2, &t);
 996       z__1.r = qraux[i__4].r * z__2.r - qraux[i__4].i * z__2.i, z__1.i =
 997         qraux[i__4].r * z__2.i + qraux[i__4].i * z__2.r;
 998       qraux[i__3].r = z__1.r, qraux[i__3].i = z__1.i;
 999       goto L140;
1000     L130:
1001       i__3 = j;
1002       i__4 = n - l;
1003       d__1 = blas_dznrm2(i__4, &x[l + 1 + j * x_dim1], 1);
1004       z__1.r = d__1, z__1.i = 0.;
1005       qraux[i__3].r = z__1.r, qraux[i__3].i = z__1.i;
1006       i__3 = j;
1007       i__4 = j;
1008       work[i__3].r = qraux[i__4].r, work[i__3].i = qraux[i__4].i;
1009     L140:
1010     L150:
1011       /* L160: */
1012       ;
1013     }
1014   L170:
1015
1016     /*              save the transformation. */
1017
1018     i__2 = l;
1019     i__3 = l + l * x_dim1;
1020     qraux[i__2].r = x[i__3].r, qraux[i__2].i = x[i__3].i;
1021     i__2 = l + l * x_dim1;
1022     z__1.r = -nrmxl.r, z__1.i = -nrmxl.i;
1023     x[i__2].r = z__1.r, x[i__2].i = z__1.i;
1024   L180:
1025   L190:
1026     /* L200: */
1027     ;
1028   }
1029   return;
1030 }
1031
1032 VOID linpack_zqrsl P13C(dcomplex *, x,
1033                         int, ldx,
1034                         int, n,
1035                         int, k,
1036                         dcomplex *, qraux,
1037                         dcomplex *, y,
1038                         dcomplex *, qy,
1039                         dcomplex *, qty,
1040                         dcomplex *, b,
1041                         dcomplex *, rsd,
1042                         dcomplex *, xb,
1043                         int, job,
1044                         int *, info)
1045 {
1046   /* System generated locals */
1047   int x_dim1, x_offset, i__1, i__2, i__3;
1048   double d__1, d__2;
1049   dcomplex z__1, z__2, z__3;
1050
1051   /* Local variables */
1052   dcomplex temp;
1053   int cqty;
1054   int i, j;
1055   dcomplex t;
1056   int cb;
1057   int jj;
1058   int cr;
1059   int ju, kp1;
1060   int cxb, cqy;
1061
1062
1063   /*     zqrsl applies the output of zqrdc to compute coordinate */
1064   /*     transformations, projections, and least squares solutions. */
1065   /*     for k .le. min(n,p), let xk be the matrix */
1066
1067   /*            xk = (x(jpvt(1)),x(jpvt(2)), ... ,x(jpvt(k))) */
1068
1069   /*     formed from columnns jpvt(1), ... ,jpvt(k) of the original */
1070   /*     n x p matrix x that was input to zqrdc (if no pivoting was */
1071   /*     done, xk consists of the first k columns of x in their */
1072   /*     original order).  zqrdc produces a factored unitary matrix q */
1073   /*     and an upper triangular matrix r such that */
1074
1075   /*              xk = q * (r) */
1076   /*                       (0) */
1077
1078   /*     this information is contained in coded form in the arrays */
1079   /*     x and qraux. */
1080
1081   /*     on entry */
1082
1083   /*        x      complex*16(ldx,p). */
1084   /*               x contains the output of zqrdc. */
1085
1086   /*        ldx    integer. */
1087   /*               ldx is the leading dimension of the array x. */
1088
1089   /*        n      integer. */
1090   /*               n is the number of rows of the matrix xk.  it must */
1091   /*               have the same value as n in zqrdc. */
1092
1093   /*        k      integer. */
1094   /*               k is the number of columns of the matrix xk.  k */
1095   /*               must nnot be greater than min(n,p), where p is the */
1096   /*               same as in the calling sequence to zqrdc. */
1097
1098   /*        qraux  complex*16(p). */
1099   /*               qraux contains the auxiliary output from zqrdc. */
1100
1101   /*        y      complex*16(n) */
1102   /*               y contains an n-vector that is to be manipulated */
1103   /*               by zqrsl. */
1104
1105   /*        job    integer. */
1106   /*               job specifies what is to be computed.  job has */
1107   /*               the decimal expansion abcde, with the following */
1108   /*               meaning. */
1109
1110   /*                    if a.ne.0, compute qy. */
1111   /*                    if b,c,d, or e .ne. 0, compute qty. */
1112   /*                    if c.ne.0, compute b. */
1113   /*                    if d.ne.0, compute rsd. */
1114   /*                    if e.ne.0, compute xb. */
1115
1116   /*               note that a request to compute b, rsd, or xb */
1117   /*               automatically triggers the computation of qty, for */
1118   /*               which an array must be provided in the calling */
1119   /*               sequence. */
1120
1121   /*     on return */
1122
1123   /*        qy     complex*16(n). */
1124   /*               qy conntains q*y, if its computation has been */
1125   /*               requested. */
1126
1127   /*        qty    complex*16(n). */
1128   /*               qty contains ctrans(q)*y, if its computation has */
1129   /*               been requested.  here ctrans(q) is the conjugate */
1130   /*               transpose of the matrix q. */
1131
1132   /*        b      complex*16(k) */
1133   /*               b contains the solution of the least squares problem */
1134
1135   /*                    minimize norm2(y - xk*b), */
1136
1137   /*               if its computation has been requested.  (note that */
1138   /*               if pivoting was requested in zqrdc, the j-th */
1139   /*               component of b will be associated with column jpvt(j) */
1140   /*               of the original matrix x that was input into zqrdc.) */
1141
1142   /*        rsd    complex*16(n). */
1143   /*               rsd contains the least squares residual y - xk*b, */
1144   /*               if its computation has been requested.  rsd is */
1145   /*               also the orthogonal projection of y onto the */
1146   /*               orthogonal complement of the column space of xk. */
1147
1148   /*        xb     complex*16(n). */
1149   /*               xb contains the least squares approximation xk*b, */
1150   /*               if its computation has been requested.  xb is also */
1151   /*               the orthogonal projection of y onto the column space */
1152   /*               of x. */
1153
1154   /*        info   integer. */
1155   /*               info is zero unless the computation of b has */
1156   /*               been requested and r is exactly singular.  in */
1157   /*               this case, info is the index of the first zero */
1158   /*               diagonal element of r and b is left unaltered. */
1159
1160   /*     the parameters qy, qty, b, rsd, and xb are not referenced */
1161   /*     if their computation is not requested and in this case */
1162   /*     can be replaced by dummy variables in the calling program. */
1163   /*     to save storage, the user may in some cases use the same */
1164   /*     array for different parameters in the calling sequence.  a */
1165   /*     frequently occuring example is when one wishes to compute */
1166   /*     any of b, rsd, or xb and does not need y or qty.  in this */
1167   /*     case one may identify y, qty, and one of b, rsd, or xb, while */
1168   /*     providing separate arrays for anything else that is to be */
1169   /*     computed.  thus the calling sequence */
1170
1171   /*          call zqrsl(x,ldx,n,k,qraux,y,dum,y,b,y,dum,110,info) */
1172
1173   /*     will result in the computation of b and rsd, with rsd */
1174   /*     overwriting y.  more generally, each item in the following */
1175   /*     list contains groups of permissible identifications for */
1176   /*     a single callinng sequence. */
1177
1178   /*          1. (y,qty,b) (rsd) (xb) (qy) */
1179
1180   /*          2. (y,qty,rsd) (b) (xb) (qy) */
1181
1182   /*          3. (y,qty,xb) (b) (rsd) (qy) */
1183
1184   /*          4. (y,qy) (qty,b) (rsd) (xb) */
1185
1186   /*          5. (y,qy) (qty,rsd) (b) (xb) */
1187
1188   /*          6. (y,qy) (qty,xb) (b) (rsd) */
1189
1190   /*     in any group the value returned in the array allocated to */
1191   /*     the group corresponds to the last member of the group. */
1192
1193   /*     linpack. this version dated 08/14/78 . */
1194   /*     g.w. stewart, university of maryland, argonne national lab. */
1195
1196   /*     zqrsl uses the following functions and subprograms. */
1197
1198   /*     blas zaxpy,zcopy,zdotc */
1199   /*     fortran dabs,min0,mod */
1200
1201   /*     internal variables */
1202
1203
1204
1205   /*     set info flag. */
1206
1207   /* Parameter adjustments */
1208   x_dim1 = ldx;
1209   x_offset = x_dim1 + 1;
1210   x -= x_offset;
1211   --qraux;
1212   --y;
1213   --qy;
1214   --qty;
1215   --b;
1216   --rsd;
1217   --xb;
1218
1219   /* Function Body */
1220   *info = 0;
1221
1222   /*     determine what is to be computed. */
1223
1224   cqy = job / 10000 != 0;
1225   cqty = job % 10000 != 0;
1226   cb = job % 1000 / 100 != 0;
1227   cr = job % 100 / 10 != 0;
1228   cxb = job % 10 != 0;
1229   /* Computing MIN */
1230   i__1 = k, i__2 = n - 1;
1231   ju = min(i__1,i__2);
1232
1233   /*     special action when n=1. */
1234
1235   if (ju != 0) {
1236     goto L40;
1237   }
1238   if (cqy) {
1239     qy[1].r = y[1].r, qy[1].i = y[1].i;
1240   }
1241   if (cqty) {
1242     qty[1].r = y[1].r, qty[1].i = y[1].i;
1243   }
1244   if (cxb) {
1245     xb[1].r = y[1].r, xb[1].i = y[1].i;
1246   }
1247   if (! cb) {
1248     goto L30;
1249   }
1250   i__1 = x_dim1 + 1;
1251   i__2 = x_dim1 + 1;
1252   z__1.r = x[i__2].r * 0. - x[i__2].i * -1., z__1.i = x[i__2].r * -1. + x[
1253                                                                           i__2].i * 0.;
1254   if ((d__1 = x[i__1].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) != 0.) {
1255     goto L10;
1256   }
1257   *info = 1;
1258   goto L20;
1259  L10:
1260   z_div(&z__1, &y[1], &x[x_dim1 + 1]);
1261   b[1].r = z__1.r, b[1].i = z__1.i;
1262  L20:
1263  L30:
1264   if (cr) {
1265     rsd[1].r = 0., rsd[1].i = 0.;
1266   }
1267   goto L250;
1268  L40:
1269
1270   /*        set up to compute qy or qty. */
1271
1272   if (cqy) {
1273     blas_zcopy(n, &y[1], 1, &qy[1], 1);
1274   }
1275   if (cqty) {
1276     blas_zcopy(n, &y[1], 1, &qty[1], 1);
1277   }
1278   if (! cqy) {
1279     goto L70;
1280   }
1281
1282   /*           compute qy. */
1283
1284   i__1 = ju;
1285   for (jj = 1; jj <= i__1; ++jj) {
1286     j = ju - jj + 1;
1287     i__2 = j;
1288     i__3 = j;
1289     z__1.r = qraux[i__3].r * 0. - qraux[i__3].i * -1., z__1.i = qraux[
1290                                                                       i__3].r * -1. + qraux[i__3].i * 0.;
1291     if ((d__1 = qraux[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) ==
1292         0.) {
1293       goto L50;
1294     }
1295     i__2 = j + j * x_dim1;
1296     temp.r = x[i__2].r, temp.i = x[i__2].i;
1297     i__2 = j + j * x_dim1;
1298     i__3 = j;
1299     x[i__2].r = qraux[i__3].r, x[i__2].i = qraux[i__3].i;
1300     i__2 = n - j + 1;
1301     blas_zdotc(&z__3, i__2, &x[j + j * x_dim1], 1, &qy[j], 1);
1302     z__2.r = -z__3.r, z__2.i = -z__3.i;
1303     z_div(&z__1, &z__2, &x[j + j * x_dim1]);
1304     t.r = z__1.r, t.i = z__1.i;
1305     i__2 = n - j + 1;
1306     blas_zaxpy(i__2, &t, &x[j + j * x_dim1], 1, &qy[j], 1);
1307     i__2 = j + j * x_dim1;
1308     x[i__2].r = temp.r, x[i__2].i = temp.i;
1309   L50:
1310     /* L60: */
1311     ;
1312   }
1313  L70:
1314   if (! cqty) {
1315     goto L100;
1316   }
1317
1318   /*           compute ctrans(q)*y. */
1319
1320   i__1 = ju;
1321   for (j = 1; j <= i__1; ++j) {
1322     i__2 = j;
1323     i__3 = j;
1324     z__1.r = qraux[i__3].r * 0. - qraux[i__3].i * -1., z__1.i = qraux[
1325                                                                       i__3].r * -1. + qraux[i__3].i * 0.;
1326     if ((d__1 = qraux[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) ==
1327         0.) {
1328       goto L80;
1329     }
1330     i__2 = j + j * x_dim1;
1331     temp.r = x[i__2].r, temp.i = x[i__2].i;
1332     i__2 = j + j * x_dim1;
1333     i__3 = j;
1334     x[i__2].r = qraux[i__3].r, x[i__2].i = qraux[i__3].i;
1335     i__2 = n - j + 1;
1336     blas_zdotc(&z__3, i__2, &x[j + j * x_dim1], 1, &qty[j], 1);
1337     z__2.r = -z__3.r, z__2.i = -z__3.i;
1338     z_div(&z__1, &z__2, &x[j + j * x_dim1]);
1339     t.r = z__1.r, t.i = z__1.i;
1340     i__2 = n - j + 1;
1341     blas_zaxpy(i__2, &t, &x[j + j * x_dim1], 1, &qty[j], 1);
1342     i__2 = j + j * x_dim1;
1343     x[i__2].r = temp.r, x[i__2].i = temp.i;
1344   L80:
1345     /* L90: */
1346     ;
1347   }
1348  L100:
1349
1350   /*        set up to compute b, rsd, or xb. */
1351
1352   if (cb) {
1353     blas_zcopy(k, &qty[1], 1, &b[1], 1);
1354   }
1355   kp1 = k + 1;
1356   if (cxb) {
1357     blas_zcopy(k, &qty[1], 1, &xb[1], 1);
1358   }
1359   if (cr && k < n) {
1360     i__1 = n - k;
1361     blas_zcopy(i__1, &qty[kp1], 1, &rsd[kp1], 1);
1362   }
1363   if (! cxb || kp1 > n) {
1364     goto L120;
1365   }
1366   i__1 = n;
1367   for (i = kp1; i <= i__1; ++i) {
1368     i__2 = i;
1369     xb[i__2].r = 0., xb[i__2].i = 0.;
1370     /* L110: */
1371   }
1372  L120:
1373   if (! cr) {
1374     goto L140;
1375   }
1376   i__1 = k;
1377   for (i = 1; i <= i__1; ++i) {
1378     i__2 = i;
1379     rsd[i__2].r = 0., rsd[i__2].i = 0.;
1380     /* L130: */
1381   }
1382  L140:
1383   if (! cb) {
1384     goto L190;
1385   }
1386
1387   /*           compute b. */
1388
1389   i__1 = k;
1390   for (jj = 1; jj <= i__1; ++jj) {
1391     j = k - jj + 1;
1392     i__2 = j + j * x_dim1;
1393     i__3 = j + j * x_dim1;
1394     z__1.r = x[i__3].r * 0. - x[i__3].i * -1., z__1.i = x[i__3].r * -1. +
1395       x[i__3].i * 0.;
1396     if ((d__1 = x[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) != 0.)
1397       {
1398         goto L150;
1399       }
1400     *info = j;
1401     /*           ......exit */
1402     goto L180;
1403   L150:
1404     i__2 = j;
1405     z_div(&z__1, &b[j], &x[j + j * x_dim1]);
1406     b[i__2].r = z__1.r, b[i__2].i = z__1.i;
1407     if (j == 1) {
1408       goto L160;
1409     }
1410     i__2 = j;
1411     z__1.r = -b[i__2].r, z__1.i = -b[i__2].i;
1412     t.r = z__1.r, t.i = z__1.i;
1413     i__2 = j - 1;
1414     blas_zaxpy(i__2, &t, &x[j * x_dim1 + 1], 1, &b[1], 1);
1415   L160:
1416     /* L170: */
1417     ;
1418   }
1419  L180:
1420  L190:
1421   if (! cr && ! cxb) {
1422     goto L240;
1423   }
1424
1425   /*           compute rsd or xb as required. */
1426
1427   i__1 = ju;
1428   for (jj = 1; jj <= i__1; ++jj) {
1429     j = ju - jj + 1;
1430     i__2 = j;
1431     i__3 = j;
1432     z__1.r = qraux[i__3].r * 0. - qraux[i__3].i * -1., z__1.i = qraux[
1433                                                                       i__3].r * -1. + qraux[i__3].i * 0.;
1434     if ((d__1 = qraux[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) ==
1435         0.) {
1436       goto L220;
1437     }
1438     i__2 = j + j * x_dim1;
1439     temp.r = x[i__2].r, temp.i = x[i__2].i;
1440     i__2 = j + j * x_dim1;
1441     i__3 = j;
1442     x[i__2].r = qraux[i__3].r, x[i__2].i = qraux[i__3].i;
1443     if (! cr) {
1444       goto L200;
1445     }
1446     i__2 = n - j + 1;
1447     blas_zdotc(&z__3, i__2, &x[j + j * x_dim1], 1, &rsd[j], 1);
1448     z__2.r = -z__3.r, z__2.i = -z__3.i;
1449     z_div(&z__1, &z__2, &x[j + j * x_dim1]);
1450     t.r = z__1.r, t.i = z__1.i;
1451     i__2 = n - j + 1;
1452     blas_zaxpy(i__2, &t, &x[j + j * x_dim1], 1, &rsd[j], 1);
1453   L200:
1454     if (! cxb) {
1455       goto L210;
1456     }
1457     i__2 = n - j + 1;
1458     blas_zdotc(&z__3, i__2, &x[j + j * x_dim1], 1, &xb[j], 1);
1459     z__2.r = -z__3.r, z__2.i = -z__3.i;
1460     z_div(&z__1, &z__2, &x[j + j * x_dim1]);
1461     t.r = z__1.r, t.i = z__1.i;
1462     i__2 = n - j + 1;
1463     blas_zaxpy(i__2, &t, &x[j + j * x_dim1], 1, &xb[j], 1);
1464   L210:
1465     i__2 = j + j * x_dim1;
1466     x[i__2].r = temp.r, x[i__2].i = temp.i;
1467   L220:
1468     /* L230: */
1469     ;
1470   }
1471  L240:
1472  L250:
1473   return;
1474 }