SUBROUTINE DQRLS (A,LDA,M,N,TOL,KR,B,X,RSD,JPVT,QRAUX,WORK, * ITASK,IND) C***BEGIN PROLOGUE DQRLS C***REVISION SEPTEMBER 11, 1987 C***CATEGORY NO. E2G1A C***KEYWORD(S) OVERDETERMINED,LEAST SQUARES,LINEAR EQUATIONS C***AUTHOR STEPHEN NASH (GEORGE MASON UNIVERSITY) C***DATE WRITTEN SEPTEMBER 11, 1987 C***PURPOSE C SOLVES AN OVERDETERMINED, UNDERDETERMINED, OR SINGULAR SYSTEM OF C LINEAR EQUATIONS IN LEAST SQUARE SENSE. THE SOLUTION IS OBTAINED C USING A QR FACTORIZATION OF THE M BY N COEFFICIENT MATRIX A. C C***DESCRIPTION C DQRLS IS USED TO SOLVE OVERDETERMINED, UNDERDETERMINED AND C SINGULAR LINEAR SYSTEMS IN A LEAST SQUARES SENSE. C THE SYSTEM IS A*X APPROXIMATES B WHERE A IS M BY N. C B IS A GIVEN M-VECTOR, AND X IS THE N-VECTOR TO BE COMPUTED. C A SOLUTION X IS FOUND WHICH MINIMIZES THE SUM OF SQUARES (2-NORM) C OF THE RESIDUAL, A*X - B . C DQRLS USES THE LINPACK SUBROUTINE SQRDC TO COMPUTE THE QR C FACTORIZATION, WITH COLUMN PIVOTING, OF AN M BY N MATRIX A . C THE NUMERICAL RANK IS DETERMINED USING THE TOLERANCE TOL. C C FOR MORE INFORMATION, SEE CHAPTER 9 OF LINPACK USERS GUIDE, C J. DONGARRA ET ALL, SIAM, 1979. C C C ON ENTRY C C A DOUBLE PRECISION (LDA,N) . C THE MATRIX WHOSE DECOMPOSITION IS TO BE COMPUTED. C C LDA INTEGER. C THE LEADING DIMENSION OF A . C C M INTEGER. C THE NUMBER OF ROWS OF A . C C N INTEGER. C THE NUMBER OF COLUMNS OF A . C C TOL DOUBLE PRECISION. C A RELATIVE TOLERANCE USED TO DETERMINE THE NUMERICAL C RANK. THE PROBLEM SHOULD BE SCALED SO THAT ALL THE C ELEMENTS OF A HAVE ROUGHLY THE SAME ABSOLUTE ACCURACY C EPS. THEN A REASONABLE VALUE FOR TOL IS ROUGHLY EPS C DIVIDED BY THE MAGNITUDE OF THE LARGEST ELEMENT. C C JPVT INTEGER(N) C C QRAUX DOUBLE PRECISION(N) C C WORK DOUBLE PRECISION(N) C THREE AUXILIARY VECTORS USED TO FACTOR THE MATRIX A. C (ARRAY WORK IS NOT REQUIRED IF ITASK .GT. 1) C C B DOUBLE PRECISION(M) C THE RIGHT HAND SIDE OF THE LINEAR SYSTEM. C C ITASK INTEGER. C IF ITASK=1, THEN DQRLS FACTORS THE MATRIX A AND C SOLVES THE LEAST SQUARES PROBLEM. C IF ITASK>1, THEN DQRLS ASSUMES THAT THE MATRIX A C WAS FACTORED WITH AN EARLIER CALL TO C DQRLS, AND ONLY SOLVES THE LEAST SQUARES C PROBLEM. C C ON RETURN C C X DOUBLE PRECISION(N) . C A LEAST SQUARES SOLUTION TO THE LINEAR SYSTEM. C C RSD DOUBLE PRECISION(M) . C THE RESIDUAL, B - A*X . RSD MAY OVERWRITE B . C C IND INTEGER C ERROR CODE: IND = 0: NO ERROR C IND = -1: N .GT. LDA (FATAL ERROR) C IND = -2: N .LT. 1 (FATAL ERROR) C IND = -3: ITASK .LT. 1 (FATAL ERROR) C C A CONTAINS THE OUTPUT FROM SQRDC. C THE TRIANGULAR MATRIX R OF THE QR FACTORIZATION IS C CONTAINED IN THE UPPER TRIANGLE AND INFORMATION NEEDED C TO RECOVER THE ORTHOGONAL MATRIX Q IS STORED BELOW C THE DIAGONAL IN A AND IN THE VECTOR QRAUX . C C KR INTEGER. C THE NUMERICAL RANK. C C JPVT THE PIVOT INFORMATION FROM SQRDC. C C COLUMNS JPVT(1),...,JPVT(KR) OF THE ORIGINAL MATRIX ARE LINEARLY C INDEPENDENT TO WITHIN THE TOLERANCE TOL AND THE REMAINING COLUMNS C ARE LINEARLY DEPENDENT. ABS(A(1,1))/ABS(A(KR,KR)) IS AN ESTIMATE C OF THE CONDITION NUMBER OF THE MATRIX OF INDEPENDENT COLUMNS, C AND OF R . THIS ESTIMATE WILL BE .LE. 1/TOL . C C USAGE.... C DQRLS CAN BE EFFICIENTLY USED TO SOLVE SEVERAL LEAST SQUARES C PROBLEMS WITH THE SAME MATRIX A. THE FIRST SYSTEM IS SOLVED C WITH ITASK = 1. THE SUBSEQUENT SYSTEMS ARE SOLVED WITH C ITASK = 2, TO AVOID THE RECOMPUTATION OF THE MATRIX FACTORS. C THE PARAMETERS KR, JPVT, AND QRAUX MUST NOT BE MODIFIED C BETWEEN CALLS TO DQRLS. C C***REFERENCE(S) C DONGARRA, ET AL, LINPACK USERS GUIDE, SIAM, 1979 C***ROUTINES CALLED DQRANK, DQRLSS, XERROR C***END PROLOGUE INTEGER LDA, M, N, ITASK, JPVT(1), KR, IND DOUBLE PRECISION A(LDA,1), B(1), X(1), RSD(1), QRAUX(1), * WORK(1), TOL CHARACTER MSG*54 C***FIRST EXECUTABLE STATEMENT IF (LDA .LT. N) GO TO 10 IF (N .LE. 0) GO TO 20 IF (ITASK .LT. 1) GO TO 30 IND = 0 C IF (ITASK .EQ. 1) C C FACTOR MATRIX C * CALL DQRANK (A, LDA, M, N, TOL, KR, JPVT, QRAUX, WORK) C C SOLVE LEAST-SQUARES PROBLEM C CALL DQRLSS (A, LDA, M, N, KR, B, X, RSD, JPVT, QRAUX) RETURN C C ERROR IN CALLING SEQUENCE C C LDA .LT. N 10 IND = -1 WRITE (MSG, '( * ''DQRLS ERROR (IND=-1) -- LDA='', I5, '' IS LESS THAN N='', * I5 )' ) LDA, N CALL XERROR(MSG(1:54), 54, -1, 0) RETURN C C N .LT. 1 20 IND = -2 WRITE (MSG, '( * ''DQRLS ERROR (IND=-2) -- N='', I5, '' IS LESS THAN 1.'') ')N CALL XERROR(MSG(1:47), 47, -2, 0) RETURN C C ITASK .LT. 1 30 IND = -3 WRITE (MSG, '( * ''DQRLS ERROR (IND=-3) -- ITASK='', I5, '' IS LESS THAN 1.'') * ') ITASK CALL XERROR(MSG(1:51), 51, -3, 0) RETURN END SUBROUTINE DQRANK(A,LDA,M,N,TOL,KR,JPVT,QRAUX,WORK) C***BEGIN PROLOGUE DQRANK C***REVISION NOVEMBER 15, 1980 C***CATEGORY NO. E2G1A C***KEYWORD(S) OVERDETERMINED,LEAST SQUARES,LINEAR EQUATIONS C***AUTHOR D. KAHANER (NBS) C***DATE WRITTEN C***PURPOSE C COMPUTES THE QR FACTORIZATION OF AN M BY N MATRIX A. THIS C ROUTINE IS USED IN CONJUNCTION WITH DQRLSS TO SOLVE LINEAR C SYSTEMS OF EQUATIONS IN A LEAST SQUARE SENSE. C***DESCRIPTION C C DQRANK IS USED IN CONJUNCTION WITH DQRLSS TO SOLVE C OVERDETERMINED, UNDERDETERMINED AND SINGULAR LINEAR SYSTEMS C IN A LEAST SQUARES SENSE. C DQRANK USES THE LINPACK SUBROUTINE DQRDC TO COMPUTE THE QR C FACTORIZATION, WITH COLUMN PIVOTING, OF AN M BY N MATRIX A . C THE NUMERICAL RANK IS DETERMINED USING THE TOLERANCE TOL. C C FOR MORE INFORMATION, SEE CHAPTER 9 OF LINPACK USERS GUIDE, C J. DONGARRA ET ALL, SIAM, 1979. C C ON ENTRY C C A DOUBLE PRECISION (LDA,N) . C THE MATRIX WHOSE DECOMPOSITION IS TO BE COMPUTED. C C LDA INTEGER. C THE LEADING DIMENSION OF A . C C M INTEGER. C THE NUMBER OF ROWS OF A . C C N INTEGER. C THE NUMBER OF COLUMNS OF A . C C TOL DOUBLE PRECISION. C A RELATIVE TOLERANCE USED TO DETERMINE THE NUMERICAL C RANK. THE PROBLEM SHOULD BE SCALED SO THAT ALL THE C ELEMENTS OF A HAVE ROUGHLY THE SAME ABSOLUTE ACCURACY C EPS. THEN A REASONABLE VALUE FOR TOL IS ROUGHLY EPS C DIVIDED BY THE MAGNITUDE OF THE LARGEST ELEMENT. C C JPVT INTEGER(N) C C QRAUX DOUBLE PRECISION(N) C C WORK DOUBLE PRECISION(N) C THREE AUXILLIARY VECTORS USED BY DQRDC . C C ON RETURN C C A CONTAINS THE OUTPUT FROM DQRDC. C THE TRIANGULAR MATRIX R OF THE QR FACTORIZATION IS C CONTAINED IN THE UPPER TRIANGLE AND INFORMATION NEEDED C TO RECOVER THE ORTHOGONAL MATRIX Q IS STORED BELOW C THE DIAGONAL IN A AND IN THE VECTOR QRAUX . C C KR INTEGER. C THE NUMERICAL RANK. C C JPVT THE PIVOT INFORMATION FROM DQRDC. C C COLUMNS JPVT(1),...,JPVT(KR) OF THE ORIGINAL MATRIX ARE LINEARLY C INDEPENDENT TO WITHIN THE TOLERANCE TOL AND THE REMAINING COLUMNS C ARE LINEARLY DEPENDENT. ABS(A(1,1))/ABS(A(KR,KR)) IS AN ESTIMATE C OF THE CONDITION NUMBER OF THE MATRIX OF INDEPENDENT COLUMNS, C AND OF R . THIS ESTIMATE WILL BE .LE. 1/TOL . C C USAGE.....SEE SUBROUTINE DQRLSS C C***REFERENCE(S) C DONGARRA, ET AL, LINPACK USERS GUIDE, SIAM, 1979 C***ROUTINES CALLED DQRDC C***END PROLOGUE INTEGER LDA,M,N,KR,JPVT(1) DOUBLE PRECISION A(LDA,1),TOL,QRAUX(1),WORK(1) INTEGER J,K C***FIRST EXECUTABLE STATEMENT DO 10 J = 1, N JPVT(J) = 0 10 CONTINUE CALL DQRDC(A,LDA,M,N,QRAUX,JPVT,WORK,1) KR = 0 K = MIN0(M,N) DO 20 J = 1, K IF (ABS(A(J,J)) .LE. TOL*ABS(A(1,1))) GO TO 30 KR = J 20 CONTINUE 30 RETURN END SUBROUTINE DQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,JOB) C***BEGIN PROLOGUE DQRDC C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 861211 (YYMMDD) C***CATEGORY NO. D5 C***KEYWORDS LIBRARY=SLATEC(LINPACK), C TYPE=DOUBLE PRECISION(SQRDC-S DQRDC-D CQRDC-C), C LINEAR ALGEBRA,MATRIX,ORTHOGONAL TRIANGULAR, C QR DECOMPOSITION C***AUTHOR STEWART, G. W., (U. OF MARYLAND) C***PURPOSE Uses Householder transformations to compute the Qr factori- C zation of N by P matrix X. Column pivoting is optional. C***DESCRIPTION C C DQRDC uses Householder transformations to compute the QR C factorization of an N by P matrix X. Column pivoting C based on the 2-norms of the reduced columns may be C performed at the user's option. C C On Entry C C X DOUBLE PRECISION(LDX,P), where LDX .GE. N. C X contains the matrix whose decomposition is to be C computed. C C LDX INTEGER. C LDX is the leading dimension of the array X. C C N INTEGER. C N is the number of rows of the matrix X. C C P INTEGER. C P is the number of columns of the matrix X. C C JPVT INTEGER(P). C JPVT contains integers that control the selection C of the pivot columns. The K-th column X(K) of X C is placed in one of three classes according to the C value of JPVT(K). C C If JPVT(K) .GT. 0, then X(K) is an initial C column. C C If JPVT(K) .EQ. 0, then X(K) is a free column. C C If JPVT(K) .LT. 0, then X(K) is a final column. C C Before the decomposition is computed, initial columns C are moved to the beginning of the array X and final C columns to the end. Both initial and final columns C are frozen in place during the computation and only C free columns are moved. At the K-th stage of the C reduction, if X(K) is occupied by a free column C it is interchanged with the free column of largest C reduced norm. JPVT is not referenced if C JOB .EQ. 0. C C WORK DOUBLE PRECISION(P). C WORK is a work array. WORK is not referenced if C JOB .EQ. 0. C C JOB INTEGER. C JOB is an integer that initiates column pivoting. C If JOB .EQ. 0, no pivoting is done. C If JOB .NE. 0, pivoting is done. C C On Return C C X X contains in its upper triangle the upper C triangular matrix R of the QR factorization. C Below its diagonal X contains information from C which the orthogonal part of the decomposition C can be recovered. Note that if pivoting has C been requested, the decomposition is not that C of the original matrix X but that of X C with its columns permuted as described by JPVT. C C QRAUX DOUBLE PRECISION(P). C QRAUX contains further information required to recover C the orthogonal part of the decomposition. C C JPVT JPVT(K) contains the index of the column of the C original matrix that has been interchanged into C the K-th column, if pivoting was requested. C C LINPACK. This version dated 08/14/78 . C G. W. Stewart, University of Maryland, Argonne National Lab. C C DQRDC uses the following functions and subprograms. C C BLAS DAXPY,DDOT,DSCAL,DSWAP,DNRM2 C Fortran DABS,DMAX1,MIN0,DSQRT C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED DAXPY,DDOT,DNRM2,DSCAL,DSWAP C***END PROLOGUE DQRDC INTEGER LDX,N,P,JOB INTEGER JPVT(1) DOUBLE PRECISION X(LDX,1),QRAUX(1),WORK(1) C INTEGER J,JP,L,LP1,LUP,MAXJ,PL,PU DOUBLE PRECISION MAXNRM,DNRM2,TT DOUBLE PRECISION DDOT,NRMXL,T LOGICAL NEGJ,SWAPJ C C***FIRST EXECUTABLE STATEMENT DQRDC PL = 1 PU = 0 IF (JOB .EQ. 0) GO TO 60 C C PIVOTING HAS BEEN REQUESTED. REARRANGE THE COLUMNS C ACCORDING TO JPVT. C DO 20 J = 1, P SWAPJ = JPVT(J) .GT. 0 NEGJ = JPVT(J) .LT. 0 JPVT(J) = J IF (NEGJ) JPVT(J) = -J IF (.NOT.SWAPJ) GO TO 10 IF (J .NE. PL) CALL DSWAP(N,X(1,PL),1,X(1,J),1) JPVT(J) = JPVT(PL) JPVT(PL) = J PL = PL + 1 10 CONTINUE 20 CONTINUE PU = P DO 50 JJ = 1, P J = P - JJ + 1 IF (JPVT(J) .GE. 0) GO TO 40 JPVT(J) = -JPVT(J) IF (J .EQ. PU) GO TO 30 CALL DSWAP(N,X(1,PU),1,X(1,J),1) JP = JPVT(PU) JPVT(PU) = JPVT(J) JPVT(J) = JP 30 CONTINUE PU = PU - 1 40 CONTINUE 50 CONTINUE 60 CONTINUE C C COMPUTE THE NORMS OF THE FREE COLUMNS. C IF (PU .LT. PL) GO TO 80 DO 70 J = PL, PU QRAUX(J) = DNRM2(N,X(1,J),1) WORK(J) = QRAUX(J) 70 CONTINUE 80 CONTINUE C C PERFORM THE HOUSEHOLDER REDUCTION OF X. C LUP = MIN0(N,P) DO 200 L = 1, LUP IF (L .LT. PL .OR. L .GE. PU) GO TO 120 C C LOCATE THE COLUMN OF LARGEST NORM AND BRING IT C INTO THE PIVOT POSITION. C MAXNRM = 0.0D0 MAXJ = L DO 100 J = L, PU IF (QRAUX(J) .LE. MAXNRM) GO TO 90 MAXNRM = QRAUX(J) MAXJ = J 90 CONTINUE 100 CONTINUE IF (MAXJ .EQ. L) GO TO 110 CALL DSWAP(N,X(1,L),1,X(1,MAXJ),1) QRAUX(MAXJ) = QRAUX(L) WORK(MAXJ) = WORK(L) JP = JPVT(MAXJ) JPVT(MAXJ) = JPVT(L) JPVT(L) = JP 110 CONTINUE 120 CONTINUE QRAUX(L) = 0.0D0 IF (L .EQ. N) GO TO 190 C C COMPUTE THE HOUSEHOLDER TRANSFORMATION FOR COLUMN L. C NRMXL = DNRM2(N-L+1,X(L,L),1) IF (NRMXL .EQ. 0.0D0) GO TO 180 IF (X(L,L) .NE. 0.0D0) NRMXL = DSIGN(NRMXL,X(L,L)) CALL DSCAL(N-L+1,1.0D0/NRMXL,X(L,L),1) X(L,L) = 1.0D0 + X(L,L) C C APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS, C UPDATING THE NORMS. C LP1 = L + 1 IF (P .LT. LP1) GO TO 170 DO 160 J = LP1, P T = -DDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L) CALL DAXPY(N-L+1,T,X(L,L),1,X(L,J),1) IF (J .LT. PL .OR. J .GT. PU) GO TO 150 IF (QRAUX(J) .EQ. 0.0D0) GO TO 150 TT = 1.0D0 - (DABS(X(L,J))/QRAUX(J))**2 TT = DMAX1(TT,0.0D0) T = TT TT = 1.0D0 + 0.05D0*TT*(QRAUX(J)/WORK(J))**2 IF (TT .EQ. 1.0D0) GO TO 130 QRAUX(J) = QRAUX(J)*DSQRT(T) GO TO 140 130 CONTINUE QRAUX(J) = DNRM2(N-L,X(L+1,J),1) WORK(J) = QRAUX(J) 140 CONTINUE 150 CONTINUE 160 CONTINUE 170 CONTINUE C C SAVE THE TRANSFORMATION. C QRAUX(L) = X(L,L) X(L,L) = -NRMXL 180 CONTINUE 190 CONTINUE 200 CONTINUE RETURN END SUBROUTINE DQRLSS(A,LDA,M,N,KR,B,X,RSD,JPVT,QRAUX) C***BEGIN PROLOGUE DQRLSS C***REVISION NOVEMBER 15, 1980 C***CATEGORY NO. E2G1A C***KEYWORD(S) OVERDETERMINED,LEAST SQUARES,LINEAR EQUATIONS C***AUTHOR D. KAHANER (NBS) C***DATE WRITTEN C***PURPOSE C SOLVES AN OVERDETERMINED, UNDERDETERMINED, OR SINGULAR SYSTEM OF C LINEAR EQUATIONS IN LEAST SQUARE SENSE. C***DESCRIPTION C C DQRLSS USES THE LINPACK SUBROUTINE DQRSL TO SOLVE AN C OVERDETERMINED, UNDERDETERMINED, OR SINGULAR LINEAR SYSTEM IN A C LEAST SQUARES SENSE. IT MUST BE PRECEEDED BY DQRANK . C THE SYSTEM IS A*X APPROXIMATES B WHERE A IS M BY N WITH C RANK KR (DETERMINED BY DQRANK), B IS A GIVEN M-VECTOR, C AND X IS THE N-VECTOR TO BE COMPUTED. A SOLUTION X WITH C AT MOST KR NONZERO COMPONENTS IS FOUND WHICH MINIMIMZES C THE 2-NORM OF THE RESIDUAL, A*X - B . C C ON ENTRY C C A,LDA,M,N,KR,JPVT,QRAUX C THE OUTPUT FROM DQRANK . C C B DOUBLE PRECISION(M) . C THE RIGHT HAND SIDE OF THE LINEAR SYSTEM. C C ON RETURN C C X DOUBLE PRECISION(N) . C A LEAST SQUARES SOLUTION TO THE LINEAR SYSTEM. C C RSD DOUBLE PRECISION(M) . C THE RESIDUAL, B - A*X . RSD MAY OVERWITE B . C C USAGE.... C ONCE THE MATRIX A HAS BEEN FORMED, DQRANK SHOULD BE C CALLED ONCE TO DECOMPOSE IT. THEN FOR EACH RIGHT HAND C SIDE, B, DQRLSS SHOULD BE CALLED ONCE TO OBTAIN THE C SOLUTION AND RESIDUAL. C C INTEGER J,K,INFO C DOUBLE PRECISION T C C***REFERENCE(S) C DONGARRA, ET AL, LINPACK USERS GUIDE, SIAM, 1979 C***ROUTINES CALLED DQRSL C***END PROLOGUE INTEGER LDA,M,N,KR,JPVT(1) DOUBLE PRECISION A(LDA,1),B(1),X(1),RSD(1),QRAUX(1) C***FIRST EXECUTABLE STATEMENT IF (KR .NE. 0) 1 CALL DQRSL(A,LDA,M,KR,QRAUX,B,RSD,RSD,X,RSD,RSD,110,INFO) DO 40 J = 1, N JPVT(J) = -JPVT(J) IF (J .GT. KR) X(J) = 0.0D0 40 CONTINUE DO 70 J = 1, N IF (JPVT(J) .GT. 0) GO TO 70 K = -JPVT(J) JPVT(J) = K 50 CONTINUE IF (K .EQ. J) GO TO 60 T = X(J) X(J) = X(K) X(K) = T JPVT(K) = -JPVT(K) K = JPVT(K) GO TO 50 60 CONTINUE 70 CONTINUE RETURN END SUBROUTINE DQRSL(X,LDX,N,K,QRAUX,Y,QY,QTY,B,RSD,XB,JOB,INFO) C***BEGIN PROLOGUE DQRSL C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 861211 (YYMMDD) C***CATEGORY NO. D9,D2A1 C***KEYWORDS LIBRARY=SLATEC(LINPACK), C TYPE=DOUBLE PRECISION(SQRSL-S DQRSL-D CQRSL-C), C LINEAR ALGEBRA,MATRIX,ORTHOGONAL TRIANGULAR,SOLVE C***AUTHOR STEWART, G. W., (U. OF MARYLAND) C***PURPOSE Applies the output of DQRDC to compute coordinate C transformations, projections, and least squares solutions. C***DESCRIPTION C C DQRSL applies the output of DQRDC to compute coordinate C transformations, projections, and least squares solutions. C For K .LE. MIN(N,P), let XK be the matrix C C XK = (X(JPVT(1)),X(JPVT(2)), ... ,X(JPVT(K))) C C formed from columnns JPVT(1), ... ,JPVT(K) of the original C N X P matrix X that was input to DQRDC (if no pivoting was C done, XK consists of the first K columns of X in their C original order). DQRDC produces a factored orthogonal matrix Q C and an upper triangular matrix R such that C C XK = Q * (R) C (0) C C This information is contained in coded form in the arrays C X and QRAUX. C C On Entry C C X DOUBLE PRECISION(LDX,P). C X contains the output of DQRDC. C C LDX INTEGER. C LDX is the leading dimension of the array X. C C N INTEGER. C N is the number of rows of the matrix XK. It must C have the same value as N in DQRDC. C C K INTEGER. C K is the number of columns of the matrix XK. K C must not be greater than MIN(N,P), where P is the C same as in the calling sequence to DQRDC. C C QRAUX DOUBLE PRECISION(P). C QRAUX contains the auxiliary output from DQRDC. C C Y DOUBLE PRECISION(N) C Y contains an N-vector that is to be manipulated C by DQRSL. C C JOB INTEGER. C JOB specifies what is to be computed. JOB has C the decimal expansion ABCDE, with the following C meaning. C C If A .NE. 0, compute QY. C If B,C,D, or E .NE. 0, compute QTY. C If C .NE. 0, compute B. C If D .NE. 0, compute RSD. C If E .NE. 0, compute XB. C C Note that a request to compute B, RSD, or XB C automatically triggers the computation of QTY, for C which an array must be provided in the calling C sequence. C C On Return C C QY DOUBLE PRECISION(N). C QY contains Q*Y, if its computation has been C requested. C C QTY DOUBLE PRECISION(N). C QTY contains TRANS(Q)*Y, if its computation has C been requested. Here TRANS(Q) is the C transpose of the matrix Q. C C B DOUBLE PRECISION(K) C B contains the solution of the least squares problem C C minimize norm2(Y - XK*B), C C if its computation has been requested. (Note that C if pivoting was requested in DQRDC, the J-th C component of B will be associated with column JPVT(J) C of the original matrix X that was input into DQRDC.) C C RSD DOUBLE PRECISION(N). C RSD contains the least squares residual Y - XK*B, C if its computation has been requested. RSD is C also the orthogonal projection of Y onto the C orthogonal complement of the column space of XK. C C XB DOUBLE PRECISION(N). C XB contains the least squares approximation XK*B, C if its computation has been requested. XB is also C the orthogonal projection of Y onto the column space C of X. C C INFO INTEGER. C INFO is zero unless the computation of B has C been requested and R is exactly singular. In C this case, INFO is the index of the first zero C diagonal element of R and B is left unaltered. C C The parameters QY, QTY, B, RSD, and XB are not referenced C if their computation is not requested and in this case C can be replaced by dummy variables in the calling program. C To save storage, the user may in some cases use the same C array for different parameters in the calling sequence. A C frequently occuring example is when one wishes to compute C any of B, RSD, or XB and does not need Y or QTY. In this C case one may identify Y, QTY, and one of B, RSD, or XB, while C providing separate arrays for anything else that is to be C computed. Thus the calling sequence C C CALL DQRSL(X,LDX,N,K,QRAUX,Y,DUM,Y,B,Y,DUM,110,INFO) C C will result in the computation of B and RSD, with RSD C overwriting Y. More generally, each item in the following C list contains groups of permissible identifications for C a single calling sequence. C C 1. (Y,QTY,B) (RSD) (XB) (QY) C C 2. (Y,QTY,RSD) (B) (XB) (QY) C C 3. (Y,QTY,XB) (B) (RSD) (QY) C C 4. (Y,QY) (QTY,B) (RSD) (XB) C C 5. (Y,QY) (QTY,RSD) (B) (XB) C C 6. (Y,QY) (QTY,XB) (B) (RSD) C C In any group the value returned in the array allocated to C the group corresponds to the last member of the group. C C LINPACK. This version dated 08/14/78 . C G. W. Stewart, University of Maryland, Argonne National Lab. C C DQRSL uses the following functions and subprograms. C C BLAS DAXPY,DCOPY,DDOT C Fortran DABS,MIN0,MOD C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED DAXPY,DCOPY,DDOT C***END PROLOGUE DQRSL INTEGER LDX,N,K,JOB,INFO DOUBLE PRECISION X(LDX,1),QRAUX(1),Y(1),QY(1),QTY(1),B(1),RSD(1), 1 XB(1) C INTEGER I,J,JJ,JU,KP1 DOUBLE PRECISION DDOT,T,TEMP LOGICAL CB,CQY,CQTY,CR,CXB C C SET INFO FLAG. C C***FIRST EXECUTABLE STATEMENT DQRSL INFO = 0 C C DETERMINE WHAT IS TO BE COMPUTED. C CQY = JOB/10000 .NE. 0 CQTY = MOD(JOB,10000) .NE. 0 CB = MOD(JOB,1000)/100 .NE. 0 CR = MOD(JOB,100)/10 .NE. 0 CXB = MOD(JOB,10) .NE. 0 JU = MIN0(K,N-1) C C SPECIAL ACTION WHEN N=1. C IF (JU .NE. 0) GO TO 40 IF (CQY) QY(1) = Y(1) IF (CQTY) QTY(1) = Y(1) IF (CXB) XB(1) = Y(1) IF (.NOT.CB) GO TO 30 IF (X(1,1) .NE. 0.0D0) GO TO 10 INFO = 1 GO TO 20 10 CONTINUE B(1) = Y(1)/X(1,1) 20 CONTINUE 30 CONTINUE IF (CR) RSD(1) = 0.0D0 GO TO 250 40 CONTINUE C C SET UP TO COMPUTE QY OR QTY. C IF (CQY) CALL DCOPY(N,Y,1,QY,1) IF (CQTY) CALL DCOPY(N,Y,1,QTY,1) IF (.NOT.CQY) GO TO 70 C C COMPUTE QY. C DO 60 JJ = 1, JU J = JU - JJ + 1 IF (QRAUX(J) .EQ. 0.0D0) GO TO 50 TEMP = X(J,J) X(J,J) = QRAUX(J) T = -DDOT(N-J+1,X(J,J),1,QY(J),1)/X(J,J) CALL DAXPY(N-J+1,T,X(J,J),1,QY(J),1) X(J,J) = TEMP 50 CONTINUE 60 CONTINUE 70 CONTINUE IF (.NOT.CQTY) GO TO 100 C C COMPUTE TRANS(Q)*Y. C DO 90 J = 1, JU IF (QRAUX(J) .EQ. 0.0D0) GO TO 80 TEMP = X(J,J) X(J,J) = QRAUX(J) T = -DDOT(N-J+1,X(J,J),1,QTY(J),1)/X(J,J) CALL DAXPY(N-J+1,T,X(J,J),1,QTY(J),1) X(J,J) = TEMP 80 CONTINUE 90 CONTINUE 100 CONTINUE C C SET UP TO COMPUTE B, RSD, OR XB. C IF (CB) CALL DCOPY(K,QTY,1,B,1) KP1 = K + 1 IF (CXB) CALL DCOPY(K,QTY,1,XB,1) IF (CR .AND. K .LT. N) CALL DCOPY(N-K,QTY(KP1),1,RSD(KP1),1) IF (.NOT.CXB .OR. KP1 .GT. N) GO TO 120 DO 110 I = KP1, N XB(I) = 0.0D0 110 CONTINUE 120 CONTINUE IF (.NOT.CR) GO TO 140 DO 130 I = 1, K RSD(I) = 0.0D0 130 CONTINUE 140 CONTINUE IF (.NOT.CB) GO TO 190 C C COMPUTE B. C DO 170 JJ = 1, K J = K - JJ + 1 IF (X(J,J) .NE. 0.0D0) GO TO 150 INFO = J C ......EXIT GO TO 180 150 CONTINUE B(J) = B(J)/X(J,J) IF (J .EQ. 1) GO TO 160 T = -B(J) CALL DAXPY(J-1,T,X(1,J),1,B,1) 160 CONTINUE 170 CONTINUE 180 CONTINUE 190 CONTINUE IF (.NOT.CR .AND. .NOT.CXB) GO TO 240 C C COMPUTE RSD OR XB AS REQUIRED. C DO 230 JJ = 1, JU J = JU - JJ + 1 IF (QRAUX(J) .EQ. 0.0D0) GO TO 220 TEMP = X(J,J) X(J,J) = QRAUX(J) IF (.NOT.CR) GO TO 200 T = -DDOT(N-J+1,X(J,J),1,RSD(J),1)/X(J,J) CALL DAXPY(N-J+1,T,X(J,J),1,RSD(J),1) 200 CONTINUE IF (.NOT.CXB) GO TO 210 T = -DDOT(N-J+1,X(J,J),1,XB(J),1)/X(J,J) CALL DAXPY(N-J+1,T,X(J,J),1,XB(J),1) 210 CONTINUE X(J,J) = TEMP 220 CONTINUE 230 CONTINUE 240 CONTINUE 250 CONTINUE RETURN END