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 / COMQR3.FOR < prev    next >
Encoding:
Text File  |  1991-04-13  |  12.2 KB  |  408 lines
  1.       SUBROUTINE COMQR3 (NM, N, LOW, IGH, ORTR, ORTI, HR, HI,
  2.      .                   WR, WI, ZR, ZI, IERR, JOB)
  3.       IMPLICIT NONE
  4. C
  5. C A TRANSLATION OF A UNITARY ANALOGUE OF THE ALGOL PROCEDURE COMLR2,
  6. C   NUM. MATH. 16, 181-204 (1970) BY PETERS AND WILKINSON.
  7. C   HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395 (1971).
  8. C   THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS
  9. C   (COMP. JOUR. 4, 332-345 (1962)) FOR THE LR ALGORITHM.
  10. C
  11. C FINDS THE EIGENVALUES AND EIGENVECTORS OF A COMPLEX UPPER HESSENBERG
  12. C   MATRIX BY THE QR METHOD.  THE EIGENVECTORS OF A COMPLEX GENERAL MATRIX
  13. C   CAN ALSO BE FOUND IF CORTH HAS BEEN USED TO REDUCE THIS GENERAL MATRIX
  14. C   TO HESSENBERG FORM.
  15. C
  16. C ON INPUT:
  17. C   NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
  18. C      AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT.
  19. C
  20. C   N IS THE ORDER OF THE MATRIX.
  21. C
  22. C   LOW & IGH ARE INTEGERS DETERMINED BY THE BALANCING SUBROUTINE CBAL.
  23. C     IF CBAL HAS NOT BEEN USED, SET LOW = 1, IGH = N.
  24. C
  25. C   ORTR & ORTI CONTAIN INFORMATION ABOUT THE UNITARY TRANSFORMATIONS USED
  26. C     IN THE REDUCTION BY CORTH, IF PERFORMED.  ONLY ELEMENTS LOW THROUGH
  27. C     IGH ARE USED.  IF THE EIGENVECTORS OF THE HESSENBERG MATRIX ARE
  28. C     DESIRED, SET ORTR(J) AND ORTI(J) TO 0.0D0 FOR THESE ELEMENTS.
  29. C
  30. C   HR & HI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY, OF THE
  31. C     COMPLEX UPPER HESSENBERG MATRIX.  THEIR LOWER TRIANGLES BELOW THE
  32. C     SUBDIAGONAL CONTAIN FURTHER INFORMATION ABOUT THE TRANSFORMATIONS
  33. C     WHICH WERE USED IN THE REDUCTION BY CORTH, IF PERFORMED.  IF THE
  34. C     EIGENVECTORS OF THE HESSENBERG MATRIX ARE DESIRED, THESE ELEMENTS
  35. C     MAY BE ARBITRARY.
  36. C
  37. C   JOB = 0  OUTPUT H = SCHUR TRIANGULAR FORM, Z NOT USED
  38. C       = 1  OUTPUT H = SCHUR FORM, Z = UNITARY SIMILARITY
  39. C       = 2  SAME AS COMQR2
  40. C       = 3  OUTPUT H = HESSENBERG FORM, Z = UNITARY SIMILARITY
  41. C
  42. C ON OUTPUT:
  43. C   ORTR, ORTI, & THE UPPER HESSENBERG PORTIONS OF HR & HI HAVE BEEN TRASHED.
  44. C
  45. C   WR & WI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY, OF THE
  46. C     EIGENVALUES.  IF AN ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE
  47. C     CORRECT FOR INDICES IERR+1,...,N.
  48. C
  49. C   ZR & ZI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY, OF THE
  50. C     EIGENVECTORS.  THE EIGENVECTORS ARE UNNORMALIZED.  IF AN ERROR EXIT
  51. C     IS MADE, NONE OF THE EIGENVECTORS HAS BEEN FOUND.
  52. C
  53. C   IERR IS SET TO:
  54. C     ZERO FOR NORMAL RETURN,
  55. C     J    IF J-TH EIGENVALUE HAS NOT BEEN DETERMINED AFTER 30*N ITERATIONS.
  56. C
  57. C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
  58. C   APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  59. C
  60.       INTEGER NM, N, LOW, IGH, IERR, JOB
  61.       DOUBLE PRECISION ORTR(IGH), ORTI(IGH), HR(NM,N), HI(NM,N),
  62.      .                 WR(N), WI(N), ZR(NM,N), ZI(NM,N)
  63. C
  64.       INTEGER I, J, K, L, M, EN, II, JJ, LL, NN, IP1, ITN, ITS,
  65.      .        LP1, ENM1, IEND
  66.       DOUBLE PRECISION SI, SR, TI, TR, XI, XR, YI, YR, ZZI, ZZR, NORM
  67. C
  68.       DOUBLE PRECISION FLOP, PYTHAG
  69. C
  70. C
  71.       IERR = 0
  72.       IF (JOB.EQ.0) GO TO 150
  73. C
  74. C ***      INITIALIZE EIGENVECTOR MATRIX
  75.       DO 101 I = 1, N
  76.         DO 100 J = 1, N
  77.           ZR(I,J) = 0.0D0
  78.           ZI(I,J) = 0.0D0
  79.           IF (I.EQ.J) ZR(I,J) = 1.0D0
  80. 100     CONTINUE
  81. 101   CONTINUE
  82. C
  83. C ***      FORM THE MATRIX OF ACCUMULATED TRANSFORMATIONS
  84. C           FROM THE INFORMATION LEFT BY CORTH
  85.       IEND = IGH-LOW-1
  86.       IF (IEND) 180, 150, 105
  87. C
  88. C ***      FOR I = IGH-1 STEP -1 UNTIL LOW+1 DO
  89. 105   CONTINUE
  90.       DO 140 II = 1, IEND
  91.         I = IGH-II
  92.         IF (ORTR(I).EQ.0.0D0 .AND. ORTI(I).EQ.0.0D0) GO TO 140
  93.         IF (HR(I,I-1).EQ.0.0D0 .AND. HI(I,I-1).EQ.0.0D0) GO TO 140
  94. C
  95. C ***      NORM BELOW IS NEGATIVE OF H FORMED IN CORTH
  96.         NORM = FLOP (HR(I,I-1)*ORTR(I)+HI(I,I-1)*ORTI(I))
  97.         IP1 = I+1
  98.         DO 110 K = IP1, IGH
  99.           ORTR(K) = HR(K,I-1)
  100.           ORTI(K) = HI(K,I-1)
  101. 110     CONTINUE
  102.         DO 130 J = I, IGH
  103.           SR = 0.0D0
  104.           SI = 0.0D0
  105.           DO 115 K = I, IGH
  106.             SR = FLOP (SR+ORTR(K)*ZR(K,J)+ORTI(K)*ZI(K,J))
  107.             SI = FLOP (SI+ORTR(K)*ZI(K,J)-ORTI(K)*ZR(K,J))
  108. 115       CONTINUE
  109.           SR = FLOP (SR/NORM)
  110.           SI = FLOP (SI/NORM)
  111.           DO 120 K = I, IGH
  112.             ZR(K,J) = FLOP (ZR(K,J)+SR*ORTR(K)-SI*ORTI(K))
  113.             ZI(K,J) = FLOP (ZI(K,J)+SR*ORTI(K)+SI*ORTR(K))
  114. 120       CONTINUE
  115. 130     CONTINUE
  116. 140   CONTINUE
  117. C
  118.       IF (JOB.EQ.3) GO TO 1001
  119. C
  120. C ***      CREATE REAL SUBDIAGONAL ELEMENTS
  121. 150   CONTINUE
  122.       L = LOW+1
  123.       DO 170 I = L, IGH
  124.         LL = MIN0 (I+1, IGH)
  125.         IF (HI(I,I-1).EQ.0.0D0) GO TO 170
  126.         NORM = PYTHAG (HR(I,I-1), HI(I,I-1))
  127.         YR = FLOP (HR(I,I-1)/NORM)
  128.         YI = FLOP (HI(I,I-1)/NORM)
  129.         HR(I,I-1) = NORM
  130.         HI(I,I-1) = 0.0D0
  131.         DO 155 J = I, N
  132.           SI = FLOP (YR*HI(I,J)-YI*HR(I,J))
  133.           HR(I,J) = FLOP (YR*HR(I,J)+YI*HI(I,J))
  134.           HI(I,J) = SI
  135. 155     CONTINUE
  136.         DO 160 J = 1, LL
  137.           SI = FLOP (YR*HI(J,I)+YI*HR(J,I))
  138.           HR(J,I) = FLOP (YR*HR(J,I)-YI*HI(J,I))
  139.           HI(J,I) = SI
  140. 160     CONTINUE
  141.         IF (JOB.EQ.0) GO TO 170
  142.         DO 165 J = LOW, IGH
  143.           SI = FLOP (YR*ZI(J,I)+YI*ZR(J,I))
  144.           ZR(J,I) = FLOP (YR*ZR(J,I)-YI*ZI(J,I))
  145.           ZI(J,I) = SI
  146. 165     CONTINUE
  147. 170   CONTINUE
  148. C
  149. C ***      STORE ROOTS ISOLATED BY CBAL
  150. 180   CONTINUE
  151.       DO 200 I = 1, N
  152.         IF (I.GE.LOW .AND. I.LE.IGH) GO TO 200
  153.         WR(I) = HR(I,I)
  154.         WI(I) = HI(I,I)
  155. 200   CONTINUE
  156.       EN = IGH
  157.       TR = 0.0D0
  158.       TI = 0.0D0
  159.       ITN = 30*N
  160. C
  161. C ***      SEARCH FOR NEXT EIGENVALUE
  162. 220   CONTINUE
  163.       IF (EN.LT.LOW) GO TO 680
  164.       ITS = 0
  165.       ENM1 = EN-1
  166. C
  167. C ***      LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
  168. C            FOR L = EN STEP -1 UNTIL LOW DO
  169. 240   CONTINUE
  170.       DO 260 LL = LOW, EN
  171.         L = EN+LOW-LL
  172.         IF (L.EQ.LOW) GO TO 300
  173.         XR = FLOP (DABS (HR(L-1,L-1))+DABS (HI(L-1,L-1))+
  174.      .             DABS (HR(L,L))+DABS (HI(L,L)))
  175.         YR = FLOP (XR+DABS (HR(L,L-1)))
  176.         IF (XR.EQ.YR) GO TO 300
  177. 260   CONTINUE
  178. C
  179. C ***      FORM SHIFT
  180. 300   CONTINUE
  181.       IF (L.EQ.EN) GO TO 660
  182.       IF (ITN.EQ.0) GO TO 1000
  183.       IF (ITS.EQ.10 .OR. ITS.EQ.20) GO TO 320
  184.       SR = HR(EN,EN)
  185.       SI = HI(EN,EN)
  186.       XR = FLOP (HR(ENM1,EN)*HR(EN,ENM1))
  187.       XI = FLOP (HI(ENM1,EN)*HR(EN,ENM1))
  188.       IF (XR.EQ.0.0D0 .AND. XI.EQ.0.0D0) GO TO 340
  189.       YR = FLOP ((HR(ENM1,ENM1)-SR)/2.0D0)
  190.       YI = FLOP ((HI(ENM1,ENM1)-SI)/2.0D0)
  191.       CALL WSQRT (YR**2-YI**2+XR, 2.0D0*YR*YI+XI, ZZR, ZZI)
  192.       IF (YR*ZZR+YI*ZZI.GE.0.0D0) GO TO 310
  193.       ZZR = -ZZR
  194.       ZZI = -ZZI
  195. 310   CONTINUE
  196.       CALL WDIV (XR, XI, YR+ZZR, YI+ZZI, ZZR, ZZI)
  197.       SR = FLOP (SR-ZZR)
  198.       SI = FLOP (SI-ZZI)
  199.       GO TO 340
  200. C
  201. C ***      FORM EXCEPTIONAL SHIFT
  202. 320   CONTINUE
  203.       SR = FLOP (DABS (HR(EN,ENM1))+DABS (HR(ENM1,EN-2)))
  204.       SI = 0.0D0
  205. 340   CONTINUE
  206.       DO 360 I = LOW, EN
  207.         HR(I,I) = FLOP (HR(I,I)-SR)
  208.         HI(I,I) = FLOP (HI(I,I)-SI)
  209. 360   CONTINUE
  210.       TR = FLOP (TR+SR)
  211.       TI = FLOP (TI+SI)
  212.       ITS = ITS+1
  213.       ITN = ITN-1
  214. C
  215. C ***      REDUCE TO TRIANGLE (ROWS)
  216.       LP1 = L+1
  217.       DO 500 I = LP1, EN
  218.         SR = HR(I,I-1)
  219.         HR(I,I-1) = 0.0D0
  220.         NORM = FLOP (DABS (HR(I-1,I-1))+DABS (HI(I-1,I-1))+DABS (SR))
  221.         NORM = FLOP (NORM*DSQRT ((HR(I-1,I-1)/NORM)**2+
  222.      .               (HI(I-1,I-1)/NORM)**2+(SR/NORM)**2))
  223.         XR = FLOP (HR(I-1,I-1)/NORM)
  224.         WR(I-1) = XR
  225.         XI = FLOP (HI(I-1,I-1)/NORM)
  226.         WI(I-1) = XI
  227.         HR(I-1,I-1) = NORM
  228.         HI(I-1,I-1) = 0.0D0
  229.         HI(I,I-1) = FLOP (SR/NORM)
  230.         DO 490 J = I, N
  231.           YR = HR(I-1,J)
  232.           YI = HI(I-1,J)
  233.           ZZR = HR(I,J)
  234.           ZZI = HI(I,J)
  235.           HR(I-1,J) = FLOP (XR*YR+XI*YI+HI(I,I-1)*ZZR)
  236.           HI(I-1,J) = FLOP (XR*YI-XI*YR+HI(I,I-1)*ZZI)
  237.           HR(I,J) = FLOP (XR*ZZR-XI*ZZI-HI(I,I-1)*YR)
  238.           HI(I,J) = FLOP (XR*ZZI+XI*ZZR-HI(I,I-1)*YI)
  239. 490     CONTINUE
  240. 500   CONTINUE
  241. C
  242.       SI = HI(EN,EN)
  243.       IF (SI.EQ.0.0D0) GO TO 540
  244.       NORM = PYTHAG (HR(EN,EN), SI)
  245.       SR = FLOP (HR(EN,EN)/NORM)
  246.       SI = FLOP (SI/NORM)
  247.       HR(EN,EN) = NORM
  248.       HI(EN,EN) = 0.0D0
  249.       IF (EN.EQ.N) GO TO 540
  250.       IP1 = EN+1
  251.       DO 520 J = IP1, N
  252.         YR = HR(EN,J)
  253.         YI = HI(EN,J)
  254.         HR(EN,J) = FLOP (SR*YR+SI*YI)
  255.         HI(EN,J) = FLOP (SR*YI-SI*YR)
  256. 520   CONTINUE
  257. C
  258. C ***      INVERSE OPERATION (COLUMNS)
  259. 540   CONTINUE
  260.       DO 600 J = LP1, EN
  261.         XR = WR(J-1)
  262.         XI = WI(J-1)
  263.         DO 580 I = 1, J
  264.           YR = HR(I,J-1)
  265.           YI = 0.0D0
  266.           ZZR = HR(I,J)
  267.           ZZI = HI(I,J)
  268.           IF (I.EQ.J) GO TO 560
  269.           YI = HI(I,J-1)
  270.           HI(I,J-1) = FLOP (XR*YI+XI*YR+HI(J,J-1)*ZZI)
  271. 560       CONTINUE
  272.           HR(I,J-1) = FLOP (XR*YR-XI*YI+HI(J,J-1)*ZZR)
  273.           HR(I,J) = FLOP (XR*ZZR+XI*ZZI-HI(J,J-1)*YR)
  274.           HI(I,J) = FLOP (XR*ZZI-XI*ZZR-HI(J,J-1)*YI)
  275. 580     CONTINUE
  276.         IF (JOB.EQ.0) GO TO 600
  277.         DO 590 I = LOW, IGH
  278.           YR = ZR(I,J-1)
  279.           YI = ZI(I,J-1)
  280.           ZZR = ZR(I,J)
  281.           ZZI = ZI(I,J)
  282.           ZR(I,J-1) = FLOP (XR*YR-XI*YI+HI(J,J-1)*ZZR)
  283.           ZI(I,J-1) = FLOP (XR*YI+XI*YR+HI(J,J-1)*ZZI)
  284.           ZR(I,J) = FLOP (XR*ZZR+XI*ZZI-HI(J,J-1)*YR)
  285.           ZI(I,J) = FLOP (XR*ZZI-XI*ZZR-HI(J,J-1)*YI)
  286. 590     CONTINUE
  287. 600   CONTINUE
  288. C
  289.       IF (SI.EQ.0.0D0) GO TO 240
  290.       DO 630 I = 1, EN
  291.         YR = HR(I,EN)
  292.         YI = HI(I,EN)
  293.         HR(I,EN) = FLOP (SR*YR-SI*YI)
  294.         HI(I,EN) = FLOP (SR*YI+SI*YR)
  295. 630   CONTINUE
  296. C
  297.       IF (JOB.EQ.0) GO TO 240
  298.       DO 640 I = LOW, IGH
  299.         YR = ZR(I,EN)
  300.         YI = ZI(I,EN)
  301.         ZR(I,EN) = FLOP (SR*YR-SI*YI)
  302.         ZI(I,EN) = FLOP (SR*YI+SI*YR)
  303. 640   CONTINUE
  304.       GO TO 240
  305. C
  306. C ***      A ROOT FOUND
  307. 660   CONTINUE
  308.       HR(EN,EN) = FLOP (HR(EN,EN)+TR)
  309.       WR(EN) = HR(EN,EN)
  310.       HI(EN,EN) = FLOP (HI(EN,EN)+TI)
  311.       WI(EN) = HI(EN,EN)
  312.       EN = ENM1
  313.       GO TO 220
  314. C
  315. C ***      ALL ROOTS FOUND.  BACKSUBSTITUTE TO FIND VECTORS OF UPPER
  316. C            TRIANGULAR FORM
  317. C
  318. 680   CONTINUE
  319.       IF (JOB.NE.2) GO TO 1001
  320.       NORM = 0.0D0
  321.       DO 721 I = 1, N
  322.         DO 720 J = I, N
  323.           TR = FLOP (DABS (HR(I,J)))+FLOP (DABS (HI(I,J)))
  324.           IF (TR.GT.NORM) NORM = TR
  325. 720     CONTINUE
  326. 721   CONTINUE
  327.       IF (N.EQ.1 .OR. NORM.EQ.0.0D0) GO TO 1001
  328. C
  329. C ***      FOR EN = N STEP -1 UNTIL 2 DO
  330.       DO 800 NN = 2, N
  331.         EN = N+2-NN
  332.         XR = WR(EN)
  333.         XI = WI(EN)
  334.         HR(EN,EN) = 1.0D0
  335.         HI(EN,EN) = 0.0D0
  336.         ENM1 = EN-1
  337. C
  338. C ***      FOR I = EN-1 STEP -1 UNTIL 1 DO
  339.         DO 780 II = 1, ENM1
  340.           I = EN-II
  341.           ZZR = 0.0D0
  342.           ZZI = 0.0D0
  343.           IP1 = I+1
  344.           DO 740 J = IP1, EN
  345.             ZZR = FLOP (ZZR+HR(I,J)*HR(J,EN)-HI(I,J)*HI(J,EN))
  346.             ZZI = FLOP (ZZI+HR(I,J)*HI(J,EN)+HI(I,J)*HR(J,EN))
  347. 740       CONTINUE
  348.           YR = FLOP (XR-WR(I))
  349.           YI = FLOP (XI-WI(I))
  350.           IF (YR.NE.0.0D0 .OR. YI.NE.0.0D0) GO TO 765
  351.           YR = NORM
  352. 760       CONTINUE
  353.           YR = FLOP (YR/100.0D0)
  354.           YI = FLOP (NORM+YR)
  355.           IF (YI.NE.NORM) GO TO 760
  356.           YI = 0.0D0
  357. 765       CONTINUE
  358.           CALL WDIV (ZZR, ZZI, YR, YI, HR(I,EN), HI(I,EN))
  359.           TR = FLOP (DABS (HR(I,EN)))+FLOP (DABS (HI(I,EN)))
  360.           IF (TR.EQ.0.0D0) GO TO 780
  361.           IF (TR+1.0D0/TR.GT.TR) GO TO 780
  362.           DO 770 J = I, EN
  363.             HR(J,EN) = FLOP (HR(J,EN)/TR)
  364.             HI(J,EN) = FLOP (HI(J,EN)/TR)
  365. 770       CONTINUE
  366. 780     CONTINUE
  367. 800   CONTINUE
  368. C
  369. C ***      END BACKSUBSTITUTION
  370.       ENM1 = N-1
  371. C
  372. C ***      VECTORS OF ISOLATED ROOTS
  373.       DO  840 I = 1, ENM1
  374.         IF (I.GE.LOW .AND. I.LE.IGH) GO TO 840
  375.         IP1 = I+1
  376.         DO 820 J = IP1, N
  377.           ZR(I,J) = HR(I,J)
  378.           ZI(I,J) = HI(I,J)
  379. 820     CONTINUE
  380. 840   CONTINUE
  381. C
  382. C ***      MULTIPLY BY TRANSFORMATION MATRIX TO GIVE VECTORS OF ORIGINAL
  383. C            FULL MATRIX.
  384. C ***      FOR J = N STEP -1 UNTIL LOW+1 DO
  385.       DO 881 JJ = LOW, ENM1
  386.         J = N+LOW-JJ
  387.         M = MIN0 (J, IGH)
  388.         DO 880 I = LOW, IGH
  389.           ZZR = 0.0D0
  390.           ZZI = 0.0D0
  391.           DO 860 K = LOW, M
  392.             ZZR = FLOP (ZZR+ZR(I,K)*HR(K,J)-ZI(I,K)*HI(K,J))
  393.             ZZI = FLOP (ZZI+ZR(I,K)*HI(K,J)+ZI(I,K)*HR(K,J))
  394. 860       CONTINUE
  395.           ZR(I,J) = ZZR
  396.           ZI(I,J) = ZZI
  397. 880     CONTINUE
  398. 881   CONTINUE
  399.       GO TO 1001
  400. C
  401. C ***      SET ERROR -- NO CONVERGENCE TO AN EIGENVALUE AFTER 30 ITERATIONS
  402. 1000  CONTINUE
  403.       IERR = EN
  404. C
  405. 1001  CONTINUE
  406.       RETURN
  407.       END
  408.