SUBROUTINE SVA (A,MDA,M,N,MDATA,B,SING,NAMES,ISCALE,D) SVA00100 C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 SVA00200 C TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974SVA00300 C SINGULAR VALUE ANALYSIS PRINTOUT SVA00400 C SVA00500 C ISCALE SET BY USER TO 1, 2, OR 3 TO SELECT COLUMN SCALING SVA00600 C OPTION. SVA00700 C 1 SUBR WILL USE IDENTITY SCALING AND IGNORE THE D()SVA00800 C ARRAY. SVA00900 C 2 SUBR WILL SCALE NONZERO COLS TO HAVE UNIT EUCLIDESVA01000 C LENGTH AND WILL STORE RECIPROCAL LENGTHS OF SVA01100 C ORIGINAL NONZERO COLS IN D(). SVA01200 C 3 USER SUPPLIES COL SCALE FACTORS IN D(). SUBR SVA01300 C WILL MULT COL J BY D(J) AND REMOVE THE SCALING SVA01400 C FROM THE SOLN AT THE END. SVA01500 C SVA01600 DIMENSION A(MDA,N), B(M), SING(01),D(N) SVA01700 C SING(3*N) SVA01800 DIMENSION NAMES(N) SVA01900 LOGICAL TEST SVA02000 DOUBLE PRECISION SB, DZERO SVA02100 DZERO=0.D0 SVA02200 ONE=1. SVA02300 ZERO=0. SVA02400 IF (M.LE.0 .OR. N.LE.0) RETURN SVA02500 NP1=N+1 SVA02600 WRITE (6,270) SVA02700 WRITE (6,260) M,N,MDATA SVA02800 GO TO (60,10,10), ISCALE SVA02900 C SVA03000 C APPLY COLUMN SCALING TO A SVA03100 C SVA03200 10 DO 50 J=1,N SVA03300 A1=D(J) SVA03400 GO TO (20,20,40), ISCALE SVA03500 20 SB=DZERO SVA03600 DO 30 I=1,M SVA03700 30 SB=SB+DBLE(A(I,J))**2 SVA03800 A1=DSQRT(SB) SVA03900 IF (A1.EQ.ZERO) A1=ONE SVA04000 A1=ONE/A1 SVA04100 D(J)=A1 SVA04200 40 DO 50 I=1,M SVA04300 50 A(I,J)=A(I,J)*A1 SVA04400 WRITE (6,280) ISCALE,(D(J),J=1,N) SVA04500 GO TO 70 SVA04600 60 CONTINUE SVA04700 WRITE (6,290) SVA04800 70 CONTINUE SVA04900 C SVA05000 C OBTAIN SING. VALUE DECOMP. OF SCALED MATRIX. SVA05100 C SVA05200 C********************************************************** SVA05300 CALL SVDRS (A,MDA,M,N,B,1,1,SING) SVA05400 C********************************************************** SVA05500 C SVA05600 C PRINT THE V MATRIX. SVA05700 C SVA05800 CALL MFEOUT (A,MDA,N,N,NAMES,1) SVA05900 IF (ISCALE.EQ.1) GO TO 90 SVA06000 C SVA06100 C REPLACE V BY D*V IN THE ARRAY A() SVA06200 C SVA06300 DO 80 I=1,N SVA06400 DO 80 J=1,N SVA06500 80 A(I,J)=D(I)*A(I,J) SVA06600 C SVA06700 C G NOW IN B ARRAY. V NOW IN A ARRAY. SVA06800 C SVA06900 C SVA07000 C OBTAIN SUMMARY OUTPUT. SVA07100 C SVA07200 90 CONTINUE SVA07300 WRITE (6,220) SVA07400 C SVA07500 C COMPUTE CUMULATIVE SUMS OF SQUARES OF COMPONENTS OF SVA07600 C G AND STORE THEM IN SING(I), I=MINMN+1,...,2*MINMN+1 SVA07700 C SVA07800 SB=DZERO SVA07900 MINMN=MIN0(M,N) SVA08000 MINMN1=MINMN+1 SVA08100 IF (M.EQ.MINMN) GO TO 110 SVA08200 DO 100 I=MINMN1,M SVA08300 100 SB=SB+DBLE(B(I))**2 SVA08400 110 SING(2*MINMN+1)=SB SVA08500 DO 120 JJ=1,MINMN SVA08600 J=MINMN+1-JJ SVA08700 SB=SB+DBLE(B(J))**2 SVA08800 JS=MINMN+J SVA08900 120 SING(JS)=SB SVA09000 A3=SING(MINMN+1) SVA09100 A4=SQRT(A3/FLOAT(MAX0(1,MDATA))) SVA09200 WRITE (6,230) A3,A4 SVA09300 C SVA09400 NSOL=0 SVA09500 C SVA09600 C SVA09700 C SVA09800 DO 160 K=1,MINMN SVA09900 IF (SING(K).EQ.ZERO) GO TO 130 SVA10000 NSOL=K SVA10100 PI=B(K)/SING(K) SVA10200 A1=ONE/SING(K) SVA10300 A2=B(K)**2 SVA10400 A3=SING(MINMN1+K) SVA10500 A4=SQRT(A3/FLOAT(MAX0(1,MDATA-K))) SVA10600 TEST=SING(K).GE.100..OR.SING(K).LT..001 SVA10700 IF (TEST) WRITE (6,240) K,SING(K),PI,A1,B(K),A2,A3,A4 SVA10800 IF (.NOT.TEST) WRITE (6,250) K,SING(K),PI,A1,B(K),A2,A3,A4 SVA10900 GO TO 140 SVA11000 130 WRITE (6,240) K,SING(K) SVA11100 PI=ZERO SVA11200 140 DO 150 I=1,N SVA11300 A(I,K)=A(I,K)*PI SVA11400 150 IF (K.GT.1) A(I,K)=A(I,K)+A(I,K-1) SVA11500 160 CONTINUE SVA11600 C SVA11700 C COMPUTE AND PRINT VALUES OF YNORM AND RNORM. SVA11800 C SVA11900 WRITE (6,300) SVA12000 J=0 SVA12100 YSQ=ZERO SVA12200 GO TO 180 SVA12300 170 J=J+1 SVA12400 YSQ=YSQ+(B(J)/SING(J))**2 SVA12500 180 YNORM=SQRT(YSQ) SVA12600 JS=MINMN1+J SVA12700 RNORM=SQRT(SING(JS)) SVA12800 YL=-1000. SVA12900 IF (YNORM.GT.0.) YL=ALOG10(YNORM) SVA13000 RL=-1000. SVA13100 IF (RNORM.GT.0.) RL=ALOG10(RNORM) SVA13200 WRITE (6,310) J,YNORM,RNORM,YL,RL SVA13300 IF (J.LT.NSOL) GO TO 170 SVA13400 C SVA13500 C COMPUTE VALUES OF XNORM AND RNORM FOR A SEQUENCE OF VALUES OF SVA13600 C THE LEVENBERG-MARQUARDT PARAMETER. SVA13700 C SVA13800 IF (SING(1).EQ.ZERO) GO TO 210 SVA13900 EL=ALOG10(SING(1))+ONE SVA14000 EL2=ALOG10(SING(NSOL))-ONE SVA14100 DEL=(EL2-EL)/20. SVA14200 TEN=10. SVA14300 ALN10=ALOG(TEN) SVA14400 WRITE (6,320) SVA14500 DO 200 IE=1,21 SVA14600 C COMPUTE ALAMB=10.**EL SVA14700 ALAMB=EXP(ALN10*EL) SVA14800 YS=0. SVA14900 JS=MINMN1+NSOL SVA15000 RS=SING(JS) SVA15100 DO 190 I=1,MINMN SVA15200 SL=SING(I)**2+ALAMB**2 SVA15300 YS=YS+(B(I)*SING(I)/SL)**2 SVA15400 RS=RS+(B(I)*(ALAMB**2)/SL)**2 SVA15500 190 CONTINUE SVA15600 YNORM=SQRT(YS) SVA15700 RNORM=SQRT(RS) SVA15800 RL=-1000. SVA15900 IF (RNORM.GT.ZERO) RL=ALOG10(RNORM) SVA16000 YL=-1000. SVA16100 IF (YNORM.GT.ZERO) YL=ALOG10(YNORM) SVA16200 WRITE (6,330) ALAMB,YNORM,RNORM,EL,YL,RL SVA16300 EL=EL+DEL SVA16400 200 CONTINUE SVA16500 C SVA16600 C PRINT CANDIDATE SOLUTIONS. SVA16700 C SVA16800 210 IF (NSOL.GE.1) CALL MFEOUT (A,MDA,N,NSOL,NAMES,2) SVA16900 RETURN SVA17000 220 FORMAT (42H0 INDEX SING. VALUE P COEF ,48HRECIP. S.SVA17100 1V. G COEF G**2 ,39H C.S.S. SVA17200 2 N.S.R.C.S.S.) SVA17300 230 FORMAT (1H ,5X,1H0,88X,1PE15.4,1PE17.4) SVA17400 240 FORMAT (1H ,I6,E12.4,1P(E15.4,4X,E15.4,4X,E15.4,4X,E15.4,4X,E15.4,SVA17500 12X,E15.4)) SVA17600 250 FORMAT (1H ,I6,F12.4,1P(E15.4,4X,E15.4,4X,E15.4,4X,E15.4,4X,E15.4,SVA17700 12X,E15.4)) SVA17800 260 FORMAT (5H0M = ,I6,8H, N = ,I4,12H, MDATA = ,I8) SVA17900 270 FORMAT (45H0SINGULAR VALUE ANALYSIS OF THE LEAST SQUARES,42H PROBLSVA18000 1EM, A*X=B, SCALED AS (A*D)*Y=B .) SVA18100 280 FORMAT (19H0SCALING OPTION NO.,I2,18H. D IS A DIAGONAL,46H MATRIXSVA18200 1 WITH THE FOLLOWING DIAGONAL ELEMENTS../(5X,10E12.4)) SVA18300 290 FORMAT (50H0SCALING OPTION NO. 1. D IS THE IDENTITY MATRIX./1X) SVA18400 300 FORMAT (6H0INDEX,13X,28H YNORM RNORM,14X,28H LOG1SVA18500 10(YNORM) LOG10(RNORM)/1X) SVA18600 310 FORMAT (1H ,I4,14X,2E14.5,14X,2F14.5) SVA18700 320 FORMAT (54H0NORMS OF SOLUTION AND RESIDUAL VECTORS FOR A RANGE OF,SVA18800 144H VALUES OF THE LEVENBERG-MARQUARDT PARAMETER,9H, LAMBDA./1H0,4XSVA18900 2,42H LAMBDA YNORM RNORM,42H LOG10(LAMBDA) SVA19000 3LOG10(YNORM) LOG10(RNORM)) SVA19100 330 FORMAT (5X,3E14.5,3F14.5) SVA19200 END SVA19300