home *** CD-ROM | disk | FTP | other *** search
/ Crawly Crypt Collection 1 / crawlyvol1.bin / apps / math / matlab / src / lib3.f < prev    next >
Text File  |  1992-04-21  |  28KB  |  800 lines
  1.       SUBROUTINE WSVDC(XR,XI,LDX,N,P,SR,SI,ER,EI,UR,UI,LDU,VR,VI,LDV,
  2.      *                 WORKR,WORKI,JOB,INFO)
  3.       INTEGER LDX,N,P,LDU,LDV,JOB,INFO
  4.       DOUBLE PRECISION XR(LDX,1),XI(LDX,1),SR(1),SI(1),ER(1),EI(1),
  5.      *                 UR(LDU,1),UI(LDU,1),VR(LDV,1),VI(LDV,1),
  6.      *                 WORKR(1),WORKI(1)
  7. C
  8. C
  9. C     WSVDC IS A SUBROUTINE TO REDUCE A DOUBLE-COMPLEX NXP MATRIX X BY
  10. C     UNITARY TRANSFORMATIONS U AND V TO DIAGONAL FORM.  THE
  11. C     DIAGONAL ELEMENTS S(I) ARE THE SINGULAR VALUES OF X.  THE
  12. C     COLUMNS OF U ARE THE CORRESPONDING LEFT SINGULAR VECTORS,
  13. C     AND THE COLUMNS OF V THE RIGHT SINGULAR VECTORS.
  14. C
  15. C     ON ENTRY
  16. C
  17. C         X         DOUBLE-COMPLEX(LDX,P), WHERE LDX.GE.N.
  18. C                   X CONTAINS THE MATRIX WHOSE SINGULAR VALUE
  19. C                   DECOMPOSITION IS TO BE COMPUTED.  X IS
  20. C                   DESTROYED BY WSVDC.
  21. C
  22. C         LDX       INTEGER.
  23. C                   LDX IS THE LEADING DIMENSION OF THE ARRAY X.
  24. C
  25. C         N         INTEGER.
  26. C                   N IS THE NUMBER OF COLUMNS OF THE MATRIX X.
  27. C
  28. C         P         INTEGER.
  29. C                   P IS THE NUMBER OF ROWS OF THE MATRIX X.
  30. C
  31. C         LDU       INTEGER.
  32. C                   LDU IS THE LEADING DIMENSION OF THE ARRAY U
  33. C                   (SEE BELOW).
  34. C
  35. C         LDV       INTEGER.
  36. C                   LDV IS THE LEADING DIMENSION OF THE ARRAY V
  37. C                   (SEE BELOW).
  38. C
  39. C         WORK      DOUBLE-COMPLEX(N).
  40. C                   WORK IS A SCRATCH ARRAY.
  41. C
  42. C         JOB       INTEGER.
  43. C                   JOB CONTROLS THE COMPUTATION OF THE SINGULAR
  44. C                   VECTORS.  IT HAS THE DECIMAL EXPANSION AB
  45. C                   WITH THE FOLLOWING MEANING
  46. C
  47. C                        A.EQ.0    DO NOT COMPUTE THE LEFT SINGULAR
  48. C                                  VECTORS.
  49. C                        A.EQ.1    RETURN THE N LEFT SINGULAR VECTORS
  50. C                                  IN U.
  51. C                        A.GE.2    RETURNS THE FIRST MIN(N,P)
  52. C                                  LEFT SINGULAR VECTORS IN U.
  53. C                        B.EQ.0    DO NOT COMPUTE THE RIGHT SINGULAR
  54. C                                  VECTORS.
  55. C                        B.EQ.1    RETURN THE RIGHT SINGULAR VECTORS
  56. C                                  IN V.
  57. C
  58. C     ON RETURN
  59. C
  60. C         S         DOUBLE-COMPLEX(MM), WHERE MM=MIN(N+1,P).
  61. C                   THE FIRST MIN(N,P) ENTRIES OF S CONTAIN THE
  62. C                   SINGULAR VALUES OF X ARRANGED IN DESCENDING
  63. C                   ORDER OF MAGNITUDE.
  64. C
  65. C         E         DOUBLE-COMPLEX(P).
  66. C                   E ORDINARILY CONTAINS ZEROS.  HOWEVER SEE THE
  67. C                   DISCUSSION OF INFO FOR EXCEPTIONS.
  68. C
  69. C         U         DOUBLE-COMPLEX(LDU,K), WHERE LDU.GE.N.
  70. C                   IF JOBA.EQ.1 THEN K.EQ.N,
  71. C                   IF JOBA.EQ.2 THEN K.EQ.MIN(N,P).
  72. C                   U CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
  73. C                   U IS NOT REFERENCED IF JOBA.EQ.0.  IF N.LE.P
  74. C                   OR IF JOBA.GT.2, THEN U MAY BE IDENTIFIED WITH X
  75. C                   IN THE SUBROUTINE CALL.
  76. C
  77. C         V         DOUBLE-COMPLEX(LDV,P), WHERE LDV.GE.P.
  78. C                   V CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
  79. C                   V IS NOT REFERENCED IF JOBB.EQ.0.  IF P.LE.N,
  80. C                   THEN V MAY BE IDENTIFIED WHTH X IN THE
  81. C                   SUBROUTINE CALL.
  82. C
  83. C         INFO      INTEGER.
  84. C                   THE SINGULAR VALUES (AND THEIR CORRESPONDING
  85. C                   SINGULAR VECTORS) S(INFO+1),S(INFO+2),...,S(M)
  86. C                   ARE CORRECT (HERE M=MIN(N,P)).  THUS IF
  87. C                   INFO.EQ.0, ALL THE SINGULAR VALUES AND THEIR
  88. C                   VECTORS ARE CORRECT.  IN ANY EVENT, THE MATRIX
  89. C                   B = CTRANS(U)*X*V IS THE BIDIAGONAL MATRIX
  90. C                   WITH THE ELEMENTS OF S ON ITS DIAGONAL AND THE
  91. C                   ELEMENTS OF E ON ITS SUPER-DIAGONAL (CTRANS(U)
  92. C                   IS THE CONJUGATE-TRANSPOSE OF U).  THUS THE
  93. C                   SINGULAR VALUES OF X AND B ARE THE SAME.
  94. C
  95. C     LINPACK. THIS VERSION DATED 07/03/79 .
  96. C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
  97. C
  98. C     WSVDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.
  99. C
  100. C     BLAS WAXPY,PYTHAG,WDOTCR,WDOTCI,WSCAL,WSWAP,WNRM2,RROTG
  101. C     FORTRAN DABS,DIMAG,DMAX1
  102. C     FORTRAN MAX0,MIN0,MOD,DSQRT
  103. C
  104. C     INTERNAL VARIABLES
  105. C
  106.       INTEGER I,ITER,J,JOBU,K,KASE,KK,L,LL,LLS,LM1,LP1,LS,LU,M,MAXIT,
  107.      *        MM,MM1,MP1,NCT,NCTP1,NCU,NRT,NRTP1
  108.       DOUBLE PRECISION PYTHAG,WDOTCR,WDOTCI,TR,TI,RR,RI
  109.       DOUBLE PRECISION B,C,CS,EL,EMM1,F,G,WNRM2,SCALE,SHIFT,SL,SM,SN,
  110.      *                 SMM1,T1,TEST,ZTEST,SMALL,FLOP
  111.       LOGICAL WANTU,WANTV
  112. C
  113.       DOUBLE PRECISION ZDUMR,ZDUMI
  114.       DOUBLE PRECISION CABS1
  115. C
  116. C     SET THE MAXIMUM NUMBER OF ITERATIONS.
  117. C
  118.       MAXIT = 75
  119. C
  120. C     SMALL NUMBER, ROUGHLY MACHINE EPSILON, USED TO AVOID UNDERFLOW
  121. C
  122.       SMALL = 1.D0/2.D0**48
  123. C
  124. C     DETERMINE WHAT IS TO BE COMPUTED.
  125. C
  126.       WANTU = .FALSE.
  127.       WANTV = .FALSE.
  128.       JOBU = MOD(JOB,100)/10
  129.       NCU = N
  130.       IF (JOBU .GT. 1) NCU = MIN0(N,P)
  131.       IF (JOBU .NE. 0) WANTU = .TRUE.
  132.       IF (MOD(JOB,10) .NE. 0) WANTV = .TRUE.
  133. C
  134. C     REDUCE X TO BIDIAGONAL FORM, STORING THE DIAGONAL ELEMENTS
  135. C     IN S AND THE SUPER-DIAGONAL ELEMENTS IN E.
  136. C
  137.       INFO = 0
  138.       NCT = MIN0(N-1,P)
  139.       NRT = MAX0(0,MIN0(P-2,N))
  140.       LU = MAX0(NCT,NRT)
  141.       IF (LU .LT. 1) GO TO 190
  142.       DO 180 L = 1, LU
  143.          LP1 = L + 1
  144.          IF (L .GT. NCT) GO TO 30
  145. C
  146. C           COMPUTE THE TRANSFORMATION FOR THE L-TH COLUMN AND
  147. C           PLACE THE L-TH DIAGONAL IN S(L).
  148. C
  149.             SR(L) = WNRM2(N-L+1,XR(L,L),XI(L,L),1)
  150.             SI(L) = 0.0D0
  151.             IF (CABS1(SR(L),SI(L)) .EQ. 0.0D0) GO TO 20
  152.                IF (CABS1(XR(L,L),XI(L,L)) .EQ. 0.0D0) GO TO 10
  153.                   CALL WSIGN(SR(L),SI(L),XR(L,L),XI(L,L),SR(L),SI(L))
  154.    10          CONTINUE
  155.                CALL WDIV(1.0D0,0.0D0,SR(L),SI(L),TR,TI)
  156.                CALL WSCAL(N-L+1,TR,TI,XR(L,L),XI(L,L),1)
  157.                XR(L,L) = FLOP(1.0D0 + XR(L,L))
  158.    20       CONTINUE
  159.             SR(L) = -SR(L)
  160.             SI(L) = -SI(L)
  161.    30    CONTINUE
  162.          IF (P .LT. LP1) GO TO 60
  163.          DO 50 J = LP1, P
  164.             IF (L .GT. NCT) GO TO 40
  165.             IF (CABS1(SR(L),SI(L)) .EQ. 0.0D0) GO TO 40
  166. C
  167. C              APPLY THE TRANSFORMATION.
  168. C
  169.                TR = -WDOTCR(N-L+1,XR(L,L),XI(L,L),1,XR(L,J),XI(L,J),1)
  170.                TI = -WDOTCI(N-L+1,XR(L,L),XI(L,L),1,XR(L,J),XI(L,J),1)
  171.                CALL WDIV(TR,TI,XR(L,L),XI(L,L),TR,TI)
  172.                CALL WAXPY(N-L+1,TR,TI,XR(L,L),XI(L,L),1,XR(L,J),
  173.      *                    XI(L,J),1)
  174.    40       CONTINUE
  175. C
  176. C           PLACE THE L-TH ROW OF X INTO  E FOR THE
  177. C           SUBSEQUENT CALCULATION OF THE ROW TRANSFORMATION.
  178. C
  179.             ER(J) = XR(L,J)
  180.             EI(J) = -XI(L,J)
  181.    50    CONTINUE
  182.    60    CONTINUE
  183.          IF (.NOT.WANTU .OR. L .GT. NCT) GO TO 80
  184. C
  185. C           PLACE THE TRANSFORMATION IN U FOR SUBSEQUENT BACK
  186. C           MULTIPLICATION.
  187. C
  188.             DO 70 I = L, N
  189.                UR(I,L) = XR(I,L)
  190.                UI(I,L) = XI(I,L)
  191.    70       CONTINUE
  192.    80    CONTINUE
  193.          IF (L .GT. NRT) GO TO 170
  194. C
  195. C           COMPUTE THE L-TH ROW TRANSFORMATION AND PLACE THE
  196. C           L-TH SUPER-DIAGONAL IN E(L).
  197. C
  198.             ER(L) = WNRM2(P-L,ER(LP1),EI(LP1),1)
  199.             EI(L) = 0.0D0
  200.             IF (CABS1(ER(L),EI(L)) .EQ. 0.0D0) GO TO 100
  201.                IF (CABS1(ER(LP1),EI(LP1)) .EQ. 0.0D0) GO TO 90
  202.                   CALL WSIGN(ER(L),EI(L),ER(LP1),EI(LP1),ER(L),EI(L))
  203.    90          CONTINUE
  204.                CALL WDIV(1.0D0,0.0D0,ER(L),EI(L),TR,TI)
  205.                CALL WSCAL(P-L,TR,TI,ER(LP1),EI(LP1),1)
  206.                ER(LP1) = FLOP(1.0D0 + ER(LP1))
  207.   100       CONTINUE
  208.             ER(L) = -ER(L)
  209.             EI(L) = +EI(L)
  210.             IF (LP1 .GT. N .OR. CABS1(ER(L),EI(L)) .EQ. 0.0D0)
  211.      *         GO TO 140
  212. C
  213. C              APPLY THE TRANSFORMATION.
  214. C
  215.                DO 110 I = LP1, N
  216.                   WORKR(I) = 0.0D0
  217.                   WORKI(I) = 0.0D0
  218.   110          CONTINUE
  219.                DO 120 J = LP1, P
  220.                   CALL WAXPY(N-L,ER(J),EI(J),XR(LP1,J),XI(LP1,J),1,
  221.      *                       WORKR(LP1),WORKI(LP1),1)
  222.   120          CONTINUE
  223.                DO 130 J = LP1, P
  224.                   CALL WDIV(-ER(J),-EI(J),ER(LP1),EI(LP1),TR,TI)
  225.                   CALL WAXPY(N-L,TR,-TI,WORKR(LP1),WORKI(LP1),1,
  226.      *                       XR(LP1,J),XI(LP1,J),1)
  227.   130          CONTINUE
  228.   140       CONTINUE
  229.             IF (.NOT.WANTV) GO TO 160
  230. C
  231. C              PLACE THE TRANSFORMATION IN V FOR SUBSEQUENT
  232. C              BACK MULTIPLICATION.
  233. C
  234.                DO 150 I = LP1, P
  235.                   VR(I,L) = ER(I)
  236.                   VI(I,L) = EI(I)
  237.   150          CONTINUE
  238.   160       CONTINUE
  239.   170    CONTINUE
  240.   180 CONTINUE
  241.   190 CONTINUE
  242. C
  243. C     SET UP THE FINAL BIDIAGONAL MATRIX OR ORDER M.
  244. C
  245.       M = MIN0(P,N+1)
  246.       NCTP1 = NCT + 1
  247.       NRTP1 = NRT + 1
  248.       IF (NCT .GE. P) GO TO 200
  249.          SR(NCTP1) = XR(NCTP1,NCTP1)
  250.          SI(NCTP1) = XI(NCTP1,NCTP1)
  251.   200 CONTINUE
  252.       IF (N .GE. M) GO TO 210
  253.          SR(M) = 0.0D0
  254.          SI(M) = 0.0D0
  255.   210 CONTINUE
  256.       IF (NRTP1 .GE. M) GO TO 220
  257.          ER(NRTP1) = XR(NRTP1,M)
  258.          EI(NRTP1) = XI(NRTP1,M)
  259.   220 CONTINUE
  260.       ER(M) = 0.0D0
  261.       EI(M) = 0.0D0
  262. C
  263. C     IF REQUIRED, GENERATE U.
  264. C
  265.       IF (.NOT.WANTU) GO TO 350
  266.          IF (NCU .LT. NCTP1) GO TO 250
  267.          DO 240 J = NCTP1, NCU
  268.             DO 230 I = 1, N
  269.                UR(I,J) = 0.0D0
  270.                UI(I,J) = 0.0D0
  271.   230       CONTINUE
  272.             UR(J,J) = 1.0D0
  273.             UI(J,J) = 0.0D0
  274.   240    CONTINUE
  275.   250    CONTINUE
  276.          IF (NCT .LT. 1) GO TO 340
  277.          DO 330 LL = 1, NCT
  278.             L = NCT - LL + 1
  279.             IF (CABS1(SR(L),SI(L)) .EQ. 0.0D0) GO TO 300
  280.                LP1 = L + 1
  281.                IF (NCU .LT. LP1) GO TO 270
  282.                DO 260 J = LP1, NCU
  283.                   TR = -WDOTCR(N-L+1,UR(L,L),UI(L,L),1,UR(L,J),
  284.      *                         UI(L,J),1)
  285.                   TI = -WDOTCI(N-L+1,UR(L,L),UI(L,L),1,UR(L,J),
  286.      *                         UI(L,J),1)
  287.                   CALL WDIV(TR,TI,UR(L,L),UI(L,L),TR,TI)
  288.                   CALL WAXPY(N-L+1,TR,TI,UR(L,L),UI(L,L),1,UR(L,J),
  289.      *                       UI(L,J),1)
  290.   260          CONTINUE
  291.   270          CONTINUE
  292.                CALL WRSCAL(N-L+1,-1.0D0,UR(L,L),UI(L,L),1)
  293.                UR(L,L) = FLOP(1.0D0 + UR(L,L))
  294.                LM1 = L - 1
  295.                IF (LM1 .LT. 1) GO TO 290
  296.                DO 280 I = 1, LM1
  297.                   UR(I,L) = 0.0D0
  298.                   UI(I,L) = 0.0D0
  299.   280          CONTINUE
  300.   290          CONTINUE
  301.             GO TO 320
  302.   300       CONTINUE
  303.                DO 310 I = 1, N
  304.                   UR(I,L) = 0.0D0
  305.                   UI(I,L) = 0.0D0
  306.   310          CONTINUE
  307.                UR(L,L) = 1.0D0
  308.                UI(L,L) = 0.0D0
  309.   320       CONTINUE
  310.   330    CONTINUE
  311.   340    CONTINUE
  312.   350 CONTINUE
  313. C
  314. C     IF IT IS REQUIRED, GENERATE V.
  315. C
  316.       IF (.NOT.WANTV) GO TO 400
  317.          DO 390 LL = 1, P
  318.             L = P - LL + 1
  319.             LP1 = L + 1
  320.             IF (L .GT. NRT) GO TO 370
  321.             IF (CABS1(ER(L),EI(L)) .EQ. 0.0D0) GO TO 370
  322.                DO 360 J = LP1, P
  323.                   TR = -WDOTCR(P-L,VR(LP1,L),VI(LP1,L),1,VR(LP1,J),
  324.      *                         VI(LP1,J),1)
  325.                   TI = -WDOTCI(P-L,VR(LP1,L),VI(LP1,L),1,VR(LP1,J),
  326.      *                         VI(LP1,J),1)
  327.                   CALL WDIV(TR,TI,VR(LP1,L),VI(LP1,L),TR,TI)
  328.                   CALL WAXPY(P-L,TR,TI,VR(LP1,L),VI(LP1,L),1,VR(LP1,J),
  329.      *                       VI(LP1,J),1)
  330.   360          CONTINUE
  331.   370       CONTINUE
  332.             DO 380 I = 1, P
  333.                VR(I,L) = 0.0D0
  334.                VI(I,L) = 0.0D0
  335.   380       CONTINUE
  336.             VR(L,L) = 1.0D0
  337.             VI(L,L) = 0.0D0
  338.   390    CONTINUE
  339.   400 CONTINUE
  340. C
  341. C     TRANSFORM S AND E SO THAT THEY ARE REAL.
  342. C
  343.       DO 420 I = 1, M
  344.             TR = PYTHAG(SR(I),SI(I))
  345.             IF (TR .EQ. 0.0D0) GO TO 405
  346.             RR = SR(I)/TR
  347.             RI = SI(I)/TR
  348.             SR(I) = TR
  349.             SI(I) = 0.0D0
  350.             IF (I .LT. M) CALL WDIV(ER(I),EI(I),RR,RI,ER(I),EI(I))
  351.             IF (WANTU) CALL WSCAL(N,RR,RI,UR(1,I),UI(1,I),1)
  352.   405    CONTINUE
  353. C     ...EXIT
  354.          IF (I .EQ. M) GO TO 430
  355.             TR = PYTHAG(ER(I),EI(I))
  356.             IF (TR .EQ. 0.0D0) GO TO 410
  357.             CALL WDIV(TR,0.0D0,ER(I),EI(I),RR,RI)
  358.             ER(I) = TR
  359.             EI(I) = 0.0D0
  360.             CALL WMUL(SR(I+1),SI(I+1),RR,RI,SR(I+1),SI(I+1))
  361.             IF (WANTV) CALL WSCAL(P,RR,RI,VR(1,I+1),VI(1,I+1),1)
  362.   410    CONTINUE
  363.   420 CONTINUE
  364.   430 CONTINUE
  365. C
  366. C     MAIN ITERATION LOOP FOR THE SINGULAR VALUES.
  367. C
  368.       MM = M
  369.       ITER = 0
  370.   440 CONTINUE
  371. C
  372. C        QUIT IF ALL THE SINGULAR VALUES HAVE BEEN FOUND.
  373. C
  374. C     ...EXIT
  375.          IF (M .EQ. 0) GO TO 700
  376. C
  377. C        IF TOO MANY ITERATIONS HAVE BEEN PERFORMED, SET
  378. C        FLAG AND RETURN.
  379. C
  380.          IF (ITER .LT. MAXIT) GO TO 450
  381.             INFO = M
  382. C     ......EXIT
  383.             GO TO 700
  384.   450    CONTINUE
  385. C
  386. C        THIS SECTION OF THE PROGRAM INSPECTS FOR
  387. C        NEGLIGIBLE ELEMENTS IN THE S AND E ARRAYS.  ON
  388. C        COMPLETION THE VARIABLE KASE IS SET AS FOLLOWS.
  389. C
  390. C           KASE = 1     IF SR(M) AND ER(L-1) ARE NEGLIGIBLE AND L.LT.M
  391. C           KASE = 2     IF SR(L) IS NEGLIGIBLE AND L.LT.M
  392. C           KASE = 3     IF ER(L-1) IS NEGLIGIBLE, L.LT.M, AND
  393. C                        SR(L), ..., SR(M) ARE NOT NEGLIGIBLE (QR STEP).
  394. C           KASE = 4     IF ER(M-1) IS NEGLIGIBLE (CONVERGENCE).
  395. C
  396.          DO 470 LL = 1, M
  397.             L = M - LL
  398. C        ...EXIT
  399.             IF (L .EQ. 0) GO TO 480
  400.             TEST = FLOP(DABS(SR(L)) + DABS(SR(L+1)))
  401.             ZTEST = FLOP(TEST + DABS(ER(L))/2.0D0)
  402.             IF (SMALL*ZTEST .NE. SMALL*TEST) GO TO 460
  403.                ER(L) = 0.0D0
  404. C        ......EXIT
  405.                GO TO 480
  406.   460       CONTINUE
  407.   470    CONTINUE
  408.   480    CONTINUE
  409.          IF (L .NE. M - 1) GO TO 490
  410.             KASE = 4
  411.          GO TO 560
  412.   490    CONTINUE
  413.             LP1 = L + 1
  414.             MP1 = M + 1
  415.             DO 510 LLS = LP1, MP1
  416.                LS = M - LLS + LP1
  417. C           ...EXIT
  418.                IF (LS .EQ. L) GO TO 520
  419.                TEST = 0.0D0
  420.                IF (LS .NE. M) TEST = FLOP(TEST + DABS(ER(LS)))
  421.                IF (LS .NE. L + 1) TEST = FLOP(TEST + DABS(ER(LS-1)))
  422.                ZTEST = FLOP(TEST + DABS(SR(LS))/2.0D0)
  423.                IF (SMALL*ZTEST .NE. SMALL*TEST) GO TO 500
  424.                   SR(LS) = 0.0D0
  425. C           ......EXIT
  426.                   GO TO 520
  427.   500          CONTINUE
  428.   510       CONTINUE
  429.   520       CONTINUE
  430.             IF (LS .NE. L) GO TO 530
  431.                KASE = 3
  432.             GO TO 550
  433.   530       CONTINUE
  434.             IF (LS .NE. M) GO TO 540
  435.                KASE = 1
  436.             GO TO 550
  437.   540       CONTINUE
  438.                KASE = 2
  439.                L = LS
  440.   550       CONTINUE
  441.   560    CONTINUE
  442.          L = L + 1
  443. C
  444. C        PERFORM THE TASK INDICATED BY KASE.
  445. C
  446.          GO TO (570, 600, 620, 650), KASE
  447. C
  448. C        DEFLATE NEGLIGIBLE SR(M).
  449. C
  450.   570    CONTINUE
  451.             MM1 = M - 1
  452.             F = ER(M-1)
  453.             ER(M-1) = 0.0D0
  454.             DO 590 KK = L, MM1
  455.                K = MM1 - KK + L
  456.                T1 = SR(K)
  457.                CALL RROTG(T1,F,CS,SN)
  458.                SR(K) = T1
  459.                IF (K .EQ. L) GO TO 580
  460.                   F = FLOP(-SN*ER(K-1))
  461.                   ER(K-1) = FLOP(CS*ER(K-1))
  462.   580          CONTINUE
  463.                IF (WANTV) CALL RROT(P,VR(1,K),1,VR(1,M),1,CS,SN)
  464.                IF (WANTV) CALL RROT(P,VI(1,K),1,VI(1,M),1,CS,SN)
  465.   590       CONTINUE
  466.          GO TO 690
  467. C
  468. C        SPLIT AT NEGLIGIBLE SR(L).
  469. C
  470.   600    CONTINUE
  471.             F = ER(L-1)
  472.             ER(L-1) = 0.0D0
  473.             DO 610 K = L, M
  474.                T1 = SR(K)
  475.                CALL RROTG(T1,F,CS,SN)
  476.                SR(K) = T1
  477.                F = FLOP(-SN*ER(K))
  478.                ER(K) = FLOP(CS*ER(K))
  479.                IF (WANTU) CALL RROT(N,UR(1,K),1,UR(1,L-1),1,CS,SN)
  480.                IF (WANTU) CALL RROT(N,UI(1,K),1,UI(1,L-1),1,CS,SN)
  481.   610       CONTINUE
  482.          GO TO 690
  483. C
  484. C        PERFORM ONE QR STEP.
  485. C
  486.   620    CONTINUE
  487. C
  488. C           CALCULATE THE SHIFT.
  489. C
  490.             SCALE = DMAX1(DABS(SR(M)),DABS(SR(M-1)),DABS(ER(M-1)),
  491.      *                    DABS(SR(L)),DABS(ER(L)))
  492.             SM = SR(M)/SCALE
  493.             SMM1 = SR(M-1)/SCALE
  494.             EMM1 = ER(M-1)/SCALE
  495.             SL = SR(L)/SCALE
  496.             EL = ER(L)/SCALE
  497.             B = FLOP(((SMM1 + SM)*(SMM1 - SM) + EMM1**2)/2.0D0)
  498.             C = FLOP((SM*EMM1)**2)
  499.             SHIFT = 0.0D0
  500.             IF (B .EQ. 0.0D0 .AND. C .EQ. 0.0D0) GO TO 630
  501.                SHIFT = FLOP(DSQRT(B**2+C))
  502.                IF (B .LT. 0.0D0) SHIFT = -SHIFT
  503.                SHIFT = FLOP(C/(B + SHIFT))
  504.   630       CONTINUE
  505.             F = FLOP((SL + SM)*(SL - SM) - SHIFT)
  506.             G = FLOP(SL*EL)
  507. C
  508. C           CHASE ZEROS.
  509. C
  510.             MM1 = M - 1
  511.             DO 640 K = L, MM1
  512.                CALL RROTG(F,G,CS,SN)
  513.                IF (K .NE. L) ER(K-1) = F
  514.                F = FLOP(CS*SR(K) + SN*ER(K))
  515.                ER(K) = FLOP(CS*ER(K) - SN*SR(K))
  516.                G = FLOP(SN*SR(K+1))
  517.                SR(K+1) = FLOP(CS*SR(K+1))
  518.                IF (WANTV) CALL RROT(P,VR(1,K),1,VR(1,K+1),1,CS,SN)
  519.                IF (WANTV) CALL RROT(P,VI(1,K),1,VI(1,K+1),1,CS,SN)
  520.                CALL RROTG(F,G,CS,SN)
  521.                SR(K) = F
  522.                F = FLOP(CS*ER(K) + SN*SR(K+1))
  523.                SR(K+1) = FLOP(-SN*ER(K) + CS*SR(K+1))
  524.                G = FLOP(SN*ER(K+1))
  525.                ER(K+1) = FLOP(CS*ER(K+1))
  526.                IF (WANTU .AND. K .LT. N)
  527.      *            CALL RROT(N,UR(1,K),1,UR(1,K+1),1,CS,SN)
  528.                IF (WANTU .AND. K .LT. N)
  529.      *            CALL RROT(N,UI(1,K),1,UI(1,K+1),1,CS,SN)
  530.   640       CONTINUE
  531.             ER(M-1) = F
  532.             ITER = ITER + 1
  533.          GO TO 690
  534. C
  535. C        CONVERGENCE
  536. C
  537.   650    CONTINUE
  538. C
  539. C           MAKE THE SINGULAR VALUE  POSITIVE
  540. C
  541.             IF (SR(L) .GE. 0.0D0) GO TO 660
  542.                SR(L) = -SR(L)
  543.              IF (WANTV) CALL WRSCAL(P,-1.0D0,VR(1,L),VI(1,L),1)
  544.   660       CONTINUE
  545. C
  546. C           ORDER THE SINGULAR VALUE.
  547. C
  548.   670       IF (L .EQ. MM) GO TO 680
  549. C           ...EXIT
  550.                IF (SR(L) .GE. SR(L+1)) GO TO 680
  551.                TR = SR(L)
  552.                SR(L) = SR(L+1)
  553.                SR(L+1) = TR
  554.                IF (WANTV .AND. L .LT. P)
  555.      *            CALL WSWAP(P,VR(1,L),VI(1,L),1,VR(1,L+1),VI(1,L+1),1)
  556.                IF (WANTU .AND. L .LT. N)
  557.      *            CALL WSWAP(N,UR(1,L),UI(1,L),1,UR(1,L+1),UI(1,L+1),1)
  558.                L = L + 1
  559.             GO TO 670
  560.   680       CONTINUE
  561.             ITER = 0
  562.             M = M - 1
  563.   690    CONTINUE
  564.       GO TO 440
  565.   700 CONTINUE
  566.       RETURN
  567.       END
  568.       SUBROUTINE WQRDC(XR,XI,LDX,N,P,QRAUXR,QRAUXI,JPVT,WORKR,WORKI,
  569.      *                 JOB)
  570.       INTEGER LDX,N,P,JOB
  571.       INTEGER JPVT(1)
  572.       DOUBLE PRECISION XR(LDX,1),XI(LDX,1),QRAUXR(1),QRAUXI(1),
  573.      *                 WORKR(1),WORKI(1)
  574. C
  575. C     WQRDC USES HOUSEHOLDER TRANSFORMATIONS TO COMPUTE THE QR
  576. C     FACTORIZATION OF AN N BY P MATRIX X.  COLUMN PIVOTING
  577. C     BASED ON THE 2-NORMS OF THE REDUCED COLUMNS MAY BE
  578. C     PERFORMED AT THE USERS OPTION.
  579. C
  580. C     ON ENTRY
  581. C
  582. C        X       DOUBLE-COMPLEX(LDX,P), WHERE LDX .GE. N.
  583. C                X CONTAINS THE MATRIX WHOSE DECOMPOSITION IS TO BE
  584. C                COMPUTED.
  585. C
  586. C        LDX     INTEGER.
  587. C                LDX IS THE LEADING DIMENSION OF THE ARRAY X.
  588. C
  589. C        N       INTEGER.
  590. C                N IS THE NUMBER OF ROWS OF THE MATRIX X.
  591. C
  592. C        P       INTEGER.
  593. C                P IS THE NUMBER OF COLUMNS OF THE MATRIX X.
  594. C
  595. C        JPVT    INTEGER(P).
  596. C                JPVT CONTAINS INTEGERS THAT CONTROL THE SELECTION
  597. C                OF THE PIVOT COLUMNS.  THE K-TH COLUMN X(K) OF X
  598. C                IS PLACED IN ONE OF THREE CLASSES ACCORDING TO THE
  599. C                VALUE OF JPVT(K).
  600. C
  601. C                   IF JPVT(K) .GT. 0, THEN X(K) IS AN INITIAL
  602. C                                      COLUMN.
  603. C
  604. C                   IF JPVT(K) .EQ. 0, THEN X(K) IS A FREE COLUMN.
  605. C
  606. C                   IF JPVT(K) .LT. 0, THEN X(K) IS A FINAL COLUMN.
  607. C
  608. C                BEFORE THE DECOMPOSITION IS COMPUTED, INITIAL COLUMNS
  609. C                ARE MOVED TO THE BEGINNING OF THE ARRAY X AND FINAL
  610. C                COLUMNS TO THE END.  BOTH INITIAL AND FINAL COLUMNS
  611. C                ARE FROZEN IN PLACE DURING THE COMPUTATION AND ONLY
  612. C                FREE COLUMNS ARE MOVED.  AT THE K-TH STAGE OF THE
  613. C                REDUCTION, IF X(K) IS OCCUPIED BY A FREE COLUMN
  614. C                IT IS INTERCHANGED WITH THE FREE COLUMN OF LARGEST
  615. C                REDUCED NORM.  JPVT IS NOT REFERENCED IF
  616. C                JOB .EQ. 0.
  617. C
  618. C        WORK    DOUBLE-COMPLEX(P).
  619. C                WORK IS A WORK ARRAY.  WORK IS NOT REFERENCED IF
  620. C                JOB .EQ. 0.
  621. C
  622. C        JOB     INTEGER.
  623. C                JOB IS AN INTEGER THAT INITIATES COLUMN PIVOTING.
  624. C                IF JOB .EQ. 0, NO PIVOTING IS DONE.
  625. C                IF JOB .NE. 0, PIVOTING IS DONE.
  626. C
  627. C     ON RETURN
  628. C
  629. C        X       X CONTAINS IN ITS UPPER TRIANGLE THE UPPER
  630. C                TRIANGULAR MATRIX R OF THE QR FACTORIZATION.
  631. C                BELOW ITS DIAGONAL X CONTAINS INFORMATION FROM
  632. C                WHICH THE UNITARY PART OF THE DECOMPOSITION
  633. C                CAN BE RECOVERED.  NOTE THAT IF PIVOTING HAS
  634. C                BEEN REQUESTED, THE DECOMPOSITION IS NOT THAT
  635. C                OF THE ORIGINAL MATRIX X BUT THAT OF X
  636. C                WITH ITS COLUMNS PERMUTED AS DESCRIBED BY JPVT.
  637. C
  638. C        QRAUX   DOUBLE-COMPLEX(P).
  639. C                QRAUX CONTAINS FURTHER INFORMATION REQUIRED TO RECOVER
  640. C                THE UNITARY PART OF THE DECOMPOSITION.
  641. C
  642. C        JPVT    JPVT(K) CONTAINS THE INDEX OF THE COLUMN OF THE
  643. C                ORIGINAL MATRIX THAT HAS BEEN INTERCHANGED INTO
  644. C                THE K-TH COLUMN, IF PIVOTING WAS REQUESTED.
  645. C
  646. C     LINPACK. THIS VERSION DATED 07/03/79 .
  647. C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
  648. C
  649. C     WQRDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.
  650. C
  651. C     BLAS WAXPY,PYTHAG,WDOTCR,WDOTCI,WSCAL,WSWAP,WNRM2
  652. C     FORTRAN DABS,DIMAG,DMAX1,MIN0
  653. C
  654. C     INTERNAL VARIABLES
  655. C
  656.       INTEGER J,JP,L,LP1,LUP,MAXJ,PL,PU
  657.       DOUBLE PRECISION MAXNRM,WNRM2,TT
  658.       DOUBLE PRECISION PYTHAG,WDOTCR,WDOTCI,NRMXLR,NRMXLI,TR,TI,FLOP
  659.       LOGICAL NEGJ,SWAPJ
  660. C
  661.       DOUBLE PRECISION ZDUMR,ZDUMI
  662.       DOUBLE PRECISION CABS1
  663. C
  664.       PL = 1
  665.       PU = 0
  666.       IF (JOB .EQ. 0) GO TO 60
  667. C
  668. C        PIVOTING HAS BEEN REQUESTED.  REARRANGE THE COLUMNS
  669. C        ACCORDING TO JPVT.
  670. C
  671.          DO 20 J = 1, P
  672.             SWAPJ = JPVT(J) .GT. 0
  673.             NEGJ = JPVT(J) .LT. 0
  674.             JPVT(J) = J
  675.             IF (NEGJ) JPVT(J) = -J
  676.             IF (.NOT.SWAPJ) GO TO 10
  677.                IF (J .NE. PL)
  678.      *            CALL WSWAP(N,XR(1,PL),XI(1,PL),1,XR(1,J),XI(1,J),1)
  679.                JPVT(J) = JPVT(PL)
  680.                JPVT(PL) = J
  681.                PL = PL + 1
  682.    10       CONTINUE
  683.    20    CONTINUE
  684.          PU = P
  685.          DO 50 JJ = 1, P
  686.             J = P - JJ + 1
  687.             IF (JPVT(J) .GE. 0) GO TO 40
  688.                JPVT(J) = -JPVT(J)
  689.                IF (J .EQ. PU) GO TO 30
  690.                   CALL WSWAP(N,XR(1,PU),XI(1,PU),1,XR(1,J),XI(1,J),1)
  691.                   JP = JPVT(PU)
  692.                   JPVT(PU) = JPVT(J)
  693.                   JPVT(J) = JP
  694.    30          CONTINUE
  695.                PU = PU - 1
  696.    40       CONTINUE
  697.    50    CONTINUE
  698.    60 CONTINUE
  699. C
  700. C     COMPUTE THE NORMS OF THE FREE COLUMNS.
  701. C
  702.       IF (PU .LT. PL) GO TO 80
  703.       DO 70 J = PL, PU
  704.          QRAUXR(J) = WNRM2(N,XR(1,J),XI(1,J),1)
  705.          QRAUXI(J) = 0.0D0
  706.          WORKR(J) = QRAUXR(J)
  707.          WORKI(J) = QRAUXI(J)
  708.    70 CONTINUE
  709.    80 CONTINUE
  710. C
  711. C     PERFORM THE HOUSEHOLDER REDUCTION OF X.
  712. C
  713.       LUP = MIN0(N,P)
  714.       DO 210 L = 1, LUP
  715.          IF (L .LT. PL .OR. L .GE. PU) GO TO 120
  716. C
  717. C           LOCATE THE COLUMN OF LARGEST NORM AND BRING IT
  718. C           INTO THE PIVOT POSITION.
  719. C
  720.             MAXNRM = 0.0D0
  721.             MAXJ = L
  722.             DO 100 J = L, PU
  723.                IF (QRAUXR(J) .LE. MAXNRM) GO TO 90
  724.                   MAXNRM = QRAUXR(J)
  725.                   MAXJ = J
  726.    90          CONTINUE
  727.   100       CONTINUE
  728.             IF (MAXJ .EQ. L) GO TO 110
  729.                CALL WSWAP(N,XR(1,L),XI(1,L),1,XR(1,MAXJ),XI(1,MAXJ),1)
  730.                QRAUXR(MAXJ) = QRAUXR(L)
  731.                QRAUXI(MAXJ) = QRAUXI(L)
  732.                WORKR(MAXJ) = WORKR(L)
  733.                WORKI(MAXJ) = WORKI(L)
  734.                JP = JPVT(MAXJ)
  735.                JPVT(MAXJ) = JPVT(L)
  736.                JPVT(L) = JP
  737.   110       CONTINUE
  738.   120    CONTINUE
  739.          QRAUXR(L) = 0.0D0
  740.          QRAUXI(L) = 0.0D0
  741.          IF (L .EQ. N) GO TO 200
  742. C
  743. C           COMPUTE THE HOUSEHOLDER TRANSFORMATION FOR COLUMN L.
  744. C
  745.             NRMXLR = WNRM2(N-L+1,XR(L,L),XI(L,L),1)
  746.             NRMXLI = 0.0D0
  747.             IF (CABS1(NRMXLR,NRMXLI) .EQ. 0.0D0) GO TO 190
  748.                IF (CABS1(XR(L,L),XI(L,L)) .EQ. 0.0D0) GO TO 130
  749.                  CALL WSIGN(NRMXLR,NRMXLI,XR(L,L),XI(L,L),NRMXLR,NRMXLI)
  750.   130          CONTINUE
  751.                CALL WDIV(1.0D0,0.0D0,NRMXLR,NRMXLI,TR,TI)
  752.                CALL WSCAL(N-L+1,TR,TI,XR(L,L),XI(L,L),1)
  753.                XR(L,L) = FLOP(1.0D0 + XR(L,L))
  754. C
  755. C              APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS,
  756. C              UPDATING THE NORMS.
  757. C
  758.                LP1 = L + 1
  759.                IF (P .LT. LP1) GO TO 180
  760.                DO 170 J = LP1, P
  761.                   TR = -WDOTCR(N-L+1,XR(L,L),XI(L,L),1,XR(L,J),
  762.      *                         XI(L,J),1)
  763.                   TI = -WDOTCI(N-L+1,XR(L,L),XI(L,L),1,XR(L,J),
  764.      *                         XI(L,J),1)
  765.                   CALL WDIV(TR,TI,XR(L,L),XI(L,L),TR,TI)
  766.                   CALL WAXPY(N-L+1,TR,TI,XR(L,L),XI(L,L),1,XR(L,J),
  767.      *                       XI(L,J),1)
  768.                   IF (J .LT. PL .OR. J .GT. PU) GO TO 160
  769.                   IF (CABS1(QRAUXR(J),QRAUXI(J)) .EQ. 0.0D0)
  770.      *               GO TO 160
  771.                      TT = 1.0D0 - (PYTHAG(XR(L,J),XI(L,J))/QRAUXR(J))**2
  772.                      TT = DMAX1(TT,0.0D0)
  773.                      TR = FLOP(TT)
  774.                      TT = FLOP(1.0D0+0.05D0*TT*(QRAUXR(J)/WORKR(J))**2)
  775.                      IF (TT .EQ. 1.0D0) GO TO 140
  776.                         QRAUXR(J) = QRAUXR(J)*DSQRT(TR)
  777.                         QRAUXI(J) = QRAUXI(J)*DSQRT(TR)
  778.                      GO TO 150
  779.   140                CONTINUE
  780.                         QRAUXR(J) = WNRM2(N-L,XR(L+1,J),XI(L+1,J),1)
  781.                         QRAUXI(J) = 0.0D0
  782.                         WORKR(J) = QRAUXR(J)
  783.                         WORKI(J) = QRAUXI(J)
  784.   150                CONTINUE
  785.   160             CONTINUE
  786.   170          CONTINUE
  787.   180          CONTINUE
  788. C
  789. C              SAVE THE TRANSFORMATION.
  790. C
  791.                QRAUXR(L) = XR(L,L)
  792.                QRAUXI(L) = XI(L,L)
  793.                XR(L,L) = -NRMXLR
  794.                XI(L,L) = -NRMXLI
  795.   190       CONTINUE
  796.   200    CONTINUE
  797.   210 CONTINUE
  798.       RETURN
  799.       END
  800.