home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 September / Simtel20_Sept92.cdr / msdos / fortran / linpklib.arc / SSVDC.FOR < prev    next >
Text File  |  1984-01-06  |  15KB  |  481 lines
  1.       SUBROUTINE SSVDC(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,JOB,INFO)
  2.       INTEGER LDX,N,P,LDU,LDV,JOB,INFO
  3.       REAL X(LDX,1),S(1),E(1),U(LDU,1),V(LDV,1),WORK(1)
  4. C
  5. C
  6. C     SSVDC IS A SUBROUTINE TO REDUCE A REAL NXP MATRIX X BY
  7. C     ORTHOGONAL TRANSFORMATIONS U AND V TO DIAGONAL FORM.  THE
  8. C     DIAGONAL ELEMENTS S(I) ARE THE SINGULAR VALUES OF X.  THE
  9. C     COLUMNS OF U ARE THE CORRESPONDING LEFT SINGULAR VECTORS,
  10. C     AND THE COLUMNS OF V THE RIGHT SINGULAR VECTORS.
  11. C
  12. C     ON ENTRY
  13. C
  14. C         X         REAL(LDX,P), WHERE LDX.GE.N.
  15. C                   X CONTAINS THE MATRIX WHOSE SINGULAR VALUE
  16. C                   DECOMPOSITION IS TO BE COMPUTED.  X IS
  17. C                   DESTROYED BY SSVDC.
  18. C
  19. C         LDX       INTEGER.
  20. C                   LDX IS THE LEADING DIMENSION OF THE ARRAY X.
  21. C
  22. C         N         INTEGER.
  23. C                   N IS THE NUMBER OF COLUMNS OF THE MATRIX X.
  24. C
  25. C         P         INTEGER.
  26. C                   P IS THE NUMBER OF ROWS OF THE MATRIX X.
  27. C
  28. C         LDU       INTEGER.
  29. C                   LDU IS THE LEADING DIMENSION OF THE ARRAY U.
  30. C                   (SEE BELOW).
  31. C
  32. C         LDV       INTEGER.
  33. C                   LDV IS THE LEADING DIMENSION OF THE ARRAY V.
  34. C                   (SEE BELOW).
  35. C
  36. C         WORK      REAL(N).
  37. C                   WORK IS A SCRATCH ARRAY.
  38. C
  39. C         JOB       INTEGER.
  40. C                   JOB CONTROLS THE COMPUTATION OF THE SINGULAR
  41. C                   VECTORS.  IT HAS THE DECIMAL EXPANSION AB
  42. C                   WITH THE FOLLOWING MEANING
  43. C
  44. C                        A.EQ.0    DO NOT COMPUTE THE LEFT SINGULAR
  45. C                                  VECTORS.
  46. C                        A.EQ.1    RETURN THE N LEFT SINGULAR VECTORS
  47. C                                  IN U.
  48. C                        A.GE.2    RETURN THE FIRST MIN(N,P) SINGULAR
  49. C                                  VECTORS IN U.
  50. C                        B.EQ.0    DO NOT COMPUTE THE RIGHT SINGULAR
  51. C                                  VECTORS.
  52. C                        B.EQ.1    RETURN THE RIGHT SINGULAR VECTORS
  53. C                                  IN V.
  54. C
  55. C     ON RETURN
  56. C
  57. C         S         REAL(MM), WHERE MM=MIN(N+1,P).
  58. C                   THE FIRST MIN(N,P) ENTRIES OF S CONTAIN THE
  59. C                   SINGULAR VALUES OF X ARRANGED IN DESCENDING
  60. C                   ORDER OF MAGNITUDE.
  61. C
  62. C         E         REAL(P).
  63. C                   E ORDINARILY CONTAINS ZEROS.  HOWEVER SEE THE
  64. C                   DISCUSSION OF INFO FOR EXCEPTIONS.
  65. C
  66. C         U         REAL(LDU,K), WHERE LDU.GE.N.  IF JOBA.EQ.1 THEN
  67. C                                   K.EQ.N, IF JOBA.GE.2 THEN
  68. C                                   K.EQ.MIN(N,P).
  69. C                   U CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
  70. C                   U IS NOT REFERENCED IF JOBA.EQ.0.  IF N.LE.P
  71. C                   OR IF JOBA.EQ.2, THEN U MAY BE IDENTIFIED WITH X
  72. C                   IN THE SUBROUTINE CALL.
  73. C
  74. C         V         REAL(LDV,P), WHERE LDV.GE.P.
  75. C                   V CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
  76. C                   V IS NOT REFERENCED IF JOB.EQ.0.  IF P.LE.N,
  77. C                   THEN V MAY BE IDENTIFIED WITH X IN THE
  78. C                   SUBROUTINE CALL.
  79. C
  80. C         INFO      INTEGER.
  81. C                   THE SINGULAR VALUES (AND THEIR CORRESPONDING
  82. C                   SINGULAR VECTORS) S(INFO+1),S(INFO+2),...,S(M)
  83. C                   ARE CORRECT (HERE M=MIN(N,P)).  THUS IF
  84. C                   INFO.EQ.0, ALL THE SINGULAR VALUES AND THEIR
  85. C                   VECTORS ARE CORRECT.  IN ANY EVENT, THE MATRIX
  86. C                   B = TRANS(U)*X*V IS THE BIDIAGONAL MATRIX
  87. C                   WITH THE ELEMENTS OF S ON ITS DIAGONAL AND THE
  88. C                   ELEMENTS OF E ON ITS SUPER-DIAGONAL (TRANS(U)
  89. C                   IS THE TRANSPOSE OF U).  THUS THE SINGULAR
  90. C                   VALUES OF X AND B ARE THE SAME.
  91. C
  92. C     LINPACK. THIS VERSION DATED 03/19/79 .
  93. C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
  94. C
  95. C     ***** USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.
  96. C
  97. C     EXTERNAL SROT
  98. C     BLAS SAXPY,SDOT,SSCAL,SSWAP,SNRM2,SROTG
  99. C     FORTRAN ABS,AMAX1,MAX0,MIN0,MOD,SQRT
  100. C
  101. C     INTERNAL VARIABLES
  102. C
  103.       INTEGER I,ITER,J,JOBU,K,KASE,KK,L,LL,LLS,LM1,LP1,LS,LU,M,MAXIT,
  104.      *        MM,MM1,MP1,NCT,NCTP1,NCU,NRT,NRTP1
  105.       REAL SDOT,T,R
  106.       REAL B,C,CS,EL,EMM1,F,G,SNRM2,SCALE,SHIFT,SL,SM,SN,SMM1,T1,TEST,
  107.      *     ZTEST
  108.       LOGICAL WANTU,WANTV
  109. C
  110. C
  111. C     SET THE MAXIMUM NUMBER OF ITERATIONS.
  112. C
  113.       MAXIT = 30
  114. C
  115. C     DETERMINE WHAT IS TO BE COMPUTED.
  116. C
  117.       WANTU = .FALSE.
  118.       WANTV = .FALSE.
  119.       JOBU = MOD(JOB,100)/10
  120.       NCU = N
  121.       IF (JOBU .GT. 1) NCU = MIN0(N,P)
  122.       IF (JOBU .NE. 0) WANTU = .TRUE.
  123.       IF (MOD(JOB,10) .NE. 0) WANTV = .TRUE.
  124. C
  125. C     REDUCE X TO BIDIAGONAL FORM, STORING THE DIAGONAL ELEMENTS
  126. C     IN S AND THE SUPER-DIAGONAL ELEMENTS IN E.
  127. C
  128.       INFO = 0
  129.       NCT = MIN0(N-1,P)
  130.       NRT = MAX0(0,MIN0(P-2,N))
  131.       LU = MAX0(NCT,NRT)
  132.       IF (LU .LT. 1) GO TO 170
  133.       DO 160 L = 1, LU
  134.          LP1 = L + 1
  135.          IF (L .GT. NCT) GO TO 20
  136. C
  137. C           COMPUTE THE TRANSFORMATION FOR THE L-TH COLUMN AND
  138. C           PLACE THE L-TH DIAGONAL IN S(L).
  139. C
  140.             S(L) = SNRM2(N-L+1,X(L,L),1)
  141.             IF (S(L) .EQ. 0.0E0) GO TO 10
  142.                IF (X(L,L) .NE. 0.0E0) S(L) = SIGN(S(L),X(L,L))
  143.                CALL SSCAL(N-L+1,1.0E0/S(L),X(L,L),1)
  144.                X(L,L) = 1.0E0 + X(L,L)
  145.    10       CONTINUE
  146.             S(L) = -S(L)
  147.    20    CONTINUE
  148.          IF (P .LT. LP1) GO TO 50
  149.          DO 40 J = LP1, P
  150.             IF (L .GT. NCT) GO TO 30
  151.             IF (S(L) .EQ. 0.0E0) GO TO 30
  152. C
  153. C              APPLY THE TRANSFORMATION.
  154. C
  155.                T = -SDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L)
  156.                CALL SAXPY(N-L+1,T,X(L,L),1,X(L,J),1)
  157.    30       CONTINUE
  158. C
  159. C           PLACE THE L-TH ROW OF X INTO  E FOR THE
  160. C           SUBSEQUENT CALCULATION OF THE ROW TRANSFORMATION.
  161. C
  162.             E(J) = X(L,J)
  163.    40    CONTINUE
  164.    50    CONTINUE
  165.          IF (.NOT.WANTU .OR. L .GT. NCT) GO TO 70
  166. C
  167. C           PLACE THE TRANSFORMATION IN U FOR SUBSEQUENT BACK
  168. C           MULTIPLICATION.
  169. C
  170.             DO 60 I = L, N
  171.                U(I,L) = X(I,L)
  172.    60       CONTINUE
  173.    70    CONTINUE
  174.          IF (L .GT. NRT) GO TO 150
  175. C
  176. C           COMPUTE THE L-TH ROW TRANSFORMATION AND PLACE THE
  177. C           L-TH SUPER-DIAGONAL IN E(L).
  178. C
  179.             E(L) = SNRM2(P-L,E(LP1),1)
  180.             IF (E(L) .EQ. 0.0E0) GO TO 80
  181.                IF (E(LP1) .NE. 0.0E0) E(L) = SIGN(E(L),E(LP1))
  182.                CALL SSCAL(P-L,1.0E0/E(L),E(LP1),1)
  183.                E(LP1) = 1.0E0 + E(LP1)
  184.    80       CONTINUE
  185.             E(L) = -E(L)
  186.             IF (LP1 .GT. N .OR. E(L) .EQ. 0.0E0) GO TO 120
  187. C
  188. C              APPLY THE TRANSFORMATION.
  189. C
  190.                DO 90 I = LP1, N
  191.                   WORK(I) = 0.0E0
  192.    90          CONTINUE
  193.                DO 100 J = LP1, P
  194.                   CALL SAXPY(N-L,E(J),X(LP1,J),1,WORK(LP1),1)
  195.   100          CONTINUE
  196.                DO 110 J = LP1, P
  197.                   CALL SAXPY(N-L,-E(J)/E(LP1),WORK(LP1),1,X(LP1,J),1)
  198.   110          CONTINUE
  199.   120       CONTINUE
  200.             IF (.NOT.WANTV) GO TO 140
  201. C
  202. C              PLACE THE TRANSFORMATION IN V FOR SUBSEQUENT
  203. C              BACK MULTIPLICATION.
  204. C
  205.                DO 130 I = LP1, P
  206.                   V(I,L) = E(I)
  207.   130          CONTINUE
  208.   140       CONTINUE
  209.   150    CONTINUE
  210.   160 CONTINUE
  211.   170 CONTINUE
  212. C
  213. C     SET UP THE FINAL BIDIAGONAL MATRIX OR ORDER M.
  214. C
  215.       M = MIN0(P,N+1)
  216.       NCTP1 = NCT + 1
  217.       NRTP1 = NRT + 1
  218.       IF (NCT .LT. P) S(NCTP1) = X(NCTP1,NCTP1)
  219.       IF (N .LT. M) S(M) = 0.0E0
  220.       IF (NRTP1 .LT. M) E(NRTP1) = X(NRTP1,M)
  221.       E(M) = 0.0E0
  222. C
  223. C     IF REQUIRED, GENERATE U.
  224. C
  225.       IF (.NOT.WANTU) GO TO 300
  226.          IF (NCU .LT. NCTP1) GO TO 200
  227.          DO 190 J = NCTP1, NCU
  228.             DO 180 I = 1, N
  229.                U(I,J) = 0.0E0
  230.   180       CONTINUE
  231.             U(J,J) = 1.0E0
  232.   190    CONTINUE
  233.   200    CONTINUE
  234.          IF (NCT .LT. 1) GO TO 290
  235.          DO 280 LL = 1, NCT
  236.             L = NCT - LL + 1
  237.             IF (S(L) .EQ. 0.0E0) GO TO 250
  238.                LP1 = L + 1
  239.                IF (NCU .LT. LP1) GO TO 220
  240.                DO 210 J = LP1, NCU
  241.                   T = -SDOT(N-L+1,U(L,L),1,U(L,J),1)/U(L,L)
  242.                   CALL SAXPY(N-L+1,T,U(L,L),1,U(L,J),1)
  243.   210          CONTINUE
  244.   220          CONTINUE
  245.                CALL SSCAL(N-L+1,-1.0E0,U(L,L),1)
  246.                U(L,L) = 1.0E0 + U(L,L)
  247.                LM1 = L - 1
  248.                IF (LM1 .LT. 1) GO TO 240
  249.                DO 230 I = 1, LM1
  250.                   U(I,L) = 0.0E0
  251.   230          CONTINUE
  252.   240          CONTINUE
  253.             GO TO 270
  254.   250       CONTINUE
  255.                DO 260 I = 1, N
  256.                   U(I,L) = 0.0E0
  257.   260          CONTINUE
  258.                U(L,L) = 1.0E0
  259.   270       CONTINUE
  260.   280    CONTINUE
  261.   290    CONTINUE
  262.   300 CONTINUE
  263. C
  264. C     IF IT IS REQUIRED, GENERATE V.
  265. C
  266.       IF (.NOT.WANTV) GO TO 350
  267.          DO 340 LL = 1, P
  268.             L = P - LL + 1
  269.             LP1 = L + 1
  270.             IF (L .GT. NRT) GO TO 320
  271.             IF (E(L) .EQ. 0.0E0) GO TO 320
  272.                DO 310 J = LP1, P
  273.                   T = -SDOT(P-L,V(LP1,L),1,V(LP1,J),1)/V(LP1,L)
  274.                   CALL SAXPY(P-L,T,V(LP1,L),1,V(LP1,J),1)
  275.   310          CONTINUE
  276.   320       CONTINUE
  277.             DO 330 I = 1, P
  278.                V(I,L) = 0.0E0
  279.   330       CONTINUE
  280.             V(L,L) = 1.0E0
  281.   340    CONTINUE
  282.   350 CONTINUE
  283. C
  284. C     MAIN ITERATION LOOP FOR THE SINGULAR VALUES.
  285. C
  286.       MM = M
  287.       ITER = 0
  288.   360 CONTINUE
  289. C
  290. C        QUIT IF ALL THE SINGULAR VALUES HAVE BEEN FOUND.
  291. C
  292. C     ...EXIT
  293.          IF (M .EQ. 0) GO TO 620
  294. C
  295. C        IF TOO MANY ITERATIONS HAVE BEEN PERFORMED, SET
  296. C        FLAG AND RETURN.
  297. C
  298.          IF (ITER .LT. MAXIT) GO TO 370
  299.             INFO = M
  300. C     ......EXIT
  301.             GO TO 620
  302.   370    CONTINUE
  303. C
  304. C        THIS SECTION OF THE PROGRAM INSPECTS FOR
  305. C        NEGLIGIBLE ELEMENTS IN THE S AND E ARRAYS.  ON
  306. C        COMPLETION THE VARIABLES KASE AND L ARE SET AS FOLLOWS.
  307. C
  308. C           KASE = 1     IF S(M) AND E(L-1) ARE NEGLIGIBLE AND L.LT.M
  309. C           KASE = 2     IF S(L) IS NEGLIGIBLE AND L.LT.M
  310. C           KASE = 3     IF E(L-1) IS NEGLIGIBLE, L.LT.M, AND
  311. C                        S(L), ..., S(M) ARE NOT NEGLIGIBLE (QR STEP).
  312. C           KASE = 4     IF E(M-1) IS NEGLIGIBLE (CONVERGENCE).
  313. C
  314.          DO 390 LL = 1, M
  315.             L = M - LL
  316. C        ...EXIT
  317.             IF (L .EQ. 0) GO TO 400
  318.             TEST = ABS(S(L)) + ABS(S(L+1))
  319.             ZTEST = TEST + ABS(E(L))
  320.             IF (ZTEST .NE. TEST) GO TO 380
  321.                E(L) = 0.0E0
  322. C        ......EXIT
  323.                GO TO 400
  324.   380       CONTINUE
  325.   390    CONTINUE
  326.   400    CONTINUE
  327.          IF (L .NE. M - 1) GO TO 410
  328.             KASE = 4
  329.          GO TO 480
  330.   410    CONTINUE
  331.             LP1 = L + 1
  332.             MP1 = M + 1
  333.             DO 430 LLS = LP1, MP1
  334.                LS = M - LLS + LP1
  335. C           ...EXIT
  336.                IF (LS .EQ. L) GO TO 440
  337.                TEST = 0.0E0
  338.                IF (LS .NE. M) TEST = TEST + ABS(E(LS))
  339.                IF (LS .NE. L + 1) TEST = TEST + ABS(E(LS-1))
  340.                ZTEST = TEST + ABS(S(LS))
  341.                IF (ZTEST .NE. TEST) GO TO 420
  342.                   S(LS) = 0.0E0
  343. C           ......EXIT
  344.                   GO TO 440
  345.   420          CONTINUE
  346.   430       CONTINUE
  347.   440       CONTINUE
  348.             IF (LS .NE. L) GO TO 450
  349.                KASE = 3
  350.             GO TO 470
  351.   450       CONTINUE
  352.             IF (LS .NE. M) GO TO 460
  353.                KASE = 1
  354.             GO TO 470
  355.   460       CONTINUE
  356.                KASE = 2
  357.                L = LS
  358.   470       CONTINUE
  359.   480    CONTINUE
  360.          L = L + 1
  361. C
  362. C        PERFORM THE TASK INDICATED BY KASE.
  363. C
  364.          GO TO (490,520,540,570), KASE
  365. C
  366. C        DEFLATE NEGLIGIBLE S(M).
  367. C
  368.   490    CONTINUE
  369.             MM1 = M - 1
  370.             F = E(M-1)
  371.             E(M-1) = 0.0E0
  372.             DO 510 KK = L, MM1
  373.                K = MM1 - KK + L
  374.                T1 = S(K)
  375.                CALL SROTG(T1,F,CS,SN)
  376.                S(K) = T1
  377.                IF (K .EQ. L) GO TO 500
  378.                   F = -SN*E(K-1)
  379.                   E(K-1) = CS*E(K-1)
  380.   500          CONTINUE
  381.                IF (WANTV) CALL SROT(P,V(1,K),1,V(1,M),1,CS,SN)
  382.   510       CONTINUE
  383.          GO TO 610
  384. C
  385. C        SPLIT AT NEGLIGIBLE S(L).
  386. C
  387.   520    CONTINUE
  388.             F = E(L-1)
  389.             E(L-1) = 0.0E0
  390.             DO 530 K = L, M
  391.                T1 = S(K)
  392.                CALL SROTG(T1,F,CS,SN)
  393.                S(K) = T1
  394.                F = -SN*E(K)
  395.                E(K) = CS*E(K)
  396.                IF (WANTU) CALL SROT(N,U(1,K),1,U(1,L-1),1,CS,SN)
  397.   530       CONTINUE
  398.          GO TO 610
  399. C
  400. C        PERFORM ONE QR STEP.
  401. C
  402.   540    CONTINUE
  403. C
  404. C           CALCULATE THE SHIFT.
  405. C
  406.             SCALE = AMAX1(ABS(S(M)),ABS(S(M-1)),ABS(E(M-1)),ABS(S(L)),
  407.      *                    ABS(E(L)))
  408.             SM = S(M)/SCALE
  409.             SMM1 = S(M-1)/SCALE
  410.             EMM1 = E(M-1)/SCALE
  411.             SL = S(L)/SCALE
  412.             EL = E(L)/SCALE
  413.             B = ((SMM1 + SM)*(SMM1 - SM) + EMM1**2)/2.0E0
  414.             C = (SM*EMM1)**2
  415.             SHIFT = 0.0E0
  416.             IF (B .EQ. 0.0E0 .AND. C .EQ. 0.0E0) GO TO 550
  417.                SHIFT = SQRT(B**2+C)
  418.                IF (B .LT. 0.0E0) SHIFT = -SHIFT
  419.                SHIFT = C/(B + SHIFT)
  420.   550       CONTINUE
  421.             F = (SL + SM)*(SL - SM) + SHIFT
  422.             G = SL*EL
  423. C
  424. C           CHASE ZEROS.
  425. C
  426.             MM1 = M - 1
  427.             DO 560 K = L, MM1
  428.                CALL SROTG(F,G,CS,SN)
  429.                IF (K .NE. L) E(K-1) = F
  430.                F = CS*S(K) + SN*E(K)
  431.                E(K) = CS*E(K) - SN*S(K)
  432.                G = SN*S(K+1)
  433.                S(K+1) = CS*S(K+1)
  434.                IF (WANTV) CALL SROT(P,V(1,K),1,V(1,K+1),1,CS,SN)
  435.                CALL SROTG(F,G,CS,SN)
  436.                S(K) = F
  437.                F = CS*E(K) + SN*S(K+1)
  438.                S(K+1) = -SN*E(K) + CS*S(K+1)
  439.                G = SN*E(K+1)
  440.                E(K+1) = CS*E(K+1)
  441.                IF (WANTU .AND. K .LT. N)
  442.      *            CALL SROT(N,U(1,K),1,U(1,K+1),1,CS,SN)
  443.   560       CONTINUE
  444.             E(M-1) = F
  445.             ITER = ITER + 1
  446.          GO TO 610
  447. C
  448. C        CONVERGENCE.
  449. C
  450.   570    CONTINUE
  451. C
  452. C           MAKE THE SINGULAR VALUE  POSITIVE.
  453. C
  454.             IF (S(L) .GE. 0.0E0) GO TO 580
  455.                S(L) = -S(L)
  456.                IF (WANTV) CALL SSCAL(P,-1.0E0,V(1,L),1)
  457.   580       CONTINUE
  458. C
  459. C           ORDER THE SINGULAR VALUE.
  460. C
  461.   590       IF (L .EQ. MM) GO TO 600
  462. C           ...EXIT
  463.                IF (S(L) .GE. S(L+1)) GO TO 600
  464.                T = S(L)
  465.                S(L) = S(L+1)
  466.                S(L+1) = T
  467.                IF (WANTV .AND. L .LT. P)
  468.      *            CALL SSWAP(P,V(1,L),1,V(1,L+1),1)
  469.                IF (WANTU .AND. L .LT. N)
  470.      *            CALL SSWAP(N,U(1,L),1,U(1,L+1),1)
  471.                L = L + 1
  472.             GO TO 590
  473.   600       CONTINUE
  474.             ITER = 0
  475.             M = M - 1
  476.   610    CONTINUE
  477.       GO TO 360
  478.   620 CONTINUE
  479.       RETURN
  480.       END
  481.