(Message inbox:44) Received: from PACIFIC-CARRIER-ANNEX.MIT.EDU by po10.MIT.EDU (5.61/4.7) id AA13142; Fri, 8 Jan 99 13:06:29 EST Received: from cyclops.ameslab.gov by MIT.EDU with SMTP id AA00872; Fri, 8 Jan 99 13:06:27 EST Received: (from wangcz@localhost) by cyclops.ameslab.gov (950413.SGI.8.6.12/950213.SGI.AUTOCF) id MAA05938 for liju99@mit.edu; Fri, 8 Jan 1999 12:06:20 -0600 Date: Fri, 8 Jan 1999 12:06:20 -0600 From: wangcz@cyclops.ameslab.gov (Cai-Zhuang Wang) Message-Id: <199901081806.MAA05938@cyclops.ameslab.gov> Apparently-To: liju99@mit.edu C ****************************************************************** C C TITLE: LINLSQ C C C LAST MODIFIED: 12 OCT., 1976. C C DESCRIPTION: C C THIS SUBROUTINE SOLVES THE WEIGHTED LEAST SQUARES PROBLEM: C M C MIN SUM ( WEIGHT(I) * (B(I)- F(X,T(I)))**2 ) C X I=1 C N C WHERE F(X,T(I))=SUM ( X(J)*PHI (T(I)) IS THE "MODEL" C J=1 J C WHICH IS LINEAR IN THE PARAMETERS X(J) (TO BE ESTIMATED) C AND T(I) IS THE VALUE OF THE INDEPENDENT VARIABLE OF THE C I-TH OBSERVATION B(I). C C THIS PROBLEM IS EQUIVALENT TO THE FOLLOWING PROBLEM: C C MIN ^^W(B-AX)^^ = MIN ^^WB-WAX^^ C X 2 X 2 C C WHERE W IS A DIAGONAL MATRIX WITH THE SQUARE ROOTS OF THE C WEIGHTS ON THE DIAGONAL ( THE USER ONLY HAS TO SUPPLY C THE 1-DIM ARRAY WEIGHT(I) OF THE WEIGHTS ) AND A IS C A M-BY-N MATRIX WITH ELEMENTS GIVEN BY A(I,J)=PHI (T(I)). C J C THE USER SHOULD FILL THE 2-DIM ARRAY A THIS WAY BEFORE C CALLING LINLSQ. C C ------------------------------------------- C C ALGORITHM: C C THE ROUTINE USES HOUSEHOLDER TRANSFORMATIONS TO REDUCE THE MATRIX C WA TO UPPER TRIANGULAR FORM : QWA=R WHERE Q IS ORTHOGONAL AND R IS C UPPER TRIANGULAR.(SUBROUTINE HECOMP). THEN IT USES BACK C SUBSTITUTION TO SOLVE FOR THE SOLUTION X FROM THE LINEAR SYSTEM C RX=P WHERE P IS THE FIRST N COMPONENTS OF THE VECTOR QWB. C (SUBROUTINE HOLVE). FOR MORE INFORMATION ON THE ALGORITHM, SEE C THE BOOK 'SOLVING LEAST SQUARES PROBLEMS' BY LAWSON AND HANSON, C PRENTICE-HALL,1974. C C ---------------------------------------------- C C HOW TO CALL: C C SINGLE RIGHT HAND SIDE B: CALL WITH MODE=1. C MULTIPLE RIGHT HAND SIDE B WITH THE SAME MATRIX A AND C THE SAME WEIGHTS: CALL WITH MODE=1 FOR THE FIRST C RIGHT HAND SIDE, AND THEN MAKE ALL SUBSEQUENT CALLS FOR C THE OTHER RIGHT HAND SIDES WITH MODE=2. C CALLS WITH MODE=2 ASSUME THAT THE MATRIX WA HAS C ALREADY BEEN DECOMPOSED USING THE HOUSEHOLDER C TRANSFORMATIONS AND THAT THE DECOMPOSITION IS STORED C OVER A, AND ONLY THE BACK SUBSTITUTION IS DONE FOR C THE NEW RIGHT HAND SIDE B. C C ----------------------------------------------- C C PARAMETERS: C C MODE VALUE 1 OR 2, ACCORDING TO THE NOTE ABOVE. C M NUMBER OF ROWS OF A (NUMBER OF OBSERVATIONS) C N NUMBER OF COLUMNS OF A (NUMBER OF PARAMETERS) C M.GE.N C A THE MATRIX A DEFINED ABOVE. ON OUTPUT, WILL C BE REPLACED BY Q AND R. C IDIM ROW DIMENSION OF A IN CALLING PROGRAM. C X N-VECTOR, ON OUTPUT WILL CONTAIN THE LEAST SQUARES C SOLUTION. C B M-VECTOR OF OBSERVATIONS. C ON OUTPUT, CONTAINS THE WEIGHTED RESIDUALS WB-WAX C WEIGHT M-VECTOR OF WEIGHTS (.GE.0). SHOULD BE SET C TO 1'S IF WEIGHTS ARE NOT DESIRED. C COVMAT N-BY-N MATRIX WHICH ON OUTPUT, CONTAINS THE C UNSCALED COVARIANCE MATRIX: C INVERSE( TRANSPOSE(WA) * WA ). C THIS MATRIX, OR A SCALAR MULTIPLE OF IT, HAS A C STATISTICAL INTERPRETATION, UNDER APPROPIATE C HYPOTHESIS, OF BEING AN ESTIMATE OF THE COVARIANCE C MATRIX FOR THE SOLUTION. (E.G. SEE PLACKETT,R.L. C (1960) PRINCIPLES OF REGRESSION ANALYSIS. OXFORD C UNIV. PRESS, N.Y. OR LAWSON,C.L. AND HANSON,R.J. C (1974) SOLVING LEAST SQUARES PROBLEMS. PRENTICE- C HALL, N.J.) C IN PARTICULAR, IF COVMAT IS MULTIPLIED BY THE C ESTIMATED VARIANCE (=RESNRM**2/(M-N)) THEN THE C DIAGONAL ELEMENTS OF COVMAT REPRESENT VARIANCES C AND THE OFF-DIAGONAL ELEMENTS COVARIANCES. C THE ACTUAL PARAMETER USED FOR COVMAT IN C THE CALLING PROGRAM CAN BE THE SAME AS THE C ACTUAL PARAMETER USED FOR A, BUT IF THAT C IS THE CASE, THE DECOMPOSITION OF A WILL BE C OVERWRITTEN BY COVMAT AND HENCE SUBSEQUENT C CALLS FOR MULTIPLE RIGHT HAND SIDES CANNOT C BE HANDLED. C ICMDIM ROW DIMENSION OF COVMAT IN CALLING PROGRAM. C RESNRM THE 2-NORM OF THE WEIGHTED RESIDUAL AT SOLUTION, C I.E. ^^WB-WAX^^ . C THE ESTIMATED VARIANCE = RESNRM**2/(M-N). C U M-VECTOR OF WORK SPACE. C IERR ERROR FLAG: C =0 ... NORMAL RETURN C =1 ... N.GT.M C =2 ... A SUSPECTED TO BE RANK DEFICIENT. C THE ABSOLUTE VALUE OF SOME DIAGONAL C ELEMENTS OF THE MATRIX R IS LESS THAN C EPS = SQRT(M*N)*MAX^WA(I,J)^*10**-8. C SUGGESTION: COMPUTE THE SINGULAR C VALUE DECOMPOSITION OF WA. C =3 ... SOME DIAGONAL ELEMENTS OF R ARE EXACTLY C ZERO AND HENCE A IS RANK DEFICIENT. C BACK SOLUTION NOT POSSIBLE. C SUGGESTION: COMPUTE THE SINGULAR VALUE C DECOMPOSITION OF WA. C =4 ... SOME WEIGHTS ARE NEGATIVE. C =5 ... IDIM.LT.M C C NOTE: ARRAYS A AND B ARE ALTERED BY LINLSQ. C IF LINLSQ IS TO BE CALLED WITH MODE=2, THEN THE ARRAYS C A,W AND U SHOULD NOT BE ALTERED BETWEEN CALLS OF LINLSQ. C ALSO,THE MOST RECENT CALL WITH MODE=1 IS ASSUMED TO HAVE C IERR=0,I.E. A NORMAL RETURN. C C PROGRAMMER: TONY CHAN, C.S.DEPT.,STANFORD UNIV. C C ************************************************************** C SUBROUTINE LINLSQ(MODE,M,N,A,IDIM,X,B,WEIGHT,COVMAT,ICMDIM, 1 RESNRM,U,IERR) INTEGER MODE,M,N,IDIM,ICMDIM,NP1 DOUBLE PRECISION A(IDIM,N),COVMAT(ICMDIM,N),X(N),B(M),U(M) DOUBLE PRECISION RESNRM,DSQRT,EPS,DABS,AMX,ABSAIJ DOUBLE PRECISION WEIGHT(M),ADIAG,SQRTWI C IERR=0 IF (N.GT.M) GO TO 10 IF (IDIM.LT.M) GO TO 90 IF (MODE.EQ.2) GO TO 80 C C COMPUTE WA. REPLACE A. C AMX=0.0D0 DO 60 I=1,M IF (WEIGHT(I).LT.0.0D0) GO TO 70 SQRTWI=DSQRT(WEIGHT(I)) DO 65 J=1,N A(I,J)=A(I,J)*SQRTWI ABSAIJ=DABS(A(I,J)) IF (ABSAIJ.GT.AMX) AMX=ABSAIJ 65 CONTINUE 60 CONTINUE C C TRANSFORM A C CALL HECOMP(IDIM,M,N,A,U) C C CALCULATE EPS FOR TESTING RANK-DEFICIENCY C EPS=DSQRT(M*N*1.0D0)*AMX*1.0D-08 C C TEST FOR SMALL DIAGONAL ELEMENTS OF R C DO 20 K=1,N ADIAG=DABS(A(K,K)) IF (ADIAG.LT.EPS) IERR=2 IF (ADIAG.EQ.0.0D0) GO TO 30 20 CONTINUE C C COMPUTE WB. REPLACE B. C 80 DO 25 I=1,M SQRTWI=DSQRT(WEIGHT(I)) 25 B(I)=B(I)*SQRTWI C C BACK SOLVE FOR X C CALL HOLVE(IDIM,M,N,A,U,B) C C SOLUTION IS NOW IN THE FIRST N COMPONENTS OF B C COPY INTO X AND COMPUTE 2-NORM OF WEIGHTED RESIDUALS C DO 40 K=1,N 40 X(K)=B(K) NP1=N+1 RESNRM=0.0D0 IF (M.EQ.N) GO TO 45 DO 50 K=NP1,M 50 RESNRM=RESNRM+B(K)**2 RESNRM=DSQRT(RESNRM) 45 CONTINUE C C COMPUTE THE WEIGHTED RESIDUALS C CALL RESID(IDIM,A,B,U,M,N) IF ( MODE .EQ. 2 ) GO TO 85 C C COMPUTE COVARIANCE MATRIX IF MODE=1. C CALL COVAR(ICMDIM,COVMAT,IDIM,A,N) 85 RETURN 10 IERR=1 RETURN 30 IERR=3 RETURN 70 IERR=4 RETURN 90 IERR=5 RETURN END SUBROUTINE COVAR(ICMDIM,COVMAT,IDIM,A,N) C C COMPUTES THE UNSCALED COVARIANCE MATRIX : C T -1 T -1 -1 -T C ((WA) *(WA)) =(R *R) =R *R C R IS STORED IN THE UPPER TRIANGULAR PART OF A. C NOTE THAT THE COVARIANCE MATRIX IS ALWAYS SYMMETRIC. C INTEGER ICMDIM,IDIM,N,NM1,IP1,JM1 DOUBLE PRECISION COVMAT(ICMDIM,N),A(IDIM,N) DOUBLE PRECISION SUM DO 10 J=1,N 10 COVMAT(J,J)=1.0D0/A(J,J) C C INVERT R AND STORE IN UPPER TRIANGULAR PART OF COVMAT C 33 IF (N.EQ.1) GO TO 70 NM1=N-1 DO 60 I=1,NM1 IP1=I+1 DO 60 J=IP1,N JM1=J-1 SUM=COVMAT(I,I)*A(I,J) IF (I.GE.JM1) GO TO 60 DO 50 L=IP1,JM1 50 SUM=SUM+COVMAT(I,L)*A(L,J) 60 COVMAT(I,J)=-SUM*COVMAT(J,J) C C FORM PRODUCT INVERSE(R) WITH ITS TRANSPOSE C 70 DO 90 I=1,N DO 90 J=I,N SUM=0.0D0 DO 80 L=J,N 80 SUM=SUM+COVMAT(I,L)*COVMAT(J,L) COVMAT(I,J)=SUM 90 COVMAT(J,I)=SUM RETURN END SUBROUTINE RESID(IDIM,A,B,U,M,N) C C COMPUTES THE WEIGHTED RESIDUAL WB-WAX BY UNWINDING THE C TRANSFORMATIONS, I.E. APPLY TRANSPOSE(Q) TO QWB (WHICH IS C STORED IN B NOW). THE WEIGHTED RESIDUALS OVERWRITE B. C INTEGER IDIM,M,N,NP1,K,KBACK,KP1 DOUBLE PRECISION A(IDIM,N),B(M),U(M) DOUBLE PRECISION SUM C C COMPUTE WEIGHTED RESIDUALS WB-WAX BY UNWINDING THE TRANSFORMATIONS C IF (M.GT.N) GO TO 80 C M=N MEANS WE ARE SOLVING A LINEAR SYSTEM AND SO RESIDUAL = 0. DO 70 I=1,N 70 B(I)=0.0D0 RETURN 80 NP1=N+1 DO 98 KBACK=1,N K=NP1-KBACK KP1=K+1 SUM=0.0D0 DO 97 I=KP1,M 97 SUM=SUM+A(I,K)*B(I) SUM=-SUM/(U(K)*A(K,K)) B(K)=-U(K)*SUM DO 98 I=KP1,M 98 B(I)=B(I)-A(I,K)*SUM RETURN END SUBROUTINE HECOMP(MDIM,M,N,A,U) INTEGER MDIM,M,N DOUBLE PRECISION A(MDIM,N),U(M) C C HOUSEHOLDER REDUCTION OF RECTANGULAR MATRIX TO UPPER C TRIANGULAR FORM. USE WITH HOLVE FOR LEAST SQUARES C SOLUTIONS OF OVERDETERMINED SYSTEMS. C C MDIM = DECLARED ROW DIMENSION OF A C M = NUMBER OF ROWS OF A C N = NUMBER OF COLUMNS OF A C A = M-BY-N MATRIX WITH M.GE.N C INPUT, MATRIX TO BE REDUCED C OUTPUT, REDUCED MATRIX AND INFORMATION ABOUT REDUCTION C U = M-VECTOR C INPUT, IGNORED C OUTPUT, INFORMATION ABOUT REDUCTION C DOUBLE PRECISION ALPHA,BETA,GAMMA,DSQRT,DSIGN C DO 6 K = 1, N C C FIND REFLECTION WHICH ZEROS A(I,K), I=K+1,...,M C ALPHA = 0.0D0 DO 1 I = K, M U(I) = A(I,K) ALPHA = ALPHA + U(I)**2 1 CONTINUE ALPHA=DSIGN(DSQRT(ALPHA),U(K)) U(K) = U(K) + ALPHA BETA = ALPHA * U(K) A(K,K) = -ALPHA IF (BETA.EQ.0.0D0 .OR. K.EQ.N) GO TO 6 C C APPLY REFLECTION TO REMAINING COLUMNS OF A C KP1 = K + 1 DO 4 J = KP1, N GAMMA = 0.0D0 DO 2 I = K, M GAMMA = GAMMA + U(I) * A(I,J) 2 CONTINUE GAMMA = GAMMA / BETA DO 3 I = K, M A(I,J) = A(I,J) - GAMMA * U(I) 3 CONTINUE 4 CONTINUE 6 CONTINUE RETURN C C TRIANGULAR RESULT STORED IN A(I,J), I.LE.J C VECTORS DEFINING REFLECTIONS STORED IN U AND REST OF A C END SUBROUTINE HOLVE(MDIM,M,N,A,U,B) INTEGER MDIM,M,N DOUBLE PRECISION A(MDIM,N),U(M),B(M) C C LEAST SQUARES SOLUTION OF OVERDETERMINED SYSTEM C FIND X THAT MINIMIZES NORM(A*X-B) C C MDIM,M,N,A,U RESULTS FROM HECOMP C B = M-VECTOR C INPUT, RIGHT HAND SIDE C OUTPUT, FIRST N COMPONENTS = THE SOLUTION, X C LAST M-N COMPONENTS = TRANSFORMED WEIGHTED RESIDUAL C DOUBLE PRECISION BETA,GAMMA,T C C APPLY REFLECTIONS TO B C DO 8 K = 1, N T = A(K,K) BETA = -U(K) * A(K,K) A(K,K) = U(K) GAMMA = 0.0D0 DO 6 I = K, M GAMMA = GAMMA + A(I,K) * B(I) 6 CONTINUE GAMMA = GAMMA / BETA DO 7 I = K, M B(I) = B(I) - GAMMA * A(I,K) 7 CONTINUE A(K,K) = T 8 CONTINUE C C BACK SUBSTITUTION C DO 11 KB = 1, N K = N + 1 - KB B(K) = B(K) / A(K,K) IF (K.EQ.1) GO TO 11 KM1 = K - 1 DO 10 I = 1, KM1 B(I) = B(I) - A(I,K) * B(K) 10 CONTINUE 11 CONTINUE RETURN END