1 /* QR Decomposition from LINPACK */
5 VOID linpack_dqrdc
P8C(double *, x
,
14 /* System generated locals */
15 int x_dim1
, x_offset
, i__1
, i__2
, i__3
;
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. */
36 /* x double precision(ldx,p), where ldx .ge. n. */
37 /* x contains the matrix whose decomposition is to be */
41 /* ldx is the leading dimension of the array x. */
44 /* n is the number of rows of the matrix x. */
47 /* p is the number of columns of the matrix x. */
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). */
55 /* if jpvt(k) .gt. 0, then x(k) is an initial */
58 /* if jpvt(k) .eq. 0, then x(k) is a free column. */
60 /* if jpvt(k) .lt. 0, then x(k) is a final column. */
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 */
72 /* work double precision(p). */
73 /* work is a work array. work is not referenced if */
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. */
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. */
92 /* qraux double precision(p). */
93 /* qraux contains further information required to recover */
94 /* the orthogonal part of the decomposition. */
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. */
100 /* linpack. this version dated 08/14/78 . */
101 /* g.w. stewart, university of maryland, argonne national lab. */
103 /* dqrdc uses the following functions and subprograms. */
105 /* blas daxpy,ddot,dscal,dswap,dnrm2 */
106 /* fortran dabs,dmax1,min0,dsqrt */
108 /* Parameter adjustments */
110 x_offset
= x_dim1
+ 1;
123 /* pivoting has been requested. rearrange the columns */
124 /* according to jpvt. */
127 for (j
= 1; j
<= i__1
; ++j
) {
138 blas_dswap(n
, &x
[pl
* x_dim1
+ 1], 1, &x
[j
* x_dim1
+ 1], 1);
149 for (jj
= 1; jj
<= i__1
; ++jj
) {
158 blas_dswap(n
, &x
[pu
* x_dim1
+ 1], 1, &x
[j
* x_dim1
+ 1], 1);
170 /* compute the norms of the free columns. */
176 for (j
= pl
; j
<= i__1
; ++j
) {
177 qraux
[j
] = blas_dnrm2(n
, &x
[j
* x_dim1
+ 1], 1);
183 /* perform the householder reduction of x. */
187 for (l
= 1; l
<= i__1
; ++l
) {
188 if (l
< pl
|| l
>= pu
) {
192 /* locate the column of largest norm and bring it */
193 /* into the pivot position. */
198 for (j
= l
; j
<= i__2
; ++j
) {
199 if (qraux
[j
] <= maxnrm
) {
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
];
215 jpvt
[maxj
] = jpvt
[l
];
224 /* compute the householder transformation for column l. */
227 nrmxl
= blas_dnrm2(i__2
, &x
[l
+ l
* x_dim1
], 1);
231 if (x
[l
+ l
* x_dim1
] != 0.) {
232 nrmxl
= d_sign(&nrmxl
, &x
[l
+ l
* x_dim1
]);
236 blas_dscal(i__2
, d__1
, &x
[l
+ l
* x_dim1
], 1);
237 x
[l
+ l
* x_dim1
] += 1.;
239 /* apply the transformation to the remaining columns, */
240 /* updating the norms. */
247 for (j
= lp1
; j
<= i__2
; ++j
) {
249 t
= -blas_ddot(i__3
, &x
[l
+ l
* x_dim1
], 1, &x
[l
+ j
* x_dim1
],
250 1) / x
[l
+ l
* x_dim1
];
252 blas_daxpy(i__3
, t
, &x
[l
+ l
* x_dim1
], 1, &x
[l
+ j
* x_dim1
], 1);
253 if (j
< pl
|| j
> pu
) {
256 if (qraux
[j
] == 0.) {
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
;
264 /* Computing 2nd power */
265 d__1
= qraux
[j
] / work
[j
];
266 tt
= tt
* .05 * (d__1
* d__1
) + 1.;
274 qraux
[j
] = blas_dnrm2(i__3
, &x
[l
+ 1 + j
* x_dim1
], 1);
283 /* save the transformation. */
285 qraux
[l
] = x
[l
+ l
* x_dim1
];
286 x
[l
+ l
* x_dim1
] = -nrmxl
;
295 VOID linpack_dqrsl
P13C(double *, x
,
309 /* System generated locals */
310 integer x_dim1
, x_offset
, i__1
, i__2
;
312 /* Local variables */
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 */
327 /* xk = (x(jpvt(1)),x(jpvt(2)), ... ,x(jpvt(k))) */
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 */
338 /* this information is contained in coded form in the arrays */
343 /* x double precision(ldx,p). */
344 /* x contains the output of dqrdc. */
347 /* ldx is the leading dimension of the array x. */
350 /* n is the number of rows of the matrix xk. it must */
351 /* have the same value as n in dqrdc. */
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. */
358 /* qraux double precision(p). */
359 /* qraux contains the auxiliary output from dqrdc. */
361 /* y double precision(n) */
362 /* y contains an n-vector that is to be manipulated */
366 /* job specifies what is to be computed. job has */
367 /* the decimal expansion abcde, with the following */
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. */
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 */
383 /* qy double precision(n). */
384 /* qy conntains q*y, if its computation has been */
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. */
392 /* b double precision(k) */
393 /* b contains the solution of the least squares problem */
395 /* minimize norm2(y - xk*b), */
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.) */
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. */
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 */
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. */
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 */
431 /* call dqrsl(x,ldx,n,k,qraux,y,dum,y,b,y,dum,110,info) */
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. */
438 /* 1. (y,qty,b) (rsd) (xb) (qy) */
440 /* 2. (y,qty,rsd) (b) (xb) (qy) */
442 /* 3. (y,qty,xb) (b) (rsd) (qy) */
444 /* 4. (y,qy) (qty,b) (rsd) (xb) */
446 /* 5. (y,qy) (qty,rsd) (b) (xb) */
448 /* 6. (y,qy) (qty,xb) (b) (rsd) */
450 /* in any group the value returned in the array allocated to */
451 /* the group corresponds to the last member of the group. */
453 /* linpack. this version dated 08/14/78 . */
454 /* g.w. stewart, university of maryland, argonne national lab. */
456 /* dqrsl uses the following functions and subprograms. */
458 /* blas daxpy,dcopy,ddot */
459 /* fortran dabs,min0,mod */
461 /* Parameter adjustments */
463 x_offset
= x_dim1
+ 1;
477 /* determine what is to be computed. */
479 cqy
= job
/ 10000 != 0;
480 cqty
= job
% 10000 != 0;
481 cb
= job
% 1000 / 100 != 0;
482 cr
= job
% 100 / 10 != 0;
485 i__1
= k
, i__2
= n
- 1;
488 /* special action when n=1. */
505 if (x
[x_dim1
+ 1] != 0.) {
511 b
[1] = y
[1] / x
[x_dim1
+ 1];
520 /* set up to compute qy or qty. */
523 blas_dcopy(n
, &y
[1], 1, &qy
[1], 1);
526 blas_dcopy(n
, &y
[1], 1, &qty
[1], 1);
535 for (jj
= 1; jj
<= i__1
; ++jj
) {
537 if (qraux
[j
] == 0.) {
540 temp
= x
[j
+ j
* x_dim1
];
541 x
[j
+ j
* x_dim1
] = qraux
[j
];
543 t
= -blas_ddot(i__2
, &x
[j
+ j
* x_dim1
], 1, &qy
[j
], 1) / x
[j
+ j
546 blas_daxpy(i__2
, t
, &x
[j
+ j
* x_dim1
], 1, &qy
[j
], 1);
547 x
[j
+ j
* x_dim1
] = temp
;
557 /* compute trans(q)*y. */
560 for (j
= 1; j
<= i__1
; ++j
) {
561 if (qraux
[j
] == 0.) {
564 temp
= x
[j
+ j
* x_dim1
];
565 x
[j
+ j
* x_dim1
] = qraux
[j
];
567 t
= -blas_ddot(i__2
, &x
[j
+ j
* x_dim1
], 1, &qty
[j
], 1) / x
[j
+
570 blas_daxpy(i__2
, t
, &x
[j
+ j
* x_dim1
], 1, &qty
[j
], 1);
571 x
[j
+ j
* x_dim1
] = temp
;
578 /* set up to compute b, rsd, or xb. */
581 blas_dcopy(k
, &qty
[1], 1, &b
[1], 1);
585 blas_dcopy(k
, &qty
[1], 1, &xb
[1], 1);
589 blas_dcopy(i__1
, &qty
[kp1
], 1, &rsd
[kp1
], 1);
591 if (! cxb
|| kp1
> n
) {
595 for (i
= kp1
; i
<= i__1
; ++i
) {
604 for (i
= 1; i
<= i__1
; ++i
) {
616 for (jj
= 1; jj
<= i__1
; ++jj
) {
618 if (x
[j
+ j
* x_dim1
] != 0.) {
625 b
[j
] /= x
[j
+ j
* x_dim1
];
631 blas_daxpy(i__2
, t
, &x
[j
* x_dim1
+ 1], 1, &b
[1], 1);
642 /* compute rsd or xb as required. */
645 for (jj
= 1; jj
<= i__1
; ++jj
) {
647 if (qraux
[j
] == 0.) {
650 temp
= x
[j
+ j
* x_dim1
];
651 x
[j
+ j
* x_dim1
] = qraux
[j
];
656 t
= -blas_ddot(i__2
, &x
[j
+ j
* x_dim1
], 1, &rsd
[j
], 1) / x
[j
+
659 blas_daxpy(i__2
, t
, &x
[j
+ j
* x_dim1
], 1, &rsd
[j
], 1);
665 t
= -blas_ddot(i__2
, &x
[j
+ j
* x_dim1
], 1, &xb
[j
], 1) / x
[j
+ j
668 blas_daxpy(i__2
, t
, &x
[j
+ j
* x_dim1
], 1, &xb
[j
], 1);
670 x
[j
+ j
* x_dim1
] = temp
;
680 VOID linpack_zqrdc
P8C(dcomplex
*, x
,
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.};
695 /* Local variables */
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. */
714 /* x complex*16(ldx,p), where ldx .ge. n. */
715 /* x contains the matrix whose decomposition is to be */
719 /* ldx is the leading dimension of the array x. */
722 /* n is the number of rows of the matrix x. */
725 /* p is the number of columns of the matrix x. */
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). */
733 /* if jpvt(k) .gt. 0, then x(k) is an initial */
736 /* if jpvt(k) .eq. 0, then x(k) is a free column. */
738 /* if jpvt(k) .lt. 0, then x(k) is a final column. */
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 */
750 /* work complex*16(p). */
751 /* work is a work array. work is not referenced if */
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. */
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. */
770 /* qraux complex*16(p). */
771 /* qraux contains further information required to recover
773 /* the unitary part of the decomposition. */
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. */
779 /* linpack. this version dated 08/14/78 . */
780 /* g.w. stewart, university of maryland, argonne national lab. */
782 /* zqrdc uses the following functions and subprograms. */
784 /* blas zaxpy,zdotc,zscal,zswap,dznrm2 */
785 /* fortran dabs,dmax1,cdabs,dcmplx,cdsqrt,min0 */
787 /* internal variables */
791 /* Parameter adjustments */
793 x_offset
= x_dim1
+ 1;
806 /* pivoting has been requested. rearrange the columns */
807 /* according to jpvt. */
810 for (j
= 1; j
<= i__1
; ++j
) {
821 blas_zswap(n
, &x
[pl
* x_dim1
+ 1], 1, &x
[j
* x_dim1
+ 1], 1);
832 for (jj
= 1; jj
<= i__1
; ++jj
) {
841 blas_zswap(n
, &x
[pu
* x_dim1
+ 1], 1, &x
[j
* x_dim1
+ 1], 1);
853 /* compute the norms of the free columns. */
859 for (j
= pl
; j
<= i__1
; ++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
;
866 work
[i__2
].r
= qraux
[i__3
].r
, work
[i__2
].i
= qraux
[i__3
].i
;
871 /* perform the householder reduction of x. */
875 for (l
= 1; l
<= i__1
; ++l
) {
876 if (l
< pl
|| l
>= pu
) {
880 /* locate the column of largest norm and bring it */
881 /* into the pivot position. */
886 for (j
= l
; j
<= i__2
; ++j
) {
888 if (qraux
[i__3
].r
<= maxnrm
) {
892 maxnrm
= qraux
[i__3
].r
;
901 blas_zswap(n
, &x
[l
* x_dim1
+ 1], 1, &x
[maxj
* x_dim1
+ 1], 1);
904 qraux
[i__2
].r
= qraux
[i__3
].r
, qraux
[i__2
].i
= qraux
[i__3
].i
;
907 work
[i__2
].r
= work
[i__3
].r
, work
[i__2
].i
= work
[i__3
].i
;
909 jpvt
[maxj
] = jpvt
[l
];
914 qraux
[i__2
].r
= 0., qraux
[i__2
].i
= 0.;
919 /* compute the householder transformation for column l. */
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. +
927 if ((d__1
= nrmxl
.r
, abs(d__1
)) + (d__2
= z__1
.r
, abs(d__2
)) == 0.) {
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. +
934 if ((d__1
= x
[i__2
].r
, abs(d__1
)) + (d__2
= z__1
.r
, abs(d__2
)) != 0.)
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
;
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
;
951 /* apply the transformation to the remaining columns, */
952 /* updating the norms. */
959 for (j
= lp1
; j
<= i__2
; ++j
) {
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
;
966 blas_zaxpy(i__3
, &t
, &x
[l
+ l
* x_dim1
], 1, &x
[l
+ j
* x_dim1
], 1);
967 if (j
< pl
|| j
> pu
) {
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
))
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
;
983 z__1
.r
= tt
, z__1
.i
= 0.;
984 t
.r
= z__1
.r
, t
.i
= z__1
.i
;
987 /* Computing 2nd power */
988 d__1
= qraux
[i__3
].r
/ work
[i__4
].r
;
989 tt
= tt
* .05 * (d__1
* d__1
) + 1.;
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
;
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
;
1008 work
[i__3
].r
= qraux
[i__4
].r
, work
[i__3
].i
= qraux
[i__4
].i
;
1016 /* save the transformation. */
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
;
1032 VOID linpack_zqrsl
P13C(dcomplex
*, x
,
1046 /* System generated locals */
1047 int x_dim1
, x_offset
, i__1
, i__2
, i__3
;
1049 dcomplex z__1
, z__2
, z__3
;
1051 /* Local variables */
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 */
1067 /* xk = (x(jpvt(1)),x(jpvt(2)), ... ,x(jpvt(k))) */
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 */
1078 /* this information is contained in coded form in the arrays */
1083 /* x complex*16(ldx,p). */
1084 /* x contains the output of zqrdc. */
1087 /* ldx is the leading dimension of the array x. */
1090 /* n is the number of rows of the matrix xk. it must */
1091 /* have the same value as n in zqrdc. */
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. */
1098 /* qraux complex*16(p). */
1099 /* qraux contains the auxiliary output from zqrdc. */
1101 /* y complex*16(n) */
1102 /* y contains an n-vector that is to be manipulated */
1106 /* job specifies what is to be computed. job has */
1107 /* the decimal expansion abcde, with the following */
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. */
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 */
1123 /* qy complex*16(n). */
1124 /* qy conntains q*y, if its computation has been */
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. */
1132 /* b complex*16(k) */
1133 /* b contains the solution of the least squares problem */
1135 /* minimize norm2(y - xk*b), */
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.) */
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. */
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 */
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. */
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 */
1171 /* call zqrsl(x,ldx,n,k,qraux,y,dum,y,b,y,dum,110,info) */
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. */
1178 /* 1. (y,qty,b) (rsd) (xb) (qy) */
1180 /* 2. (y,qty,rsd) (b) (xb) (qy) */
1182 /* 3. (y,qty,xb) (b) (rsd) (qy) */
1184 /* 4. (y,qy) (qty,b) (rsd) (xb) */
1186 /* 5. (y,qy) (qty,rsd) (b) (xb) */
1188 /* 6. (y,qy) (qty,xb) (b) (rsd) */
1190 /* in any group the value returned in the array allocated to */
1191 /* the group corresponds to the last member of the group. */
1193 /* linpack. this version dated 08/14/78 . */
1194 /* g.w. stewart, university of maryland, argonne national lab. */
1196 /* zqrsl uses the following functions and subprograms. */
1198 /* blas zaxpy,zcopy,zdotc */
1199 /* fortran dabs,min0,mod */
1201 /* internal variables */
1205 /* set info flag. */
1207 /* Parameter adjustments */
1209 x_offset
= x_dim1
+ 1;
1222 /* determine what is to be computed. */
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;
1230 i__1
= k
, i__2
= n
- 1;
1231 ju
= min(i__1
,i__2
);
1233 /* special action when n=1. */
1239 qy
[1].r
= y
[1].r
, qy
[1].i
= y
[1].i
;
1242 qty
[1].r
= y
[1].r
, qty
[1].i
= y
[1].i
;
1245 xb
[1].r
= y
[1].r
, xb
[1].i
= y
[1].i
;
1252 z__1
.r
= x
[i__2
].r
* 0. - x
[i__2
].i
* -1., z__1
.i
= x
[i__2
].r
* -1. + x
[
1254 if ((d__1
= x
[i__1
].r
, abs(d__1
)) + (d__2
= z__1
.r
, abs(d__2
)) != 0.) {
1260 z_div(&z__1
, &y
[1], &x
[x_dim1
+ 1]);
1261 b
[1].r
= z__1
.r
, b
[1].i
= z__1
.i
;
1265 rsd
[1].r
= 0., rsd
[1].i
= 0.;
1270 /* set up to compute qy or qty. */
1273 blas_zcopy(n
, &y
[1], 1, &qy
[1], 1);
1276 blas_zcopy(n
, &y
[1], 1, &qty
[1], 1);
1285 for (jj
= 1; jj
<= i__1
; ++jj
) {
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
)) ==
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
;
1299 x
[i__2
].r
= qraux
[i__3
].r
, x
[i__2
].i
= qraux
[i__3
].i
;
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
;
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
;
1318 /* compute ctrans(q)*y. */
1321 for (j
= 1; j
<= i__1
; ++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
)) ==
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
;
1334 x
[i__2
].r
= qraux
[i__3
].r
, x
[i__2
].i
= qraux
[i__3
].i
;
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
;
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
;
1350 /* set up to compute b, rsd, or xb. */
1353 blas_zcopy(k
, &qty
[1], 1, &b
[1], 1);
1357 blas_zcopy(k
, &qty
[1], 1, &xb
[1], 1);
1361 blas_zcopy(i__1
, &qty
[kp1
], 1, &rsd
[kp1
], 1);
1363 if (! cxb
|| kp1
> n
) {
1367 for (i
= kp1
; i
<= i__1
; ++i
) {
1369 xb
[i__2
].r
= 0., xb
[i__2
].i
= 0.;
1377 for (i
= 1; i
<= i__1
; ++i
) {
1379 rsd
[i__2
].r
= 0., rsd
[i__2
].i
= 0.;
1390 for (jj
= 1; jj
<= i__1
; ++jj
) {
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. +
1396 if ((d__1
= x
[i__2
].r
, abs(d__1
)) + (d__2
= z__1
.r
, abs(d__2
)) != 0.)
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
;
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
;
1414 blas_zaxpy(i__2
, &t
, &x
[j
* x_dim1
+ 1], 1, &b
[1], 1);
1421 if (! cr
&& ! cxb
) {
1425 /* compute rsd or xb as required. */
1428 for (jj
= 1; jj
<= i__1
; ++jj
) {
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
)) ==
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
;
1442 x
[i__2
].r
= qraux
[i__3
].r
, x
[i__2
].i
= qraux
[i__3
].i
;
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
;
1452 blas_zaxpy(i__2
, &t
, &x
[j
+ j
* x_dim1
], 1, &rsd
[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
;
1463 blas_zaxpy(i__2
, &t
, &x
[j
+ j
* x_dim1
], 1, &xb
[j
], 1);
1465 i__2
= j
+ j
* x_dim1
;
1466 x
[i__2
].r
= temp
.r
, x
[i__2
].i
= temp
.i
;