SUBROUTINE SVDRS (A,MDA,MM,NN,B,MDB,NB,S) C C INPUT.. C C A(*,*),MDA,M,N THE ARRAY A(*,*) IS DOUBLY SUBSCRIPTED WITH C FIRST DIMENSIONING PARAMETER EQUAL TO MDA. C THE ARRAY A(*,*) INITIALLY CONTAINS THE M BY C N MATRIX A. EITHER M.GE.N OR M.LT.N IS C PERMITTED. THE FIRST DIMENSIONING PARAMETER C MDA MUST SATISFY MDA.GE.MAX(M,N). THE C CONDITION MDA.LT.MAX(M,N) IS CONSIDERED AN C ERROR. C C B(*),MDB,NB NB DENOTES THE NUMBER OF COLUMN VECTORS IN THE C MATRIX B. IF NB = 0 THE ARRAY B(*) WILL NOT C BE REFERENCED BY THIS SUBROUTINE NOR WILL ITS C CONTENTS BE MODIFIED. IF NB.GE.2, THE ARRAY B C SHOULD BE DOUBLE SUBSCRIPTED WITH FIRST C DIMENSIONING PARAMETER EQUAL TO MDB.GE.M. IF C NB = 1, THEN B WILL BE USED AS A SINGLY C SUBSCRIPTED ARRAY OF LENGTH M. IN THIS LATTER C CASE THE VALUE OF MDB IS ARBITRARY BUT IT C SHOULD BE SET TO A DEFINED VALUE SUCH AS C MDB = M. THE CONTENTS OF THE ARRAY B(*) MUST C CONTAIN THE MATRIX B. THE CONDITION NB.GT.1 C AND. MDB.LT.M IS CONSIDERED AN ERROR. C C S(*) THE ARRAY S(*) IS USED AS 3*N LOCATIONS OF C TEMPORARY WORKING SPACE BY THE SUBROUTINE. C C OUTPUT.. C C A(*,*) ON OUTPUT THE ARRAY A(*,*) CONTAINS THE N BY N C ORTHOGONAL MATRIX V. C C B(*) ON OUTPUT THE ARRAY B(*),(OR B(*,*)), CONTAINS C THE PRODUCT OF THE TRANSPOSE OF THE M BY M C MATRIX U AND THE RIGHTHAND SIDE B. C C S(*) ON OUTPUT THE FIRST N LOCATIONS OF S(*) C CONTAIN THE DIAGONAL TERMS OF THE M BY N C DIAGONAL MATRIX S. THESE NONNEGATIVE TERMS C ARE ORDERED FROM LARGEST TO SMALLEST. IF C M.LT.N, ENTRIES M+1,...,N OF S(*) ARE ZERO. C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) C DIMENSION A(MDA,NN),B(MDB,NB),S(NN,3) C ZERO=0. ONE=1. C C BEGIN.. SPECIAL FOR ZERO ROWS AND COLS. C C PACK THE NONZERO COLS TO THE LEFT C N=NN IF (N.LE.0.OR.MM.LE.0) RETURN IF(MAX(MM,NN).GT.MDA) RETURN C IF(NB.GT.1.AND.MM.GT.MDB) RETURN C J=N 10 CONTINUE DO 20 I=1,MM IF (A(I,J)) 50,20,50 20 CONTINUE C C COL J IS ZERO. EXCHANGE IT WITH COL N. C IF (J.EQ.N) GO TO 40 DO 30 I=1,MM 30 A(I,J)=A(I,N) 40 CONTINUE A(1,N)=J N=N-1 50 CONTINUE J=J-1 IF (J.GE.1) GO TO 10 C IF N=0 THEN A IS ENTIRELY ZERO AND SVD C COMPUTATION CAN BE SKIPPED NS=0 IF (N.EQ.0) GO TO 240 C PACK NONZERO ROWS TO THE TOP C QUIT PACKING IF FIND N NONZERO ROWS I=1 M=MM 60 IF (I.GT.N.OR.I.GE.M) GO TO 150 IF (A(I,I)) 90,70,90 70 DO 80 J=1,N IF (A(I,J)) 90,80,90 80 CONTINUE GO TO 100 90 I=I+1 GO TO 60 C ROW I IS ZERO C EXCHANGE ROWS I AND M 100 IF(NB.LE.0) GO TO 115 DO 110 J=1,NB T=B(I,J) B(I,J)=B(M,J) 110 B(M,J)=T 115 DO 120 J=1,N 120 A(I,J)=A(M,J) IF (M.GT.N) GO TO 140 DO 130 J=1,N 130 A(M,J)=ZERO 140 CONTINUE C EXCHANGE IS FINISHED M=M-1 GO TO 60 C 150 CONTINUE C END.. SPECIAL FOR ZERO ROWS AND COLUMNS C BEGIN.. SVD ALGORITHM C METHOD.. C (1) REDUCE THE MATRIX TO UPPER BIDIAGONAL FORM WITH C HOUSEHOLDER TRANSFORMATIONS. C H(N)...H(1)AQ(1)...Q(N-2) = (D**T,0)**T C WHERE D IS UPPER BIDIAGONAL. C C (2) APPLY H(N)...H(1) TO B. HERE H(N)...H(1)*B REPLACES B C IN STORAGE. C C (3) THE MATRIX PRODUCT W= Q(1)...Q(N-2) OVERWRITES THE FIRST C N ROWS OF A IN STORAGE. C C (4) AN SVD FOR D IS COMPUTED. HERE K ROTATIONS RI AND PI ARE C COMPUTED SO THAT C RK...R1*D*P1**(T)...PK**(T) = DIAG(S1,...,SM) C TO WORKING ACCURACY. THE SI ARE NONNEGATIVE AND NONINCREASING. C HERE RK...R1*B OVERWRITES B IN STORAGE WHILE C A*P1**(T)...PK**(T) OVERWRITES A IN STORAGE. C C (5) IT FOLLOWS THAT,WITH THE PROPER DEFINITIONS, C U**(T)*B OVERWRITES B, WHILE V OVERWRITES THE FIRST N ROW AND C COLUMNS OF A. C L=MIN(M,N) C THE FOLLOWING LOOP REDUCES A TO UPPER BIDIAGONAL AND C ALSO APPLIES THE PREMULTIPLYING TRANSFORMATIONS TO B. C DO 170 J=1,L IF (J.GE.M) GO TO 160 CALL H12 (1,J,J+1,M,A(1,J),1,T,A(1,J+1),1,MDA,N-J) CALL H12 (2,J,J+1,M,A(1,J),1,T,B,1,MDB,NB) 160 IF (J.GE.N-1) GO TO 170 CALL H12 (1,J+1,J+2,N,A(J,1),MDA,S(J,3),A(J+1,1),MDA,1,M-J) 170 CONTINUE C C COPY THE BIDIAGONAL MATRIX INTO THE ARRAY S() FOR QRBD. C IF (N.EQ.1) GO TO 190 DO 180 J=2,N S(J,1)=A(J,J) 180 S(J,2)=A(J-1,J) 190 S(1,1)=A(1,1) C NS=N IF (M.GE.N) GO TO 200 NS=M+1 S(NS,1)=ZERO S(NS,2)=A(M,M+1) 200 CONTINUE C C CONSTRUCT THE EXPLICIT N BY N PRODUCT MATRIX, W=Q1*Q2*...*QL*I C IN THE ARRAY A(). C DO 230 K=1,N I=N+1-K IF (I.GT. MIN(M,N-2)) GO TO 210 CALL H12 (2,I+1,I+2,N,A(I,1),MDA,S(I,3),A(1,I+1),1,MDA,N-I) 210 DO 220 J=1,N 220 A(I,J)=ZERO 230 A(I,I)=ONE C C COMPUTE THE SVD OF THE BIDIAGONAL MATRIX C CALL QRBD (IPASS,S(1,1),S(1,2),NS,A,MDA,N,B,MDB,NB) C GO TO (240,310), IPASS 240 CONTINUE IF (NS.GE.N) GO TO 260 NSP1=NS+1 DO 250 J=NSP1,N 250 S(J,1)=ZERO 260 CONTINUE IF (N.EQ.NN) RETURN NP1=N+1 C MOVE RECORD OF PERMUTATIONS C AND STORE ZEROS DO 280 J=NP1,NN S(J,1)=A(1,J) DO 270 I=1,N 270 A(I,J)=ZERO 280 CONTINUE C PERMUTE ROWS AND SET ZERO SINGULAR VALUES. DO 300 K=NP1,NN I=S(K,1) S(K,1)=ZERO DO 290 J=1,NN A(K,J)=A(I,J) 290 A(I,J)=ZERO A(I,K)=ONE 300 CONTINUE C END.. SPECIAL FOR ZERO ROWS AND COLUMNS RETURN 310 CONTINUE C RETURN END