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