/* QR Decomposition from LINPACK */
#include "linalg.h"
VOID linpack_dqrdc P8C(double *, x,
int, ldx,
int, n,
int, p,
double *, qraux,
int *, jpvt,
double *, work,
int, job)
{
/* System generated locals */
int x_dim1, x_offset, i__1, i__2, i__3;
double d__1, d__2;
/* Local variables */
int negj;
int maxj;
int j, l;
double t;
int swapj;
double nrmxl;
int jj, jp, pl, pu;
double tt, maxnrm;
int lp1, lup;
/* dqrdc uses householder transformations to compute the qr */
/* factorization of an n by p matrix x. column pivoting */
/* based on the 2-norms of the reduced columns may be */
/* performed at the users option. */
/* on entry */
/* x double precision(ldx,p), where ldx .ge. n. */
/* x contains the matrix whose decomposition is to be */
/* computed. */
/* ldx integer. */
/* ldx is the leading dimension of the array x. */
/* n integer. */
/* n is the number of rows of the matrix x. */
/* p integer. */
/* p is the number of columns of the matrix x. */
/* jpvt integer(p). */
/* jpvt contains integers that control the selection */
/* of the pivot columns. the k-th column x(k) of x */
/* is placed in one of three classes according to the */
/* value of jpvt(k). */
/* if jpvt(k) .gt. 0, then x(k) is an initial */
/* column. */
/* if jpvt(k) .eq. 0, then x(k) is a free column. */
/* if jpvt(k) .lt. 0, then x(k) is a final column. */
/* before the decomposition is computed, initial columns */
/* are moved to the beginning of the array x and final */
/* columns to the end. both initial and final columns */
/* are frozen in place during the computation and only */
/* free columns are moved. at the k-th stage of the */
/* reduction, if x(k) is occupied by a free column */
/* it is interchanged with the free column of largest */
/* reduced norm. jpvt is not referenced if */
/* job .eq. 0. */
/* work double precision(p). */
/* work is a work array. work is not referenced if */
/* job .eq. 0. */
/* job integer. */
/* job is an integer that initiates column pivoting. */
/* if job .eq. 0, no pivoting is done. */
/* if job .ne. 0, pivoting is done. */
/* on return */
/* x x contains in its upper triangle the upper */
/* triangular matrix r of the qr factorization. */
/* below its diagonal x contains information from */
/* which the orthogonal part of the decomposition */
/* can be recovered. note that if pivoting has */
/* been requested, the decomposition is not that */
/* of the original matrix x but that of x */
/* with its columns permuted as described by jpvt. */
/* qraux double precision(p). */
/* qraux contains further information required to recover */
/* the orthogonal part of the decomposition. */
/* jpvt jpvt(k) contains the index of the column of the */
/* original matrix that has been interchanged into */
/* the k-th column, if pivoting was requested. */
/* linpack. this version dated 08/14/78 . */
/* g.w. stewart, university of maryland, argonne national lab. */
/* dqrdc uses the following functions and subprograms. */
/* blas daxpy,ddot,dscal,dswap,dnrm2 */
/* fortran dabs,dmax1,min0,dsqrt */
/* Parameter adjustments */
x_dim1 = ldx;
x_offset = x_dim1 + 1;
x -= x_offset;
--qraux;
--jpvt;
--work;
/* Function Body */
pl = 1;
pu = 0;
if (job == 0) {
goto L60;
}
/* pivoting has been requested. rearrange the columns */
/* according to jpvt. */
i__1 = p;
for (j = 1; j <= i__1; ++j) {
swapj = jpvt[j] > 0;
negj = jpvt[j] < 0;
jpvt[j] = j;
if (negj) {
jpvt[j] = -j;
}
if (! swapj) {
goto L10;
}
if (j != pl) {
blas_dswap(n, &x[pl * x_dim1 + 1], 1, &x[j * x_dim1 + 1], 1);
}
jpvt[j] = jpvt[pl];
jpvt[pl] = j;
++pl;
L10:
/* L20: */
;
}
pu = p;
i__1 = p;
for (jj = 1; jj <= i__1; ++jj) {
j = p - jj + 1;
if (jpvt[j] >= 0) {
goto L40;
}
jpvt[j] = -jpvt[j];
if (j == pu) {
goto L30;
}
blas_dswap(n, &x[pu * x_dim1 + 1], 1, &x[j * x_dim1 + 1], 1);
jp = jpvt[pu];
jpvt[pu] = jpvt[j];
jpvt[j] = jp;
L30:
--pu;
L40:
/* L50: */
;
}
L60:
/* compute the norms of the free columns. */
if (pu < pl) {
goto L80;
}
i__1 = pu;
for (j = pl; j <= i__1; ++j) {
qraux[j] = blas_dnrm2(n, &x[j * x_dim1 + 1], 1);
work[j] = qraux[j];
/* L70: */
}
L80:
/* perform the householder reduction of x. */
lup = min(n,p);
i__1 = lup;
for (l = 1; l <= i__1; ++l) {
if (l < pl || l >= pu) {
goto L120;
}
/* locate the column of largest norm and bring it */
/* into the pivot position. */
maxnrm = 0.;
maxj = l;
i__2 = pu;
for (j = l; j <= i__2; ++j) {
if (qraux[j] <= maxnrm) {
goto L90;
}
maxnrm = qraux[j];
maxj = j;
L90:
/* L100: */
;
}
if (maxj == l) {
goto L110;
}
blas_dswap(n, &x[l * x_dim1 + 1], 1, &x[maxj * x_dim1 + 1], 1);
qraux[maxj] = qraux[l];
work[maxj] = work[l];
jp = jpvt[maxj];
jpvt[maxj] = jpvt[l];
jpvt[l] = jp;
L110:
L120:
qraux[l] = 0.;
if (l == n) {
goto L190;
}
/* compute the householder transformation for column l. */
i__2 = n - l + 1;
nrmxl = blas_dnrm2(i__2, &x[l + l * x_dim1], 1);
if (nrmxl == 0.) {
goto L180;
}
if (x[l + l * x_dim1] != 0.) {
nrmxl = d_sign(&nrmxl, &x[l + l * x_dim1]);
}
i__2 = n - l + 1;
d__1 = 1. / nrmxl;
blas_dscal(i__2, d__1, &x[l + l * x_dim1], 1);
x[l + l * x_dim1] += 1.;
/* apply the transformation to the remaining columns, */
/* updating the norms. */
lp1 = l + 1;
if (p < lp1) {
goto L170;
}
i__2 = p;
for (j = lp1; j <= i__2; ++j) {
i__3 = n - l + 1;
t = -blas_ddot(i__3, &x[l + l * x_dim1], 1, &x[l + j * x_dim1],
1) / x[l + l * x_dim1];
i__3 = n - l + 1;
blas_daxpy(i__3, t, &x[l + l * x_dim1], 1, &x[l + j * x_dim1], 1);
if (j < pl || j > pu) {
goto L150;
}
if (qraux[j] == 0.) {
goto L150;
}
/* Computing 2nd power */
d__2 = (d__1 = x[l + j * x_dim1], abs(d__1)) / qraux[j];
tt = 1. - d__2 * d__2;
tt = max(tt,0.);
t = tt;
/* Computing 2nd power */
d__1 = qraux[j] / work[j];
tt = tt * .05 * (d__1 * d__1) + 1.;
if (tt == 1.) {
goto L130;
}
qraux[j] *= sqrt(t);
goto L140;
L130:
i__3 = n - l;
qraux[j] = blas_dnrm2(i__3, &x[l + 1 + j * x_dim1], 1);
work[j] = qraux[j];
L140:
L150:
/* L160: */
;
}
L170:
/* save the transformation. */
qraux[l] = x[l + l * x_dim1];
x[l + l * x_dim1] = -nrmxl;
L180:
L190:
/* L200: */
;
}
return;
}
VOID linpack_dqrsl P13C(double *, x,
int, ldx,
int, n,
int, k,
double *, qraux,
double *, y,
double *, qy,
double *, qty,
double *, b,
double *, rsd,
double *, xb,
int, job,
int *, info)
{
/* System generated locals */
integer x_dim1, x_offset, i__1, i__2;
/* Local variables */
double temp;
int cqty;
int i, j;
double t;
int cb;
int jj;
int cr;
int ju, kp1;
int cxb, cqy;
/* dqrsl applies the output of dqrdc to compute coordinate */
/* transformations, projections, and least squares solutions. */
/* for k .le. min(n,p), let xk be the matrix */
/* xk = (x(jpvt(1)),x(jpvt(2)), ... ,x(jpvt(k))) */
/* formed from columnns jpvt(1), ... ,jpvt(k) of the original */
/* n x p matrix x that was input to dqrdc (if no pivoting was */
/* done, xk consists of the first k columns of x in their */
/* original order). dqrdc produces a factored orthogonal matrix q */
/* and an upper triangular matrix r such that */
/* xk = q * (r) */
/* (0) */
/* this information is contained in coded form in the arrays */
/* x and qraux. */
/* on entry */
/* x double precision(ldx,p). */
/* x contains the output of dqrdc. */
/* ldx integer. */
/* ldx is the leading dimension of the array x. */
/* n integer. */
/* n is the number of rows of the matrix xk. it must */
/* have the same value as n in dqrdc. */
/* k integer. */
/* k is the number of columns of the matrix xk. k */
/* must nnot be greater than min(n,p), where p is the */
/* same as in the calling sequence to dqrdc. */
/* qraux double precision(p). */
/* qraux contains the auxiliary output from dqrdc. */
/* y double precision(n) */
/* y contains an n-vector that is to be manipulated */
/* by dqrsl. */
/* job integer. */
/* job specifies what is to be computed. job has */
/* the decimal expansion abcde, with the following */
/* meaning. */
/* if a.ne.0, compute qy. */
/* if b,c,d, or e .ne. 0, compute qty. */
/* if c.ne.0, compute b. */
/* if d.ne.0, compute rsd. */
/* if e.ne.0, compute xb. */
/* note that a request to compute b, rsd, or xb */
/* automatically triggers the computation of qty, for */
/* which an array must be provided in the calling */
/* sequence. */
/* on return */
/* qy double precision(n). */
/* qy conntains q*y, if its computation has been */
/* requested. */
/* qty double precision(n). */
/* qty contains trans(q)*y, if its computation has */
/* been requested. here trans(q) is the */
/* transpose of the matrix q. */
/* b double precision(k) */
/* b contains the solution of the least squares problem */
/* minimize norm2(y - xk*b), */
/* if its computation has been requested. (note that */
/* if pivoting was requested in dqrdc, the j-th */
/* component of b will be associated with column jpvt(j) */
/* of the original matrix x that was input into dqrdc.) */
/* rsd double precision(n). */
/* rsd contains the least squares residual y - xk*b, */
/* if its computation has been requested. rsd is */
/* also the orthogonal projection of y onto the */
/* orthogonal complement of the column space of xk. */
/* xb double precision(n). */
/* xb contains the least squares approximation xk*b, */
/* if its computation has been requested. xb is also */
/* the orthogonal projection of y onto the column space */
/* of x. */
/* info integer. */
/* info is zero unless the computation of b has */
/* been requested and r is exactly singular. in */
/* this case, info is the index of the first zero */
/* diagonal element of r and b is left unaltered. */
/* the parameters qy, qty, b, rsd, and xb are not referenced */
/* if their computation is not requested and in this case */
/* can be replaced by dummy variables in the calling program. */
/* to save storage, the user may in some cases use the same */
/* array for different parameters in the calling sequence. a */
/* frequently occuring example is when one wishes to compute */
/* any of b, rsd, or xb and does not need y or qty. in this */
/* case one may identify y, qty, and one of b, rsd, or xb, while */
/* providing separate arrays for anything else that is to be */
/* computed. thus the calling sequence */
/* call dqrsl(x,ldx,n,k,qraux,y,dum,y,b,y,dum,110,info) */
/* will result in the computation of b and rsd, with rsd */
/* overwriting y. more generally, each item in the following */
/* list contains groups of permissible identifications for */
/* a single callinng sequence. */
/* 1. (y,qty,b) (rsd) (xb) (qy) */
/* 2. (y,qty,rsd) (b) (xb) (qy) */
/* 3. (y,qty,xb) (b) (rsd) (qy) */
/* 4. (y,qy) (qty,b) (rsd) (xb) */
/* 5. (y,qy) (qty,rsd) (b) (xb) */
/* 6. (y,qy) (qty,xb) (b) (rsd) */
/* in any group the value returned in the array allocated to */
/* the group corresponds to the last member of the group. */
/* linpack. this version dated 08/14/78 . */
/* g.w. stewart, university of maryland, argonne national lab. */
/* dqrsl uses the following functions and subprograms. */
/* blas daxpy,dcopy,ddot */
/* fortran dabs,min0,mod */
/* Parameter adjustments */
x_dim1 = ldx;
x_offset = x_dim1 + 1;
x -= x_offset;
--qraux;
--y;
--qy;
--qty;
--b;
--rsd;
--xb;
/* set info flag. */
*info = 0;
/* determine what is to be computed. */
cqy = job / 10000 != 0;
cqty = job % 10000 != 0;
cb = job % 1000 / 100 != 0;
cr = job % 100 / 10 != 0;
cxb = job % 10 != 0;
/* Computing MIN */
i__1 = k, i__2 = n - 1;
ju = min(i__1,i__2);
/* special action when n=1. */
if (ju != 0) {
goto L40;
}
if (cqy) {
qy[1] = y[1];
}
if (cqty) {
qty[1] = y[1];
}
if (cxb) {
xb[1] = y[1];
}
if (! cb) {
goto L30;
}
if (x[x_dim1 + 1] != 0.) {
goto L10;
}
*info = 1;
goto L20;
L10:
b[1] = y[1] / x[x_dim1 + 1];
L20:
L30:
if (cr) {
rsd[1] = 0.;
}
goto L250;
L40:
/* set up to compute qy or qty. */
if (cqy) {
blas_dcopy(n, &y[1], 1, &qy[1], 1);
}
if (cqty) {
blas_dcopy(n, &y[1], 1, &qty[1], 1);
}
if (! cqy) {
goto L70;
}
/* compute qy. */
i__1 = ju;
for (jj = 1; jj <= i__1; ++jj) {
j = ju - jj + 1;
if (qraux[j] == 0.) {
goto L50;
}
temp = x[j + j * x_dim1];
x[j + j * x_dim1] = qraux[j];
i__2 = n - j + 1;
t = -blas_ddot(i__2, &x[j + j * x_dim1], 1, &qy[j], 1) / x[j + j
* x_dim1];
i__2 = n - j + 1;
blas_daxpy(i__2, t, &x[j + j * x_dim1], 1, &qy[j], 1);
x[j + j * x_dim1] = temp;
L50:
/* L60: */
;
}
L70:
if (! cqty) {
goto L100;
}
/* compute trans(q)*y. */
i__1 = ju;
for (j = 1; j <= i__1; ++j) {
if (qraux[j] == 0.) {
goto L80;
}
temp = x[j + j * x_dim1];
x[j + j * x_dim1] = qraux[j];
i__2 = n - j + 1;
t = -blas_ddot(i__2, &x[j + j * x_dim1], 1, &qty[j], 1) / x[j +
j * x_dim1];
i__2 = n - j + 1;
blas_daxpy(i__2, t, &x[j + j * x_dim1], 1, &qty[j], 1);
x[j + j * x_dim1] = temp;
L80:
/* L90: */
;
}
L100:
/* set up to compute b, rsd, or xb. */
if (cb) {
blas_dcopy(k, &qty[1], 1, &b[1], 1);
}
kp1 = k + 1;
if (cxb) {
blas_dcopy(k, &qty[1], 1, &xb[1], 1);
}
if (cr && k < n) {
i__1 = n - k;
blas_dcopy(i__1, &qty[kp1], 1, &rsd[kp1], 1);
}
if (! cxb || kp1 > n) {
goto L120;
}
i__1 = n;
for (i = kp1; i <= i__1; ++i) {
xb[i] = 0.;
/* L110: */
}
L120:
if (! cr) {
goto L140;
}
i__1 = k;
for (i = 1; i <= i__1; ++i) {
rsd[i] = 0.;
/* L130: */
}
L140:
if (! cb) {
goto L190;
}
/* compute b. */
i__1 = k;
for (jj = 1; jj <= i__1; ++jj) {
j = k - jj + 1;
if (x[j + j * x_dim1] != 0.) {
goto L150;
}
*info = j;
/* ......exit */
goto L180;
L150:
b[j] /= x[j + j * x_dim1];
if (j == 1) {
goto L160;
}
t = -b[j];
i__2 = j - 1;
blas_daxpy(i__2, t, &x[j * x_dim1 + 1], 1, &b[1], 1);
L160:
/* L170: */
;
}
L180:
L190:
if (! cr && ! cxb) {
goto L240;
}
/* compute rsd or xb as required. */
i__1 = ju;
for (jj = 1; jj <= i__1; ++jj) {
j = ju - jj + 1;
if (qraux[j] == 0.) {
goto L220;
}
temp = x[j + j * x_dim1];
x[j + j * x_dim1] = qraux[j];
if (! cr) {
goto L200;
}
i__2 = n - j + 1;
t = -blas_ddot(i__2, &x[j + j * x_dim1], 1, &rsd[j], 1) / x[j +
j * x_dim1];
i__2 = n - j + 1;
blas_daxpy(i__2, t, &x[j + j * x_dim1], 1, &rsd[j], 1);
L200:
if (! cxb) {
goto L210;
}
i__2 = n - j + 1;
t = -blas_ddot(i__2, &x[j + j * x_dim1], 1, &xb[j], 1) / x[j + j
* x_dim1];
i__2 = n - j + 1;
blas_daxpy(i__2, t, &x[j + j * x_dim1], 1, &xb[j], 1);
L210:
x[j + j * x_dim1] = temp;
L220:
/* L230: */
;
}
L240:
L250:
return;
}
VOID linpack_zqrdc P8C(dcomplex *, x,
int, ldx,
int, n,
int, p,
dcomplex *, qraux,
int *, jpvt,
dcomplex *, work,
int, job)
{
/* System generated locals */
int x_dim1, x_offset, i__1, i__2, i__3, i__4;
double d__1, d__2, d__3, d__4;
dcomplex z__1, z__2, z__3;
static dcomplex c_b28 = {1.,0.};
/* Local variables */
int negj;
int maxj, j, l;
dcomplex t;
int swapj;
dcomplex nrmxl;
int jj, jp, pl, pu;
double tt;
double maxnrm;
int lp1, lup;
/* zqrdc uses householder transformations to compute the qr */
/* factorization of an n by p matrix x. column pivoting */
/* based on the 2-norms of the reduced columns may be */
/* performed at the users option. */
/* on entry */
/* x complex*16(ldx,p), where ldx .ge. n. */
/* x contains the matrix whose decomposition is to be */
/* computed. */
/* ldx integer. */
/* ldx is the leading dimension of the array x. */
/* n integer. */
/* n is the number of rows of the matrix x. */
/* p integer. */
/* p is the number of columns of the matrix x. */
/* jpvt integer(p). */
/* jpvt contains integers that control the selection */
/* of the pivot columns. the k-th column x(k) of x */
/* is placed in one of three classes according to the */
/* value of jpvt(k). */
/* if jpvt(k) .gt. 0, then x(k) is an initial */
/* column. */
/* if jpvt(k) .eq. 0, then x(k) is a free column. */
/* if jpvt(k) .lt. 0, then x(k) is a final column. */
/* before the decomposition is computed, initial columns */
/* are moved to the beginning of the array x and final */
/* columns to the end. both initial and final columns */
/* are frozen in place during the computation and only */
/* free columns are moved. at the k-th stage of the */
/* reduction, if x(k) is occupied by a free column */
/* it is interchanged with the free column of largest */
/* reduced norm. jpvt is not referenced if */
/* job .eq. 0. */
/* work complex*16(p). */
/* work is a work array. work is not referenced if */
/* job .eq. 0. */
/* job integer. */
/* job is an integer that initiates column pivoting. */
/* if job .eq. 0, no pivoting is done. */
/* if job .ne. 0, pivoting is done. */
/* on return */
/* x x contains in its upper triangle the upper */
/* triangular matrix r of the qr factorization. */
/* below its diagonal x contains information from */
/* which the unitary part of the decomposition */
/* can be recovered. note that if pivoting has */
/* been requested, the decomposition is not that */
/* of the original matrix x but that of x */
/* with its columns permuted as described by jpvt. */
/* qraux complex*16(p). */
/* qraux contains further information required to recover
*/
/* the unitary part of the decomposition. */
/* jpvt jpvt(k) contains the index of the column of the */
/* original matrix that has been interchanged into */
/* the k-th column, if pivoting was requested. */
/* linpack. this version dated 08/14/78 . */
/* g.w. stewart, university of maryland, argonne national lab. */
/* zqrdc uses the following functions and subprograms. */
/* blas zaxpy,zdotc,zscal,zswap,dznrm2 */
/* fortran dabs,dmax1,cdabs,dcmplx,cdsqrt,min0 */
/* internal variables */
/* Parameter adjustments */
x_dim1 = ldx;
x_offset = x_dim1 + 1;
x -= x_offset;
--qraux;
--jpvt;
--work;
/* Function Body */
pl = 1;
pu = 0;
if (job == 0) {
goto L60;
}
/* pivoting has been requested. rearrange the columns */
/* according to jpvt. */
i__1 = p;
for (j = 1; j <= i__1; ++j) {
swapj = jpvt[j] > 0;
negj = jpvt[j] < 0;
jpvt[j] = j;
if (negj) {
jpvt[j] = -j;
}
if (! swapj) {
goto L10;
}
if (j != pl) {
blas_zswap(n, &x[pl * x_dim1 + 1], 1, &x[j * x_dim1 + 1], 1);
}
jpvt[j] = jpvt[pl];
jpvt[pl] = j;
++pl;
L10:
/* L20: */
;
}
pu = p;
i__1 = p;
for (jj = 1; jj <= i__1; ++jj) {
j = p - jj + 1;
if (jpvt[j] >= 0) {
goto L40;
}
jpvt[j] = -jpvt[j];
if (j == pu) {
goto L30;
}
blas_zswap(n, &x[pu * x_dim1 + 1], 1, &x[j * x_dim1 + 1], 1);
jp = jpvt[pu];
jpvt[pu] = jpvt[j];
jpvt[j] = jp;
L30:
--pu;
L40:
/* L50: */
;
}
L60:
/* compute the norms of the free columns. */
if (pu < pl) {
goto L80;
}
i__1 = pu;
for (j = pl; j <= i__1; ++j) {
i__2 = j;
d__1 = blas_dznrm2(n, &x[j * x_dim1 + 1], 1);
z__1.r = d__1, z__1.i = 0.;
qraux[i__2].r = z__1.r, qraux[i__2].i = z__1.i;
i__2 = j;
i__3 = j;
work[i__2].r = qraux[i__3].r, work[i__2].i = qraux[i__3].i;
/* L70: */
}
L80:
/* perform the householder reduction of x. */
lup = min(n,p);
i__1 = lup;
for (l = 1; l <= i__1; ++l) {
if (l < pl || l >= pu) {
goto L120;
}
/* locate the column of largest norm and bring it */
/* into the pivot position. */
maxnrm = 0.;
maxj = l;
i__2 = pu;
for (j = l; j <= i__2; ++j) {
i__3 = j;
if (qraux[i__3].r <= maxnrm) {
goto L90;
}
i__3 = j;
maxnrm = qraux[i__3].r;
maxj = j;
L90:
/* L100: */
;
}
if (maxj == l) {
goto L110;
}
blas_zswap(n, &x[l * x_dim1 + 1], 1, &x[maxj * x_dim1 + 1], 1);
i__2 = maxj;
i__3 = l;
qraux[i__2].r = qraux[i__3].r, qraux[i__2].i = qraux[i__3].i;
i__2 = maxj;
i__3 = l;
work[i__2].r = work[i__3].r, work[i__2].i = work[i__3].i;
jp = jpvt[maxj];
jpvt[maxj] = jpvt[l];
jpvt[l] = jp;
L110:
L120:
i__2 = l;
qraux[i__2].r = 0., qraux[i__2].i = 0.;
if (l == n) {
goto L190;
}
/* compute the householder transformation for column l. */
i__2 = n - l + 1;
d__1 = blas_dznrm2(i__2, &x[l + l * x_dim1], 1);
z__1.r = d__1, z__1.i = 0.;
nrmxl.r = z__1.r, nrmxl.i = z__1.i;
z__1.r = nrmxl.r * 0. - nrmxl.i * -1., z__1.i = nrmxl.r * -1. +
nrmxl.i * 0.;
if ((d__1 = nrmxl.r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) == 0.) {
goto L180;
}
i__2 = l + l * x_dim1;
i__3 = l + l * x_dim1;
z__1.r = x[i__3].r * 0. - x[i__3].i * -1., z__1.i = x[i__3].r * -1. +
x[i__3].i * 0.;
if ((d__1 = x[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) != 0.)
{
d__3 = z_abs(&nrmxl);
i__4 = l + l * x_dim1;
d__4 = z_abs(&x[l + l * x_dim1]);
z__3.r = x[i__4].r / d__4, z__3.i = x[i__4].i / d__4;
z__2.r = d__3 * z__3.r, z__2.i = d__3 * z__3.i;
nrmxl.r = z__2.r, nrmxl.i = z__2.i;
}
i__2 = n - l + 1;
z_div(&z__1, &c_b28, &nrmxl);
blas_zscal(i__2, &z__1, &x[l + l * x_dim1], 1);
i__2 = l + l * x_dim1;
i__3 = l + l * x_dim1;
z__1.r = x[i__3].r + 1., z__1.i = x[i__3].i + 0.;
x[i__2].r = z__1.r, x[i__2].i = z__1.i;
/* apply the transformation to the remaining columns, */
/* updating the norms. */
lp1 = l + 1;
if (p < lp1) {
goto L170;
}
i__2 = p;
for (j = lp1; j <= i__2; ++j) {
i__3 = n - l + 1;
blas_zdotc(&z__3, i__3, &x[l + l * x_dim1], 1, &x[l + j * x_dim1], 1);
z__2.r = -z__3.r, z__2.i = -z__3.i;
z_div(&z__1, &z__2, &x[l + l * x_dim1]);
t.r = z__1.r, t.i = z__1.i;
i__3 = n - l + 1;
blas_zaxpy(i__3, &t, &x[l + l * x_dim1], 1, &x[l + j * x_dim1], 1);
if (j < pl || j > pu) {
goto L150;
}
i__3 = j;
i__4 = j;
z__1.r = qraux[i__4].r * 0. - qraux[i__4].i * -1., z__1.i = qraux[
i__4].r * -1. + qraux[i__4].i * 0.;
if ((d__1 = qraux[i__3].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2))
== 0.) {
goto L150;
}
i__3 = j;
/* Computing 2nd power */
d__1 = z_abs(&x[l + j * x_dim1]) / qraux[i__3].r;
tt = 1. - d__1 * d__1;
tt = max(tt,0.);
z__1.r = tt, z__1.i = 0.;
t.r = z__1.r, t.i = z__1.i;
i__3 = j;
i__4 = j;
/* Computing 2nd power */
d__1 = qraux[i__3].r / work[i__4].r;
tt = tt * .05 * (d__1 * d__1) + 1.;
if (tt == 1.) {
goto L130;
}
i__3 = j;
i__4 = j;
z_sqrt(&z__2, &t);
z__1.r = qraux[i__4].r * z__2.r - qraux[i__4].i * z__2.i, z__1.i =
qraux[i__4].r * z__2.i + qraux[i__4].i * z__2.r;
qraux[i__3].r = z__1.r, qraux[i__3].i = z__1.i;
goto L140;
L130:
i__3 = j;
i__4 = n - l;
d__1 = blas_dznrm2(i__4, &x[l + 1 + j * x_dim1], 1);
z__1.r = d__1, z__1.i = 0.;
qraux[i__3].r = z__1.r, qraux[i__3].i = z__1.i;
i__3 = j;
i__4 = j;
work[i__3].r = qraux[i__4].r, work[i__3].i = qraux[i__4].i;
L140:
L150:
/* L160: */
;
}
L170:
/* save the transformation. */
i__2 = l;
i__3 = l + l * x_dim1;
qraux[i__2].r = x[i__3].r, qraux[i__2].i = x[i__3].i;
i__2 = l + l * x_dim1;
z__1.r = -nrmxl.r, z__1.i = -nrmxl.i;
x[i__2].r = z__1.r, x[i__2].i = z__1.i;
L180:
L190:
/* L200: */
;
}
return;
}
VOID linpack_zqrsl P13C(dcomplex *, x,
int, ldx,
int, n,
int, k,
dcomplex *, qraux,
dcomplex *, y,
dcomplex *, qy,
dcomplex *, qty,
dcomplex *, b,
dcomplex *, rsd,
dcomplex *, xb,
int, job,
int *, info)
{
/* System generated locals */
int x_dim1, x_offset, i__1, i__2, i__3;
double d__1, d__2;
dcomplex z__1, z__2, z__3;
/* Local variables */
dcomplex temp;
int cqty;
int i, j;
dcomplex t;
int cb;
int jj;
int cr;
int ju, kp1;
int cxb, cqy;
/* zqrsl applies the output of zqrdc to compute coordinate */
/* transformations, projections, and least squares solutions. */
/* for k .le. min(n,p), let xk be the matrix */
/* xk = (x(jpvt(1)),x(jpvt(2)), ... ,x(jpvt(k))) */
/* formed from columnns jpvt(1), ... ,jpvt(k) of the original */
/* n x p matrix x that was input to zqrdc (if no pivoting was */
/* done, xk consists of the first k columns of x in their */
/* original order). zqrdc produces a factored unitary matrix q */
/* and an upper triangular matrix r such that */
/* xk = q * (r) */
/* (0) */
/* this information is contained in coded form in the arrays */
/* x and qraux. */
/* on entry */
/* x complex*16(ldx,p). */
/* x contains the output of zqrdc. */
/* ldx integer. */
/* ldx is the leading dimension of the array x. */
/* n integer. */
/* n is the number of rows of the matrix xk. it must */
/* have the same value as n in zqrdc. */
/* k integer. */
/* k is the number of columns of the matrix xk. k */
/* must nnot be greater than min(n,p), where p is the */
/* same as in the calling sequence to zqrdc. */
/* qraux complex*16(p). */
/* qraux contains the auxiliary output from zqrdc. */
/* y complex*16(n) */
/* y contains an n-vector that is to be manipulated */
/* by zqrsl. */
/* job integer. */
/* job specifies what is to be computed. job has */
/* the decimal expansion abcde, with the following */
/* meaning. */
/* if a.ne.0, compute qy. */
/* if b,c,d, or e .ne. 0, compute qty. */
/* if c.ne.0, compute b. */
/* if d.ne.0, compute rsd. */
/* if e.ne.0, compute xb. */
/* note that a request to compute b, rsd, or xb */
/* automatically triggers the computation of qty, for */
/* which an array must be provided in the calling */
/* sequence. */
/* on return */
/* qy complex*16(n). */
/* qy conntains q*y, if its computation has been */
/* requested. */
/* qty complex*16(n). */
/* qty contains ctrans(q)*y, if its computation has */
/* been requested. here ctrans(q) is the conjugate */
/* transpose of the matrix q. */
/* b complex*16(k) */
/* b contains the solution of the least squares problem */
/* minimize norm2(y - xk*b), */
/* if its computation has been requested. (note that */
/* if pivoting was requested in zqrdc, the j-th */
/* component of b will be associated with column jpvt(j) */
/* of the original matrix x that was input into zqrdc.) */
/* rsd complex*16(n). */
/* rsd contains the least squares residual y - xk*b, */
/* if its computation has been requested. rsd is */
/* also the orthogonal projection of y onto the */
/* orthogonal complement of the column space of xk. */
/* xb complex*16(n). */
/* xb contains the least squares approximation xk*b, */
/* if its computation has been requested. xb is also */
/* the orthogonal projection of y onto the column space */
/* of x. */
/* info integer. */
/* info is zero unless the computation of b has */
/* been requested and r is exactly singular. in */
/* this case, info is the index of the first zero */
/* diagonal element of r and b is left unaltered. */
/* the parameters qy, qty, b, rsd, and xb are not referenced */
/* if their computation is not requested and in this case */
/* can be replaced by dummy variables in the calling program. */
/* to save storage, the user may in some cases use the same */
/* array for different parameters in the calling sequence. a */
/* frequently occuring example is when one wishes to compute */
/* any of b, rsd, or xb and does not need y or qty. in this */
/* case one may identify y, qty, and one of b, rsd, or xb, while */
/* providing separate arrays for anything else that is to be */
/* computed. thus the calling sequence */
/* call zqrsl(x,ldx,n,k,qraux,y,dum,y,b,y,dum,110,info) */
/* will result in the computation of b and rsd, with rsd */
/* overwriting y. more generally, each item in the following */
/* list contains groups of permissible identifications for */
/* a single callinng sequence. */
/* 1. (y,qty,b) (rsd) (xb) (qy) */
/* 2. (y,qty,rsd) (b) (xb) (qy) */
/* 3. (y,qty,xb) (b) (rsd) (qy) */
/* 4. (y,qy) (qty,b) (rsd) (xb) */
/* 5. (y,qy) (qty,rsd) (b) (xb) */
/* 6. (y,qy) (qty,xb) (b) (rsd) */
/* in any group the value returned in the array allocated to */
/* the group corresponds to the last member of the group. */
/* linpack. this version dated 08/14/78 . */
/* g.w. stewart, university of maryland, argonne national lab. */
/* zqrsl uses the following functions and subprograms. */
/* blas zaxpy,zcopy,zdotc */
/* fortran dabs,min0,mod */
/* internal variables */
/* set info flag. */
/* Parameter adjustments */
x_dim1 = ldx;
x_offset = x_dim1 + 1;
x -= x_offset;
--qraux;
--y;
--qy;
--qty;
--b;
--rsd;
--xb;
/* Function Body */
*info = 0;
/* determine what is to be computed. */
cqy = job / 10000 != 0;
cqty = job % 10000 != 0;
cb = job % 1000 / 100 != 0;
cr = job % 100 / 10 != 0;
cxb = job % 10 != 0;
/* Computing MIN */
i__1 = k, i__2 = n - 1;
ju = min(i__1,i__2);
/* special action when n=1. */
if (ju != 0) {
goto L40;
}
if (cqy) {
qy[1].r = y[1].r, qy[1].i = y[1].i;
}
if (cqty) {
qty[1].r = y[1].r, qty[1].i = y[1].i;
}
if (cxb) {
xb[1].r = y[1].r, xb[1].i = y[1].i;
}
if (! cb) {
goto L30;
}
i__1 = x_dim1 + 1;
i__2 = x_dim1 + 1;
z__1.r = x[i__2].r * 0. - x[i__2].i * -1., z__1.i = x[i__2].r * -1. + x[
i__2].i * 0.;
if ((d__1 = x[i__1].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) != 0.) {
goto L10;
}
*info = 1;
goto L20;
L10:
z_div(&z__1, &y[1], &x[x_dim1 + 1]);
b[1].r = z__1.r, b[1].i = z__1.i;
L20:
L30:
if (cr) {
rsd[1].r = 0., rsd[1].i = 0.;
}
goto L250;
L40:
/* set up to compute qy or qty. */
if (cqy) {
blas_zcopy(n, &y[1], 1, &qy[1], 1);
}
if (cqty) {
blas_zcopy(n, &y[1], 1, &qty[1], 1);
}
if (! cqy) {
goto L70;
}
/* compute qy. */
i__1 = ju;
for (jj = 1; jj <= i__1; ++jj) {
j = ju - jj + 1;
i__2 = j;
i__3 = j;
z__1.r = qraux[i__3].r * 0. - qraux[i__3].i * -1., z__1.i = qraux[
i__3].r * -1. + qraux[i__3].i * 0.;
if ((d__1 = qraux[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) ==
0.) {
goto L50;
}
i__2 = j + j * x_dim1;
temp.r = x[i__2].r, temp.i = x[i__2].i;
i__2 = j + j * x_dim1;
i__3 = j;
x[i__2].r = qraux[i__3].r, x[i__2].i = qraux[i__3].i;
i__2 = n - j + 1;
blas_zdotc(&z__3, i__2, &x[j + j * x_dim1], 1, &qy[j], 1);
z__2.r = -z__3.r, z__2.i = -z__3.i;
z_div(&z__1, &z__2, &x[j + j * x_dim1]);
t.r = z__1.r, t.i = z__1.i;
i__2 = n - j + 1;
blas_zaxpy(i__2, &t, &x[j + j * x_dim1], 1, &qy[j], 1);
i__2 = j + j * x_dim1;
x[i__2].r = temp.r, x[i__2].i = temp.i;
L50:
/* L60: */
;
}
L70:
if (! cqty) {
goto L100;
}
/* compute ctrans(q)*y. */
i__1 = ju;
for (j = 1; j <= i__1; ++j) {
i__2 = j;
i__3 = j;
z__1.r = qraux[i__3].r * 0. - qraux[i__3].i * -1., z__1.i = qraux[
i__3].r * -1. + qraux[i__3].i * 0.;
if ((d__1 = qraux[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) ==
0.) {
goto L80;
}
i__2 = j + j * x_dim1;
temp.r = x[i__2].r, temp.i = x[i__2].i;
i__2 = j + j * x_dim1;
i__3 = j;
x[i__2].r = qraux[i__3].r, x[i__2].i = qraux[i__3].i;
i__2 = n - j + 1;
blas_zdotc(&z__3, i__2, &x[j + j * x_dim1], 1, &qty[j], 1);
z__2.r = -z__3.r, z__2.i = -z__3.i;
z_div(&z__1, &z__2, &x[j + j * x_dim1]);
t.r = z__1.r, t.i = z__1.i;
i__2 = n - j + 1;
blas_zaxpy(i__2, &t, &x[j + j * x_dim1], 1, &qty[j], 1);
i__2 = j + j * x_dim1;
x[i__2].r = temp.r, x[i__2].i = temp.i;
L80:
/* L90: */
;
}
L100:
/* set up to compute b, rsd, or xb. */
if (cb) {
blas_zcopy(k, &qty[1], 1, &b[1], 1);
}
kp1 = k + 1;
if (cxb) {
blas_zcopy(k, &qty[1], 1, &xb[1], 1);
}
if (cr && k < n) {
i__1 = n - k;
blas_zcopy(i__1, &qty[kp1], 1, &rsd[kp1], 1);
}
if (! cxb || kp1 > n) {
goto L120;
}
i__1 = n;
for (i = kp1; i <= i__1; ++i) {
i__2 = i;
xb[i__2].r = 0., xb[i__2].i = 0.;
/* L110: */
}
L120:
if (! cr) {
goto L140;
}
i__1 = k;
for (i = 1; i <= i__1; ++i) {
i__2 = i;
rsd[i__2].r = 0., rsd[i__2].i = 0.;
/* L130: */
}
L140:
if (! cb) {
goto L190;
}
/* compute b. */
i__1 = k;
for (jj = 1; jj <= i__1; ++jj) {
j = k - jj + 1;
i__2 = j + j * x_dim1;
i__3 = j + j * x_dim1;
z__1.r = x[i__3].r * 0. - x[i__3].i * -1., z__1.i = x[i__3].r * -1. +
x[i__3].i * 0.;
if ((d__1 = x[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) != 0.)
{
goto L150;
}
*info = j;
/* ......exit */
goto L180;
L150:
i__2 = j;
z_div(&z__1, &b[j], &x[j + j * x_dim1]);
b[i__2].r = z__1.r, b[i__2].i = z__1.i;
if (j == 1) {
goto L160;
}
i__2 = j;
z__1.r = -b[i__2].r, z__1.i = -b[i__2].i;
t.r = z__1.r, t.i = z__1.i;
i__2 = j - 1;
blas_zaxpy(i__2, &t, &x[j * x_dim1 + 1], 1, &b[1], 1);
L160:
/* L170: */
;
}
L180:
L190:
if (! cr && ! cxb) {
goto L240;
}
/* compute rsd or xb as required. */
i__1 = ju;
for (jj = 1; jj <= i__1; ++jj) {
j = ju - jj + 1;
i__2 = j;
i__3 = j;
z__1.r = qraux[i__3].r * 0. - qraux[i__3].i * -1., z__1.i = qraux[
i__3].r * -1. + qraux[i__3].i * 0.;
if ((d__1 = qraux[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) ==
0.) {
goto L220;
}
i__2 = j + j * x_dim1;
temp.r = x[i__2].r, temp.i = x[i__2].i;
i__2 = j + j * x_dim1;
i__3 = j;
x[i__2].r = qraux[i__3].r, x[i__2].i = qraux[i__3].i;
if (! cr) {
goto L200;
}
i__2 = n - j + 1;
blas_zdotc(&z__3, i__2, &x[j + j * x_dim1], 1, &rsd[j], 1);
z__2.r = -z__3.r, z__2.i = -z__3.i;
z_div(&z__1, &z__2, &x[j + j * x_dim1]);
t.r = z__1.r, t.i = z__1.i;
i__2 = n - j + 1;
blas_zaxpy(i__2, &t, &x[j + j * x_dim1], 1, &rsd[j], 1);
L200:
if (! cxb) {
goto L210;
}
i__2 = n - j + 1;
blas_zdotc(&z__3, i__2, &x[j + j * x_dim1], 1, &xb[j], 1);
z__2.r = -z__3.r, z__2.i = -z__3.i;
z_div(&z__1, &z__2, &x[j + j * x_dim1]);
t.r = z__1.r, t.i = z__1.i;
i__2 = n - j + 1;
blas_zaxpy(i__2, &t, &x[j + j * x_dim1], 1, &xb[j], 1);
L210:
i__2 = j + j * x_dim1;
x[i__2].r = temp.r, x[i__2].i = temp.i;
L220:
/* L230: */
;
}
L240:
L250:
return;
}
syntax highlighted by Code2HTML, v. 0.9.1