SUBROUTINE SQRLS (A,LDA,M,N,TOL,KR,B,X,RSD,JPVT,QRAUX,WORK, * ITASK,IND) C***BEGIN PROLOGUE SQRLS C***DATE WRITTEN 870911 (YYMMDD) C***REVISION DATE 871016 (YYMMDD) C***CATEGORY NO. D9 C***KEYWORDS LEAST SQUARES,OVERDETERMINED,LINEAR EQUATIONS C***AUTHOR STEPHEN NASH (GEORGE MASOM UNIVERSITY) C***PURPOSE C***PURPOSE SQRLS solves an overdetermined, underdetermined or singular C system of linear equations in least square sense. The C solution is obtained using a QR factorization of the C M by N coefficient matrix A. C***DESCRIPTION C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C C SQRLS IS USED TO SOLVE IN A LEAST SQUARES SENSE C OVERDETERMINED, UNDERDETERMINED AND SINGULAR LINEAR SYSTEMS . 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 MINIMIMZES THE SUM OF SQUARES (2-NORM) C OF THE RESIDUAL, A*X - B . C THE NUMERICAL RANK OF A IS DETERMINED USING THE TOLERANCE TOL. C C SQRLS USES THE LINPACK SUBROUTINE SQRDC TO COMPUTE THE QR C FACTORIZATION, WITH COLUMN PIVOTING, OF AN M BY N MATRIX A . C FOR MORE INFORMATION, SEE CHAPTER 9 OF THE REFERENCE BELOW. C C C ON ENTRY C C A REAL (LDA,N) . C THE MATRIX WHOSE DECOMPOSITION IS TO BE COMPUTED. C IN A LEAST SQUARES DATA FITTING PROBLEM, A(I,J) IS THE C VALUE OF THE J-TH BASIS (MODEL) FUNCTION AT THE I-TH C DATA POINT. 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 REAL. 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 QRAUX REAL(N) C WORK REAL(N) C THREE AUXILIARY ARRAYS USED TO FACTOR THE MATRIX A. C (NOT REQUIRED IF ITASK .GT. 1) C C B REAL(M) C THE RIGHT HAND SIDE OF THE LINEAR SYSTEM. C IN A LEAST SQUARES DATA FITTING PROBLEM B(I) CONTAINS THE C VALUE OF I-TH OBSERVATION. C C ITASK INTEGER. C IF ITASK=1, THEN SQRLS FACTORS THE MATRIX A AND C SOLVES THE LEAST SQUARES PROBLEM. C IF ITASK=2, THEN SQRLS ASSUMES THAT THE MATRIX A C WAS FACTORED WITH AN EARLIER CALL TO C SQRLS, AND ONLY SOLVES THE LEAST SQUARES C PROBLEM. C C ON RETURN C C X REAL(N) . C A LEAST SQUARES SOLUTION TO THE LINEAR SYSTEM. C C RSD REAL(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 SQRLS 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 SQRLS. C C***REFERENCES DONGARRA, ET AL, LINPACK USERS GUIDE, SIAM, 1979 C***ROUTINES CALLED SQRANK,SQRLSS,XERROR C***END PROLOGUE SQRLS INTEGER LDA, M, N, ITASK, JPVT(1), KR, IND REAL 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 SQRANK (A, LDA, M, N, TOL, KR, JPVT, QRAUX, WORK) C C SOLVE LEAST-SQUARES PROBLEM C CALL SQRLSS (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, '( * ''SQRLS 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, '( * ''SQRLS 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, '( * ''SQRLS ERROR (IND=-3) -- ITASK='', I5, '' IS LESS THAN 1.'') * ') ITASK CALL XERROR(MSG(1:51), 51, -3, 0) RETURN END SUBROUTINE SQRANK(A,LDA,M,N,TOL,KR,JPVT,QRAUX,WORK) C***BEGIN PROLOGUE SQRANK C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C***END PROLOGUE INTEGER LDA,M,N,KR,JPVT(1) REAL 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 SQRDC(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 SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,JOB) C***BEGIN PROLOGUE SQRDC C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C***END PROLOGUE SQRDC INTEGER LDX,N,P,JOB INTEGER JPVT(1) REAL X(LDX,1),QRAUX(1),WORK(1) C INTEGER J,JP,L,LP1,LUP,MAXJ,PL,PU REAL MAXNRM,SNRM2,TT REAL SDOT,NRMXL,T LOGICAL NEGJ,SWAPJ C C***FIRST EXECUTABLE STATEMENT SQRDC 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 SSWAP(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 SSWAP(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) = SNRM2(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.0E0 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 SSWAP(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.0E0 IF (L .EQ. N) GO TO 190 C C COMPUTE THE HOUSEHOLDER TRANSFORMATION FOR COLUMN L. C NRMXL = SNRM2(N-L+1,X(L,L),1) IF (NRMXL .EQ. 0.0E0) GO TO 180 IF (X(L,L) .NE. 0.0E0) NRMXL = SIGN(NRMXL,X(L,L)) CALL SSCAL(N-L+1,1.0E0/NRMXL,X(L,L),1) X(L,L) = 1.0E0 + 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 = -SDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L) CALL SAXPY(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.0E0) GO TO 150 TT = 1.0E0 - (ABS(X(L,J))/QRAUX(J))**2 TT = AMAX1(TT,0.0E0) T = TT TT = 1.0E0 + 0.05E0*TT*(QRAUX(J)/WORK(J))**2 IF (TT .EQ. 1.0E0) GO TO 130 QRAUX(J) = QRAUX(J)*SQRT(T) GO TO 140 130 CONTINUE QRAUX(J) = SNRM2(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 SQRLSS(A,LDA,M,N,KR,B,X,RSD,JPVT,QRAUX) C***BEGIN PROLOGUE SQRLSS C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C***END PROLOGUE INTEGER LDA,M,N,KR,JPVT(1) REAL A(LDA,1),B(1),X(1),RSD(1),QRAUX(1) C***FIRST EXECUTABLE STATEMENT IF (KR .NE. 0) 1 CALL SQRSL(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.0 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 SQRSL(X,LDX,N,K,QRAUX,Y,QY,QTY,B,RSD,XB,JOB,INFO) C***BEGIN PROLOGUE SQRSL C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C***END PROLOGUE SQRSL INTEGER LDX,N,K,JOB,INFO REAL X(LDX,1),QRAUX(1),Y(1),QY(1),QTY(1),B(1),RSD(1),XB(1) C INTEGER I,J,JJ,JU,KP1 REAL SDOT,T,TEMP LOGICAL CB,CQY,CQTY,CR,CXB C C SET INFO FLAG. C C***FIRST EXECUTABLE STATEMENT SQRSL 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.0E0) 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.0E0 GO TO 250 40 CONTINUE C C SET UP TO COMPUTE QY OR QTY. C IF (CQY) CALL SCOPY(N,Y,1,QY,1) IF (CQTY) CALL SCOPY(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.0E0) GO TO 50 TEMP = X(J,J) X(J,J) = QRAUX(J) T = -SDOT(N-J+1,X(J,J),1,QY(J),1)/X(J,J) CALL SAXPY(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.0E0) GO TO 80 TEMP = X(J,J) X(J,J) = QRAUX(J) T = -SDOT(N-J+1,X(J,J),1,QTY(J),1)/X(J,J) CALL SAXPY(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 SCOPY(K,QTY,1,B,1) KP1 = K + 1 IF (CXB) CALL SCOPY(K,QTY,1,XB,1) IF (CR .AND. K .LT. N) CALL SCOPY(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.0E0 110 CONTINUE 120 CONTINUE IF (.NOT.CR) GO TO 140 DO 130 I = 1, K RSD(I) = 0.0E0 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.0E0) 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 SAXPY(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.0E0) GO TO 220 TEMP = X(J,J) X(J,J) = QRAUX(J) IF (.NOT.CR) GO TO 200 T = -SDOT(N-J+1,X(J,J),1,RSD(J),1)/X(J,J) CALL SAXPY(N-J+1,T,X(J,J),1,RSD(J),1) 200 CONTINUE IF (.NOT.CXB) GO TO 210 T = -SDOT(N-J+1,X(J,J),1,XB(J),1)/X(J,J) CALL SAXPY(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