C SUBROUTINE SVDRS (A,MDA,MM,NN,B,MDB,NB,S) SVD00100 C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1974 MAR 1 SVD00200 C TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974SVD00300 C SINGULAR VALUE DECOMPOSITION ALSO TREATING RIGHT SIDE VECTOR.SVD00400 C SVD00500 C THE ARRAY S OCCUPIES 3*N CELLS. SVD00600 C A OCCUPIES M*N CELLS SVD00700 C B OCCUPIES M*NB CELLS. SVD00800 C SVD00900 C SPECIAL SINGULAR VALUE DECOMPOSITION SUBROUTINE. SVD01000 C WE HAVE THE M X N MATRIX A AND THE SYSTEM A*X=B TO SOLVE. SVD01100 C EITHER M .GE. N OR M .LT. N IS PERMITTED. SVD01200 C THE SINGULAR VALUE DECOMPOSITION SVD01300 C A = U*S*V**(T) IS MADE IN SUCH A WAY THAT ONE GETS SVD01400 C (1) THE MATRIX V IN THE FIRST N ROWS AND COLUMNS OF A. SVD01500 C (2) THE DIAGONAL MATRIX OF ORDERED SINGULAR VALUES IN SVD01600 C THE FIRST N CELLS OF THE ARRAY S(IP), IP .GE. 3*N. SVD01700 C (3) THE MATRIX PRODUCT U**(T)*B=G GETS PLACED BACK IN B. SVD01800 C (4) THE USER MUST COMPLETE THE SOLUTION AND DO HIS OWN SVD01900 C SINGULAR VALUE ANALYSIS. SVD02000 C ******* SVD02100 C GIVE SPECIAL SVD02200 C TREATMENT TO ROWS AND COLUMNS WHICH ARE ENTIRELY ZERO. THIS SVD02300 C CAUSES CERTAIN ZERO SING. VALS. TO APPEAR AS EXACT ZEROS RATHER SVD02400 C THAN AS ABOUT ETA TIMES THE LARGEST SING. VAL. IT SIMILARLY SVD02500 C CLEANS UP THE ASSOCIATED COLUMNS OF U AND V. SVD02600 C METHOD.. SVD02700 C 1. EXCHANGE COLS OF A TO PACK NONZERO COLS TO THE LEFT. SVD02800 C SET N = NO. OF NONZERO COLS. SVD02900 C USE LOCATIONS A(1,NN),A(1,NN-1),...,A(1,N+1) TO RECORD THE SVD03000 C COL PERMUTATIONS. SVD03100 C 2. EXCHANGE ROWS OF A TO PACK NONZERO ROWS TO THE TOP. SVD03200 C QUIT PACKING IF FIND N NONZERO ROWS. MAKE SAME ROW EXCHANGES SVD03300 C IN B. SET M SO THAT ALL NONZERO ROWS OF THE PERMUTED A SVD03400 C ARE IN FIRST M ROWS. IF M .LE. N THEN ALL M ROWS ARE SVD03500 C NONZERO. IF M .GT. N THEN THE FIRST N ROWS ARE KNOWN SVD03600 C TO BE NONZERO,AND ROWS N+1 THRU M MAY BE ZERO OR NONZERO. SVD03700 C 3. APPLY ORIGINAL ALGORITHM TO THE M BY N PROBLEM. SVD03800 C 4. MOVE PERMUTATION RECORD FROM A(,) TO S(I),I=N+1,...,NN. SVD03900 C 5. BUILD V UP FROM N BY N TO NN BY NN BY PLACING ONES ON SVD04000 C THE DIAGONAL AND ZEROS ELSEWHERE. THIS IS ONLY PARTLY DONE SVD04100 C EXPLICITLY. IT IS COMPLETED DURING STEP 6. SVD04200 C 6. EXCHANGE ROWS OF V TO COMPENSATE FOR COL EXCHANGES OF STEP 2. SVD04300 C 7. PLACE ZEROS IN S(I),I=N+1,NN TO REPRESENT ZERO SING VALS. SVD04400 C SVD04500 SUBROUTINE SVDRS (A,MDA,MM,NN,B,MDB,NB,S) SVD04600 DIMENSION A(MDA,NN),B(MDB,NB),S(NN,3) SVD04700 ZERO=0. SVD04800 ONE=1. SVD04900 C SVD05000 C BEGIN.. SPECIAL FOR ZERO ROWS AND COLS. SVD05100 C SVD05200 C PACK THE NONZERO COLS TO THE LEFT SVD05300 C SVD05400 N=NN SVD05500 IF (N.LE.0.OR.MM.LE.0) RETURN SVD05600 J=N SVD05700 10 CONTINUE SVD05800 DO 20 I=1,MM SVD05900 IF (A(I,J)) 50,20,50 SVD06000 20 CONTINUE SVD06100 C SVD06200 C COL J IS ZERO. EXCHANGE IT WITH COL N. SVD06300 C SVD06400 IF (J.EQ.N) GO TO 40 SVD06500 DO 30 I=1,MM SVD06600 30 A(I,J)=A(I,N) SVD06700 40 CONTINUE SVD06800 A(1,N)=J SVD06900 N=N-1 SVD07000 50 CONTINUE SVD07100 J=J-1 SVD07200 IF (J.GE.1) GO TO 10 SVD07300 C IF N=0 THEN A IS ENTIRELY ZERO AND SVD SVD07400 C COMPUTATION CAN BE SKIPPED SVD07500 NS=0 SVD07600 IF (N.EQ.0) GO TO 240 SVD07700 C PACK NONZERO ROWS TO THE TOP SVD07800 C QUIT PACKING IF FIND N NONZERO ROWS SVD07900 I=1 SVD08000 M=MM SVD08100 60 IF (I.GT.N.OR.I.GE.M) GO TO 150 SVD08200 IF (A(I,I)) 90,70,90 SVD08300 70 DO 80 J=1,N SVD08400 IF (A(I,J)) 90,80,90 SVD08500 80 CONTINUE SVD08600 GO TO 100 SVD08700 90 I=I+1 SVD08800 GO TO 60 SVD08900 C ROW I IS ZERO SVD09000 C EXCHANGE ROWS I AND M SVD09100 100 IF(NB.LE.0) GO TO 115 SVD09200 DO 110 J=1,NB SVD09300 T=B(I,J) SVD09400 B(I,J)=B(M,J) SVD09500 110 B(M,J)=T SVD09600 115 DO 120 J=1,N SVD09700 120 A(I,J)=A(M,J) SVD09800 IF (M.GT.N) GO TO 140 SVD09900 DO 130 J=1,N SVD10000 130 A(M,J)=ZERO SVD10100 140 CONTINUE SVD10200 C EXCHANGE IS FINISHED SVD10300 M=M-1 SVD10400 GO TO 60 SVD10500 C SVD10600 150 CONTINUE SVD10700 C END.. SPECIAL FOR ZERO ROWS AND COLUMNS SVD10800 C BEGIN.. SVD ALGORITHM SVD10900 C METHOD.. SVD11000 C (1) REDUCE THE MATRIX TO UPPER BIDIAGONAL FORM WITH SVD11100 C HOUSEHOLDER TRANSFORMATIONS. SVD11200 C H(N)...H(1)AQ(1)...Q(N-2) = (D**T,0)**T SVD11300 C WHERE D IS UPPER BIDIAGONAL. SVD11400 C SVD11500 C (2) APPLY H(N)...H(1) TO B. HERE H(N)...H(1)*B REPLACES B SVD11600 C IN STORAGE. SVD11700 C SVD11800 C (3) THE MATRIX PRODUCT W= Q(1)...Q(N-2) OVERWRITES THE FIRST SVD11900 C N ROWS OF A IN STORAGE. SVD12000 C SVD12100 C (4) AN SVD FOR D IS COMPUTED. HERE K ROTATIONS RI AND PI ARE SVD12200 C COMPUTED SO THAT SVD12300 C RK...R1*D*P1**(T)...PK**(T) = DIAG(S1,...,SM) SVD12400 C TO WORKING ACCURACY. THE SI ARE NONNEGATIVE AND NONINCREASING. SVD12500 C HERE RK...R1*B OVERWRITES B IN STORAGE WHILE SVD12600 C A*P1**(T)...PK**(T) OVERWRITES A IN STORAGE. SVD12700 C SVD12800 C (5) IT FOLLOWS THAT,WITH THE PROPER DEFINITIONS, SVD12900 C U**(T)*B OVERWRITES B, WHILE V OVERWRITES THE FIRST N ROW AND SVD13000 C COLUMNS OF A. SVD13100 C SVD13200 L=MIN0(M,N) SVD13300 C THE FOLLOWING LOOP REDUCES A TO UPPER BIDIAGONAL AND SVD13400 C ALSO APPLIES THE PREMULTIPLYING TRANSFORMATIONS TO B. SVD13500 C SVD13600 DO 170 J=1,L SVD13700 IF (J.GE.M) GO TO 160 SVD13800 CALL H12 (1,J,J+1,M,A(1,J),1,T,A(1,J+1),1,MDA,N-J) SVD13900 CALL H12 (2,J,J+1,M,A(1,J),1,T,B,1,MDB,NB) SVD14000 160 IF (J.GE.N-1) GO TO 170 SVD14100 CALL H12 (1,J+1,J+2,N,A(J,1),MDA,S(J,3),A(J+1,1),MDA,1,M-J) SVD14200 170 CONTINUE SVD14300 C SVD14400 C COPY THE BIDIAGONAL MATRIX INTO THE ARRAY S() FOR QRBD. SVD14500 C SVD14600 IF (N.EQ.1) GO TO 190 SVD14700 DO 180 J=2,N SVD14800 S(J,1)=A(J,J) SVD14900 180 S(J,2)=A(J-1,J) SVD15000 190 S(1,1)=A(1,1) SVD15100 C SVD15200 NS=N SVD15300 IF (M.GE.N) GO TO 200 SVD15400 NS=M+1 SVD15500 S(NS,1)=ZERO SVD15600 S(NS,2)=A(M,M+1) SVD15700 200 CONTINUE SVD15800 C SVD15900 C CONSTRUCT THE EXPLICIT N BY N PRODUCT MATRIX, W=Q1*Q2*...*QL*I SVD16000 C IN THE ARRAY A(). SVD16100 C SVD16200 DO 230 K=1,N SVD16300 I=N+1-K SVD16400 IF(I.GT.MIN0(M,N-2)) GO TO 210 SVD16500 CALL H12 (2,I+1,I+2,N,A(I,1),MDA,S(I,3),A(1,I+1),1,MDA,N-I) SVD16600 210 DO 220 J=1,N SVD16700 220 A(I,J)=ZERO SVD16800 230 A(I,I)=ONE SVD16900 C SVD17000 C COMPUTE THE SVD OF THE BIDIAGONAL MATRIX SVD17100 C SVD17200 CALL QRBD (IPASS,S(1,1),S(1,2),NS,A,MDA,N,B,MDB,NB) SVD17300 C SVD17400 GO TO (240,310), IPASS SVD17500 240 CONTINUE SVD17600 IF (NS.GE.N) GO TO 260 SVD17700 NSP1=NS+1 SVD17800 DO 250 J=NSP1,N SVD17900 250 S(J,1)=ZERO SVD18000 260 CONTINUE SVD18100 IF (N.EQ.NN) RETURN SVD18200 NP1=N+1 SVD18300 C MOVE RECORD OF PERMUTATIONS SVD18400 C AND STORE ZEROS SVD18500 DO 280 J=NP1,NN SVD18600 S(J,1)=A(1,J) SVD18700 DO 270 I=1,N SVD18800 270 A(I,J)=ZERO SVD18900 280 CONTINUE SVD19000 C PERMUTE ROWS AND SET ZERO SINGULAR VALUES.SVD19100 DO 300 K=NP1,NN SVD19200 I=S(K,1) SVD19300 S(K,1)=ZERO SVD19400 DO 290 J=1,NN SVD19500 A(K,J)=A(I,J) SVD19600 290 A(I,J)=ZERO SVD19700 A(I,K)=ONE SVD19800 300 CONTINUE SVD19900 C END.. SPECIAL FOR ZERO ROWS AND COLUMNS SVD20000 RETURN SVD20100 310 WRITE (6,320) SVD20200 STOP SVD20300 320 FORMAT (49H CONVERGENCE FAILURE IN QR BIDIAGONAL SVD ROUTINE) SVD20400 END SVD20500 C SUBROUTINE QRBD (IPASS,Q,E,NN,V,MDV,NRV,C,MDC,NCC) QBD00100 C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 QBD00200 C TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974QBD00300 C QR ALGORITHM FOR SINGULAR VALUES OF A BIDIAGONAL MATRIX. QBD00400 C QBD00500 C THE BIDIAGONAL MATRIX QBD00600 C QBD00700 C (Q1,E2,0... ) QBD00800 C ( Q2,E3,0... ) QBD00900 C D= ( . ) QBD01000 C ( . 0) QBD01100 C ( .EN) QBD01200 C ( 0,QN) QBD01300 C QBD01400 C IS PRE AND POST MULTIPLIED BY QBD01500 C ELEMENTARY ROTATION MATRICES QBD01600 C RI AND PI SO THAT QBD01700 C QBD01800 C RK...R1*D*P1**(T)...PK**(T) = DIAG(S1,...,SN) QBD01900 C QBD02000 C TO WITHIN WORKING ACCURACY. QBD02100 C QBD02200 C 1. EI AND QI OCCUPY E(I) AND Q(I) AS INPUT. QBD02300 C QBD02400 C 2. RM...R1*C REPLACES 'C' IN STORAGE AS OUTPUT. QBD02500 C QBD02600 C 3. V*P1**(T)...PM**(T) REPLACES 'V' IN STORAGE AS OUTPUT. QBD02700 C QBD02800 C 4. SI OCCUPIES Q(I) AS OUTPUT. QBD02900 C QBD03000 C 5. THE SI'S ARE NONINCREASING AND NONNEGATIVE. QBD03100 C QBD03200 C THIS CODE IS BASED ON THE PAPER AND 'ALGOL' CODE.. QBD03300 C REF.. QBD03400 C 1. REINSCH,C.H. AND GOLUB,G.H. 'SINGULAR VALUE DECOMPOSITION QBD03500 C AND LEAST SQUARES SOLUTIONS' (NUMER. MATH.), VOL. 14,(1970). QBD03600 C QBD03700 SUBROUTINE QRBD (IPASS,Q,E,NN,V,MDV,NRV,C,MDC,NCC) QBD03800 LOGICAL WNTV ,HAVERS,FAIL QBD03900 DIMENSION Q(NN),E(NN),V(MDV,NN),C(MDC,NCC) QBD04000 ZERO=0. QBD04100 ONE=1. QBD04200 TWO=2. QBD04300 C QBD04400 N=NN QBD04500 IPASS=1 QBD04600 IF (N.LE.0) RETURN QBD04700 N10=10*N QBD04800 WNTV=NRV.GT.0 QBD04900 HAVERS=NCC.GT.0 QBD05000 FAIL=.FALSE. QBD05100 NQRS=0 QBD05200 E(1)=ZERO QBD05300 DNORM=ZERO QBD05400 DO 10 J=1,N QBD05500 10 DNORM=AMAX1(ABS(Q(J))+ABS(E(J)),DNORM) QBD05600 DO 200 KK=1,N QBD05700 K=N+1-KK QBD05800 C QBD05900 C TEST FOR SPLITTING OR RANK DEFICIENCIES.. QBD06000 C FIRST MAKE TEST FOR LAST DIAGONAL TERM, Q(K), BEING SMALL. QBD06100 20 IF(K.EQ.1) GO TO 50 QBD06200 IF(DIFF(DNORM+Q(K),DNORM)) 50,25,50 QBD06300 C QBD06400 C SINCE Q(K) IS SMALL WE WILL MAKE A SPECIAL PASS TO QBD06500 C TRANSFORM E(K) TO ZERO. QBD06600 C QBD06700 25 CS=ZERO QBD06800 SN=-ONE QBD06900 DO 40 II=2,K QBD07000 I=K+1-II QBD07100 F=-SN*E(I+1) QBD07200 E(I+1)=CS*E(I+1) QBD07300 CALL G1 (Q(I),F,CS,SN,Q(I)) QBD07400 C TRANSFORMATION CONSTRUCTED TO ZERO POSITION (I,K). QBD07500 C QBD07600 IF (.NOT.WNTV) GO TO 40 QBD07700 DO 30 J=1,NRV QBD07800 30 CALL G2 (CS,SN,V(J,I),V(J,K)) QBD07900 C ACCUMULATE RT. TRANSFORMATIONS IN V. QBD08000 C QBD08100 40 CONTINUE QBD08200 C QBD08300 C THE MATRIX IS NOW BIDIAGONAL, AND OF LOWER ORDER QBD08400 C SINCE E(K) .EQ. ZERO.. QBD08500 C QBD08600 50 DO 60 LL=1,K QBD08700 L=K+1-LL QBD08800 IF(DIFF(DNORM+E(L),DNORM)) 55,100,55 QBD08900 55 IF(DIFF(DNORM+Q(L-1),DNORM)) 60,70,60 QBD09000 60 CONTINUE QBD09100 C THIS LOOP CAN'T COMPLETE SINCE E(1) = ZERO. QBD09200 C QBD09300 GO TO 100 QBD09400 C QBD09500 C CANCELLATION OF E(L), L.GT.1. QBD09600 70 CS=ZERO QBD09700 SN=-ONE QBD09800 DO 90 I=L,K QBD09900 F=-SN*E(I) QBD10000 E(I)=CS*E(I) QBD10100 IF(DIFF(DNORM+F,DNORM)) 75,100,75 QBD10200 75 CALL G1 (Q(I),F,CS,SN,Q(I)) QBD10300 IF (.NOT.HAVERS) GO TO 90 QBD10400 DO 80 J=1,NCC QBD10500 80 CALL G2 (CS,SN,C(I,J),C(L-1,J)) QBD10600 90 CONTINUE QBD10700 C QBD10800 C TEST FOR CONVERGENCE.. QBD10900 100 Z=Q(K) QBD11000 IF (L.EQ.K) GO TO 170 QBD11100 C QBD11200 C SHIFT FROM BOTTOM 2 BY 2 MINOR OF B**(T)*B. QBD11300 X=Q(L) QBD11400 Y=Q(K-1) QBD11500 G=E(K-1) QBD11600 H=E(K) QBD11700 F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(TWO*H*Y) QBD11800 G=SQRT(ONE+F**2) QBD11900 IF (F.LT.ZERO) GO TO 110 QBD12000 T=F+G QBD12100 GO TO 120 QBD12200 110 T=F-G QBD12300 120 F=((X-Z)*(X+Z)+H*(Y/T-H))/X QBD12400 C QBD12500 C NEXT QR SWEEP.. QBD12600 CS=ONE QBD12700 SN=ONE QBD12800 LP1=L+1 QBD12900 DO 160 I=LP1,K QBD13000 G=E(I) QBD13100 Y=Q(I) QBD13200 H=SN*G QBD13300 G=CS*G QBD13400 CALL G1 (F,H,CS,SN,E(I-1)) QBD13500 F=X*CS+G*SN QBD13600 G=-X*SN+G*CS QBD13700 H=Y*SN QBD13800 Y=Y*CS QBD13900 IF (.NOT.WNTV) GO TO 140 QBD14000 C QBD14100 C ACCUMULATE ROTATIONS (FROM THE RIGHT) IN 'V' QBD14200 DO 130 J=1,NRV QBD14300 130 CALL G2 (CS,SN,V(J,I-1),V(J,I)) QBD14400 140 CALL G1 (F,H,CS,SN,Q(I-1)) QBD14500 F=CS*G+SN*Y QBD14600 X=-SN*G+CS*Y QBD14700 IF (.NOT.HAVERS) GO TO 160 QBD14800 DO 150 J=1,NCC QBD14900 150 CALL G2 (CS,SN,C(I-1,J),C(I,J)) QBD15000 C APPLY ROTATIONS FROM THE LEFT TO QBD15100 C RIGHT HAND SIDES IN 'C'.. QBD15200 C QBD15300 160 CONTINUE QBD15400 E(L)=ZERO QBD15500 E(K)=F QBD15600 Q(K)=X QBD15700 NQRS=NQRS+1 QBD15800 IF (NQRS.LE.N10) GO TO 20 QBD15900 C RETURN TO 'TEST FOR SPLITTING'. QBD16000 C QBD16100 FAIL=.TRUE. QBD16200 C .. QBD16300 C CUTOFF FOR CONVERGENCE FAILURE. 'NQRS' WILL BE 2*N USUALLY. QBD16400 170 IF (Z.GE.ZERO) GO TO 190 QBD16500 Q(K)=-Z QBD16600 IF (.NOT.WNTV) GO TO 190 QBD16700 DO 180 J=1,NRV QBD16800 180 V(J,K)=-V(J,K) QBD16900 190 CONTINUE QBD17000 C CONVERGENCE. Q(K) IS MADE NONNEGATIVE.. QBD17100 C QBD17200 200 CONTINUE QBD17300 IF (N.EQ.1) RETURN QBD17400 DO 210 I=2,N QBD17500 IF (Q(I).GT.Q(I-1)) GO TO 220 QBD17600 210 CONTINUE QBD17700 IF (FAIL) IPASS=2 QBD17800 RETURN QBD17900 C .. QBD18000 C EVERY SINGULAR VALUE IS IN ORDER.. QBD18100 220 DO 270 I=2,N QBD18200 T=Q(I-1) QBD18300 K=I-1 QBD18400 DO 230 J=I,N QBD18500 IF (T.GE.Q(J)) GO TO 230 QBD18600 T=Q(J) QBD18700 K=J QBD18800 230 CONTINUE QBD18900 IF (K.EQ.I-1) GO TO 270 QBD19000 Q(K)=Q(I-1) QBD19100 Q(I-1)=T QBD19200 IF (.NOT.HAVERS) GO TO 250 QBD19300 DO 240 J=1,NCC QBD19400 T=C(I-1,J) QBD19500 C(I-1,J)=C(K,J) QBD19600 240 C(K,J)=T QBD19700 250 IF (.NOT.WNTV) GO TO 270 QBD19800 DO 260 J=1,NRV QBD19900 T=V(J,I-1) QBD20000 V(J,I-1)=V(J,K) QBD20100 260 V(J,K)=T QBD20200 270 CONTINUE QBD20300 C END OF ORDERING ALGORITHM. QBD20400 C QBD20500 IF (FAIL) IPASS=2 QBD20600 RETURN QBD20700 END QBD20800