home *** CD-ROM | disk | FTP | other *** search
- 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
