SUBROUTINE SVA (A,MDA,M,N,MDATA,B,SING,NAMES,ISCALE,D) C C INPUT.. C C A(*,*),MDA,M,N THE ARRAY A(*,*) INITIALLY CONTAINS THE M BY C N MATRIX A OF THE LEAST SQUARES PROBLEM AX = C B. EITHER M.GE.N OR M.LT.N IS PERMITTED. THE C FIRST DIMENSIONING PARAMETER OF THE ARRAY C A(*,*) IS MDA, WHICH MUST SATISFY MDA .GE. C MAX(M,N). THE CONDITION MDA.LT.MAX(M,N) IS C CONSIDERED AN ERROR. C C MDATA THE NUMBER OF ROWS IN THE ORIGINAL LEAST C SQUARES PROBLEM, WHICH MAY HAVE HAD MORE ROWS C THAN (A B). THIS NUMBER IS USED ONLY IN C COMPUTING THE NUMBERS SG(J), J=0,...,K. C C B(*) THE ARRAY B(*) INITIALLY CONTAINS THE M-VECTOR C B OF THE LEAST SQUARES PROBLEM AX = B. C C SING(*) THE LOCATIONS, (SING(I), I=1,...,3*N), ARE C USED AS TEMPORARY WORKING SPACE. C C NAMES(*) IN LOCATION NAMES(I) THE USER MAY STORE AN C ALPHA-NUMERIC NAME AS AN IDENTIFIER FOR THE C I-TH COMPONENT OF THE SOLUTION VECTOR, (I = C 1,...,N). THESE NAMES WILL BE PRINTED WITH C THE APPROPRIATE COMPONENTS OF THE V MATRIX C AND WITH THE CANDIDATE SOLUTIONS. EACH NAME C IS PRINTED WITH AN A6 FORMAT. IF NAMES(1) HAS C THE SAME CONTENTS AS THE QUANTITY IBLANK, C DEFINED BY C C DATA IBLANK/' '/ C C THEN THE REMAINING WORDS OF THE ARRAY NAMES(*) C WILL BE IGNORED AND NO IDENTIFYING NAMES WILL C BE PRINTED. C C ISCALE,D(*) THE USER SETS ISCALE = 1,2, OR 3. C IF ISCALE = 1, THE SUBROUTINE FUNCTIONS AS C THOUGH THE SCALING MATRIX D WERE AN N BY N C IDENTITY MATRIX. IN THIS CASE THE SUBROUTINE C MAKES NO REFERENCE TO THE ARRAY D(*). C C IF ISCALE = 2, THE SUBROUTINE COMPUTES THE C LENGTH OF THE J-TH COLUMN VECTOR OF A AND SETS C D(J) EQUAL TO THE RECIPROCAL OF THIS LENGTH, C IF NONZERO, AND D(J) = 1, IF ZERO. IN THIS C CASE ALL NONZERO COLUMNS OF THE SCALED MATRIX C C = AD WILL HAVE UNIT EUCLIDEAN LENGTH. C C IF ISCALE = 3, THE SUBROUTINE WILL USE THE C USER-SUPPLIED VALUES D(J),J = 1,...,N AS THE C DIAGONAL ELEMENTS OF THE DIAGONAL SCALING C MATRIX D. IN THIS CASE THE USER MUST ASSIGN C VALUES TO D(J), J = 1,...,N, AND THE C SUBROUTINE WILL NOT MODIFY THESE VALUES. FOR C EXAMPLE, THE USER MIGHT SET D(J) EQUAL TO THE C A PRIORI STANDARD DEVIATION OF THE J-TH C COMPONENT OF THE SOLUTION VECTOR IF SUCH C INFORMATION IS KNOWN. AS A FURTHER EXAMPLE, C THE USER MIGHT SET D(J) EQUAL TO THE C RECIPROCAL OF THE MEAN OF THE M COMPONENTS OF C THE J-TH COLUMN VECTOR S UNCERTAINTY. THIS C TYPE OF SCALING HAS THE EFFECT OF MAKING THE C UNCERTAINTY OF EVERY COLUMN THE SAME MAGNITUDE C IN THE MATRIX C. C C OUTPUT.. C C A(*,*) ON OUTPUT, THE J-TH CANDIDATE SOLUTION WILL BE C STORED IN THE FIRST N LOCATIONS OF THE J-TH C COLUMN OF A(*,*) FOR J = 1,...,MIN(M,N). C C B(*) ON OUTPUT THE M-VECTOR G IS STORED IN THE C ARRAY B(*). C C SING(*) ON OUTPUT THE SINGULAR VALUES OF THE SCALED C MATRIX C ARE STORED IN SING(I),I = 1,..., C MIN(M,N). C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) C DIMENSION A(MDA,N), B(M), SING(01),D(N) C SING(3*N) CHARACTER*6 NAMES(N) LOGICAL TEST C DZERO=0.D0 ONE=1. ZERO=0. IF (M.LE.0 .OR. N.LE.0) RETURN IF(MAX(M,N).GT.MDA) RETURN NP1=N+1 WRITE (6,270) WRITE (6,260) M,N,MDATA GO TO (60,10,10), ISCALE C C APPLY COLUMN SCALING TO A C 10 DO 50 J=1,N A1=D(J) GO TO (20,20,40), ISCALE 20 SB=DZERO DO 30 I=1,M 30 SB=SB+A(I,J)**2 A1=SQRT(SB) IF (A1.EQ.ZERO) A1=ONE A1=ONE/A1 D(J)=A1 40 DO 50 I=1,M 50 A(I,J)=A(I,J)*A1 WRITE (6,280) ISCALE,(D(J),J=1,N) GO TO 70 60 CONTINUE WRITE (6,290) 70 CONTINUE C C OBTAIN SING. VALUE DECOMP. OF SCALED MATRIX. C C********************************************************** CALL SVDRS (A,MDA,M,N,B,1,1,SING) C********************************************************** C C PRINT THE V MATRIX. C CALL MFEOUT (A,MDA,N,N,NAMES,1) IF (ISCALE.EQ.1) GO TO 90 C C REPLACE V BY D*V IN THE ARRAY A() C DO 80 I=1,N DO 80 J=1,N 80 A(I,J)=D(I)*A(I,J) C C G NOW IN B ARRAY. V NOW IN A ARRAY. C C C OBTAIN SUMMARY OUTPUT. C 90 CONTINUE WRITE (6,220) C C COMPUTE CUMULATIVE SUMS OF SQUARES OF COMPONENTS OF C G AND STORE THEM IN SING(I), I=MINMN+1,...,2*MINMN+1 C SB=DZERO MINMN=MIN(M,N) MINMN1=MINMN+1 IF (M.EQ.MINMN) GO TO 110 DO 100 I=MINMN1,M 100 SB=SB+B(I)**2 110 SING(2*MINMN+1)=SB DO 120 JJ=1,MINMN J=MINMN+1-JJ SB=SB+B(J)**2 JS=MINMN+J 120 SING(JS)=SB A3=SING(MINMN+1) A4=SQRT(A3/DBLE(MAX(1,MDATA))) WRITE (6,230) A3,A4 C NSOL=0 C C C DO 160 K=1,MINMN IF (SING(K).EQ.ZERO) GO TO 130 NSOL=K PI=B(K)/SING(K) A1=ONE/SING(K) A2=B(K)**2 A3=SING(MINMN1+K) A4=SQRT(A3/DBLE(MAX(1,MDATA-K))) TEST=SING(K).GE.100..OR.SING(K).LT..001 IF (TEST) WRITE (6,240) K,SING(K),PI,A1,B(K),A2,A3,A4 IF (.NOT.TEST) WRITE(6,250)K,SING(K),PI,A1,B(K),A2,A3,A4 GO TO 140 130 WRITE (6,240) K,SING(K) PI=ZERO 140 DO 150 I=1,N A(I,K)=A(I,K)*PI 150 IF (K.GT.1) A(I,K)=A(I,K)+A(I,K-1) 160 CONTINUE C C COMPUTE AND PRINT VALUES OF YNORM AND RNORM. C WRITE (6,300) J=0 YSQ=ZERO GO TO 180 170 J=J+1 YSQ=YSQ+(B(J)/SING(J))**2 180 YNORM=SQRT(YSQ) JS=MINMN1+J RNORM=SQRT(SING(JS)) YL=-1000. IF (YNORM.GT.0.) YL=LOG10(YNORM) RL=-1000. IF (RNORM.GT.0.) RL=LOG10(RNORM) WRITE (6,310) J,YNORM,RNORM,YL,RL IF (J.LT.NSOL) GO TO 170 C C COMPUTE VALUES OF XNORM AND RNORM FOR A SEQUENCE OF VALUES OF C THE LEVENBERG-MARQUARDT PARAMETER. C IF (SING(1).EQ.ZERO) GO TO 210 EL=LOG10(SING(1))+ONE EL2=LOG10(SING(NSOL))-ONE DEL=(EL2-EL)/20. TEN=10. ALN10=LOG(TEN) WRITE (6,320) DO 200 IE=1,21 C COMPUTE ALAMB=10.**EL ALAMB=EXP(ALN10*EL) YS=0. JS=MINMN1+NSOL RS=SING(JS) DO 190 I=1,MINMN SL=SING(I)**2+ALAMB**2 YS=YS+(B(I)*SING(I)/SL)**2 RS=RS+(B(I)*(ALAMB**2)/SL)**2 190 CONTINUE YNORM=SQRT(YS) RNORM=SQRT(RS) RL=-1000. IF (RNORM.GT.ZERO) RL=LOG10(RNORM) YL=-1000. IF (YNORM.GT.ZERO) YL=LOG10(YNORM) WRITE (6,330) ALAMB,YNORM,RNORM,EL,YL,RL EL=EL+DEL 200 CONTINUE C C PRINT CANDIDATE SOLUTIONS. C 210 IF (NSOL.GE.1) CALL MFEOUT (A,MDA,N,NSOL,NAMES,2) RETURN 220 FORMAT (42H0 INDEX SING. VALUE P COEF ,48HRECIP. S. 1V. G COEF G**2 ,39H C.S.S. 2 N.S.R.C.S.S.) 230 FORMAT (1H ,5X,1H0,88X,1PE15.4,1PE17.4) 240 FORMAT (1H ,I6,E12.4,1P(E15.4,4X,E15.4,4X,E15.4,4X,E15.4,4X,E15.4, 12X,E15.4)) 250 FORMAT (1H ,I6,F12.4,1P(E15.4,4X,E15.4,4X,E15.4,4X,E15.4,4X,E15.4, 12X,E15.4)) 260 FORMAT (5H0M = ,I6,8H, N = ,I4,12H, MDATA = ,I8) 270 FORMAT (45H0SINGULAR VALUE ANALYSIS OF THE LEAST SQUARES,42H PROBL 1EM, A*X=B, SCALED AS (A*D)*Y=B .) 280 FORMAT (19H0SCALING OPTION NO.,I2,18H. D IS A DIAGONAL,46H MATRIX 1 WITH THE FOLLOWING DIAGONAL ELEMENTS../(5X,10E12.4)) 290 FORMAT (50H0SCALING OPTION NO. 1. D IS THE IDENTITY MATRIX./1X) 300 FORMAT (6H0INDEX,13X,28H YNORM RNORM,14X,28H LOG1 10(YNORM) LOG10(RNORM)/1X) 310 FORMAT (1H ,I4,14X,2E14.5,14X,2F14.5) 320 FORMAT (54H0NORMS OF SOLUTION AND RESIDUAL VECTORS FOR A RANGE OF, 144H VALUES OF THE LEVENBERG-MARQUARDT PARAMETER,9H, LAMBDA./1H0,4X 2,42H LAMBDA YNORM RNORM,42H LOG10(LAMBDA) 3LOG10(YNORM) LOG10(RNORM)) 330 FORMAT (5X,3E14.5,3F14.5) END