C TEST ALGS. 1 & 2 J.C. NASH JULY 1978 C USES FRANK MATRIX COLUMNS LOGICAL ESVD,NTOL INTEGER N,ND,IPOS(10),NVAR,MD,I,J,K,YPOS,M REAL A(30,10),D(30,11),G(30),X(10),Z(10),Y(30),Q,V(10,10),EPS EXTERNAL FRANKM C I/O CHANNELS NIN=15 NOUT=16 ND=10 MD=30 ND1=ND+1 1 READ(NIN,900)M,NVAR 900 FORMAT(10I4) WRITE(NOUT,950)M,NVAR 950 FORMAT('0TESTS USING DATA MATRIX',I4,' BY',I4) IF(M.LE.0.OR.NVAR.LE.0)STOP CALL A3PREP(M,NVAR,D,MD,FRANKM) WRITE(NOUT,952) 952 FORMAT('0D MATRIX') CALL OUT(D,MD,M,NVAR,NOUT) 11 READ(NIN,900)IPOS WRITE(NOUT,951)IPOS 951 FORMAT('0COL. #S OF INDEPENDENT VARIABLES'/10I4) N=0 20 N=N+1 IF(IPOS(N).LE.0)GOTO 30 K=IPOS(N) DO 25 J=1,M A(J,N)=D(J,K) 25 CONTINUE GOTO 20 30 N=N-1 IF(N.EQ.0)GOTO 1 ESVD=.FALSE. WRITE(NOUT,953) 953 FORMAT('0A MATRIX') CALL OUT(A,MD,M,N,NOUT) 35 READ(NIN,900)YPOS WRITE(NOUT,954)YPOS 954 FORMAT('0DEPENDENT VARIABLE = COL.',I4) IF(YPOS.LE.0)GOTO 11 DO 40 I=1,M Y(I)=D(I,YPOS) 40 CONTINUE WRITE(NOUT,955)(Y(I),I=1,M) 955 FORMAT(1H ,5E16.8) NTOL=.FALSE. 50 READ(NIN,902)Q 902 FORMAT(E16.8) WRITE(NOUT,956)Q 956 FORMAT('0SING. VALS. .LE.',E16.8,' ARE PRESUMED ZERO') IF(Q.LT.0.0)GOTO 35 C IBM MACHINE PRECISION VALUE EPS=16.0**(-5) CALL A2LSVD(M,N,A,MD,EPS,V,ND,Z,NOUT,Y,G,X,Q,ESVD,NTOL) WRITE(NOUT,957)(J,X(J),J=1,N) 957 FORMAT(' X(',I3,')=',1PE16.8) GOTO 50 END SUBROUTINE OUT(A,NA,N,NP,NOUT) INTEGER NA,N,NOUT,I,J REAL A(NA,NP) DO 20 I=1,N WRITE(NOUT,951)I 951 FORMAT(' ROW',I3) WRITE(NOUT,952)(A(I,J),J=1,NP) 952 FORMAT(1H ,1P5E16.8) 20 CONTINUE RETURN END SUBROUTINE FRANKM(M,N,A,NA) INTEGER M,N,NA,I,J C INPUTS FRANK MATRIX M BY N INTO A REAL A(NA,N) DO 20 I=1,M DO 10 J=1,N A(I,J)=AMIN0(I,J) 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE A3PREP(M,N1,A,NA,AIN) C PREPARE A3 TEST C MATRIX M BY N=N1-1 IS INPUT VIA SUBROUTINE AIN C COL. N1 IS SET TO SUM OF OTHER COLS. - UNIT SOLUTION ELEMENTS C BUT ONLY IF M=N - OTHERWISE SIMPLY INPUT MATRIX C NA = FIRST DIMENSION OF A INTEGER M,N1,NA,N,J,I REAL A(NA,N1),S N=N1-1 CALL AIN(M,N,A,NA) IF(M.NE.N)RETURN DO 40 I=1,N S=0.0 DO 30 J=1,N S=S+A(I,J) 30 CONTINUE A(I,N1)=S 40 CONTINUE RETURN END SUBROUTINE A1SVD(M,N,A,NA,EPS,V,NV,Z,IPR) C ALGORITHM 1 SINGULAR VALUE DECOMPOSITION BY COLUMN ORTHOGONA- C LISATION VIA PLANE ROTATIONS C J.C.NASH FEBRUARY 1980 C M BY N MATRIX A IS DECOMPOSED TO U*Z*VT C A = ARRAY CONTAINING A (INPUT), U (OUTPUT) C NA = FIRST DIMENSION OF A C EPS = MACHINE PRECISION C V = ARRAY IN WHICH ORTHOGAONAL MATRIX V IS ACCUMULATED C NV = FIRST DIMENSION OF V C Z = VECTOR OF SINGULAR VALUES C IPR = PRINT CHANNEL IF IPR.GT.0 THEN PRINTING C STEP 0 INTEGER M,N,J1,N1,COUNT REAL A(NA,N),V(NV,N),Z(N),EPS,TOL,P,Q,R,VV,C,S C UNDERFLOW AVOIDANCE STRATEGY REAL SMALL SMALL=1.0E-36 C ABOVE IS VALUE FOR IBM TOL=N*N*EPS*EPS DO 6 I=1,N DO 4 J=1,N V(I,J)=0.0 4 CONTINUE V(I,I)=1.0 6 CONTINUE N1=N-1 C STEP 1 10 COUNT=N*(N-1)/2 C STEP 2 DO 140 J=1,N1 C STEP 3 J1=J+1 DO 130 K=J1,N C STEP 4 P=0.0 Q=0.0 R=0.0 C STEP 5 DO 55 I=1,M IF(ABS(A(I,J)).GT.SMALL.AND.ABS(A(I,K)).GT.SMALL) # P=P+A(I,J)*A(I,K) IF(ABS(A(I,J)).GT.SMALL)Q=Q+A(I,J)**2 IF(ABS(A(I,K)).GT.SMALL)R=R+A(I,K)**2 C P=P+A(I,J)*A(I,K) C Q=Q+A(I,J)**2 C R=R+A(I,K)**2 55 CONTINUE C STEP 6 IF(Q.GE.R)GOTO 70 C=0.0 S=1.0 GOTO 90 C STEP 7 70 IF(R.LE.TOL)GOTO 120 IF((P*P)/(Q*R).LT.TOL)GOTO 120 C STEP 8 Q=Q-R VV=SQRT(4.0*P**2+Q**2) C=SQRT((VV+Q)/(2.0*VV)) S=P/(VV*C) C CHANGE SYMBOLS VV TO W ????????????????????????????????????????????? C STEP 9 90 DO 95 I=1,M R=A(I,J) A(I,J)=R*C+A(I,K)*S A(I,K)=-R*S+A(I,K)*C 95 CONTINUE C STEP 10 DO 105 I=1,N R=V(I,J) V(I,J)=R*C+V(I,K)*S V(I,K)=-R*S+V(I,K)*C 105 CONTINUE C STEP 11 GOTO 130 120 COUNT=COUNT-1 C STEP 13 130 CONTINUE C STEP 14 140 CONTINUE C STEP 15 IF(IPR.GT.0)WRITE(IPR,964)COUNT 964 FORMAT(1H ,I4,10H ROTATIONS) IF(COUNT.GT.0)GOTO 10 C STEP 16 DO 220 J=1,N C STEP 17 Q=0.0 C STEP 18 DO 185 I=1,M Q=Q+A(I,J)**2 185 CONTINUE C STEP 19 Q=SQRT(Q) Z(J)=Q IF(IPR.GT.0)WRITE(IPR,965)J,Q 965 FORMAT( 4H SV(,I3,2H)=,1PE16.8) C STEP 20 IF(Q.LT.TOL)GOTO 220 C STEP 21 DO 215 I=1,M A(I,J)=A(I,J)/Q 215 CONTINUE C STEP 22 220 CONTINUE RETURN END SUBROUTINE A2LSVD(M,N,A,NA,EPS,V,NV,Z,IPR,Y,G,X,Q,ESVD,NTOL) C J.C.NASH FEBRUARY 1980 C SAME COMMENTS AS SUBN A1SVD EXCEPT FOR C G = WORKING VECTOR IN N ELEMENTS C Y = VECTOR CONTAINING M VALUES OF DEPENDENT VARIABLE C X = SOLUTION VECTOR C Q = TOLERANCE FOR SINGULAR VALUES. THOSE .LE. Q TAKEN AS ZERO. C ESVD = LOGICAL FLAG SET .TRUE. IF SVD ALREADY EXISTS IN A,Z,V C NTOL = LOGICAL FLAG SET .TRUE. IF ONLY NEW TOLERANCE Q. LOGICAL ESVD,NTOL INTEGER M,N,IPR,I,J REAL A(NA,N),V(NV,N),Z(N),Y(M),G(N),X(N),S,Q C STEP 1 IF(NTOL)GOTO 41 IF(.NOT.ESVD)CALL A1SVD(M,N,A,NA,EPS,V,NV,Z,IPR) IF(IPR.GT.0)WRITE(IPR,965)(J,Z(J),J=1,N) 965 FORMAT(16H SINGULAR VALUE(,I3,2H)=,1PE16.8) C STEP 2 VIA SUBROUTINE CALL C ALTERNATIVE WITHOUT G C NO STEP 3 C STEP 3 UT*Y=G DO 36 J=1,N S=0.0 DO 34 I=1,M S=S+A(I,J)*Y(I) 34 CONTINUE G(J)=S 36 CONTINUE C STEP 4 41 IF(Q.LT.0.0)STOP C STEP 5 DO 56 J=1,N S=0.0 DO 54 I=1,N IF(Z(I).GT.Q)S=S+V(J,I)*G(I)/Z(I) 54 CONTINUE X(J)=S 56 CONTINUE C STEP 6 C NEW TOLERANCE VIA NEW CALL RETURN END 4 5 1 2 0 3 0.0 E 00 -1.0 E 00 4 0.0 E 00 -1.0 E 00 5 0.0 E 00 -1.0 E 00 -1 1 2 3 4 5 0 1 0.0 E 00 1.0 E-04 -1.0 E 00 -1 -1 20 10 1 2 3 0 4 0.0 E 00 1.0 E-04 -1.0 E 00 -1 -1 0