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 / CORTH.FOR < prev    next >
Encoding:
Text File  |  1991-04-13  |  3.9 KB  |  129 lines
  1.       SUBROUTINE CORTH (NM, N, LOW, IGH, AR, AI, ORTR, ORTI)
  2.       IMPLICIT NONE
  3. C
  4. C A TRANSLATION OF A COMPLEX ANALOGUE OF THE ALGOL PROCEDURE ORTHES,
  5. C   NUM. MATH. 12, 349-368 (1968) BY MARTIN AND WILKINSON.
  6. C   HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358 (1971).
  7. C
  8. C GIVEN A COMPLEX GENERAL MATRIX, IT REDUCES A SUBMATRIX SITUATED IN ROWS
  9. C   AND COLUMNS LOW THROUGH IGH TO UPPER HESSENBERG FORM BY UNITARY
  10. C   SIMILARITY TRANSFORMATIONS.
  11. C
  12. C ON INPUT:
  13. C   NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
  14. C      AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT.
  15. C
  16. C   N  IS THE ORDER OF THE MATRIX.
  17. C
  18. C   LOW & IGH ARE INTEGERS DETERMINED BY THE BALANCING SUBROUTINE CBAL.
  19. C      IF CBAL HAS NOT BEEN USED, SET LOW = 1, IGH = N.
  20. C
  21. C   AR & AI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY, OF THE
  22. C      COMPLEX INPUT MATRIX.
  23. C
  24. C ON OUTPUT:
  25. C   AR & AI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY, OF THE
  26. C      HESSENBERG MATRIX.  INFORMATION ABOUT THE UNITARY TRANSFORMATIONS
  27. C      USED IN THE REDUCTION IS STORED IN THE REMAINING TRIANGLES UNDER
  28. C      THE HESSENBERG MATRIX.
  29. C
  30. C   ORTR & ORTI CONTAIN FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
  31. C      ONLY ELEMENTS LOW THROUGH IGH ARE USED.
  32. C
  33. C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
  34. C   APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  35. C
  36.       INTEGER NM, N, LOW, IGH
  37.       DOUBLE PRECISION AR(NM,N), AI(NM,N), ORTR(IGH), ORTI(IGH)
  38. C
  39.       INTEGER I, J, M, II, JJ, LA, MP, KP1
  40.       DOUBLE PRECISION F, G, H, FI, FR, SCALE
  41. C
  42.       DOUBLE PRECISION FLOP, PYTHAG
  43. C
  44. C
  45.       LA = IGH-1
  46.       KP1 = LOW+1
  47.       IF (LA.LT.KP1) GO TO 200
  48. C
  49.       DO 180 M = KP1, LA
  50.         H = 0.0D0
  51.         ORTR(M) = 0.0D0
  52.         ORTI(M) = 0.0D0
  53.         SCALE = 0.0D0
  54. C
  55. C ***      SCALE COLUMN (ALGOL TOL THEN NOT NEEDED)
  56.         DO 90 I = M, IGH
  57.           SCALE = FLOP (SCALE+DABS (AR(I,M-1))+DABS (AI(I,M-1)))
  58. 90      CONTINUE
  59.         IF (SCALE.EQ.0.0D0) GO TO 180
  60.         MP = M+IGH
  61. C
  62. C ***      FOR I = IGH STEP -1 UNTIL M DO
  63.         DO 100 II = M, IGH
  64.           I = MP-II
  65.           ORTR(I) = FLOP (AR(I,M-1)/SCALE)
  66.           ORTI(I) = FLOP (AI(I,M-1)/SCALE)
  67.           H = FLOP (H+ORTR(I)*ORTR(I)+ORTI(I)*ORTI(I))
  68. 100     CONTINUE
  69.         G = FLOP (DSQRT (H))
  70.         F = PYTHAG (ORTR(M), ORTI(M))
  71.         IF (F.EQ.0.0D0) GO TO 103
  72.         H = FLOP (H+F*G)
  73.         G = FLOP (G/F)
  74.         ORTR(M) = FLOP ((1.0D0+G)*ORTR(M))
  75.         ORTI(M) = FLOP ((1.0D0+G)*ORTI(M))
  76.         GO TO 105
  77. C
  78. 103     CONTINUE
  79.         ORTR(M) = G
  80.         AR(M,M-1) = SCALE
  81. C
  82. C ***      FORM (I-(U*UT)/H)*A
  83. 105     CONTINUE
  84.         DO 130 J = M, N
  85.           FR = 0.0D0
  86.           FI = 0.0D0
  87. C
  88. C ***      FOR I = IGH STEP -1 UNTIL M DO
  89.           DO 110 II = M, IGH
  90.             I = MP-II
  91.             FR = FLOP (FR+ORTR(I)*AR(I,J)+ORTI(I)*AI(I,J))
  92.             FI = FLOP (FI+ORTR(I)*AI(I,J)-ORTI(I)*AR(I,J))
  93. 110       CONTINUE
  94.           FR = FLOP (FR/H)
  95.           FI = FLOP (FI/H)
  96.           DO 120 I = M, IGH
  97.             AR(I,J) = FLOP (AR(I,J)-FR*ORTR(I)+FI*ORTI(I))
  98.             AI(I,J) = FLOP (AI(I,J)-FR*ORTI(I)-FI*ORTR(I))
  99. 120       CONTINUE
  100. 130     CONTINUE
  101. C
  102. C ***      FORM (I-(U*UT)/H)*A*(I-(U*UT)/H)
  103.         DO 160 I = 1, IGH
  104.           FR = 0.0D0
  105.           FI = 0.0D0
  106. C
  107. C ***      FOR J = IGH STEP -1 UNTIL M DO
  108.           DO 140 JJ = M, IGH
  109.             J = MP-JJ
  110.             FR = FLOP (FR+ORTR(J)*AR(I,J)-ORTI(J)*AI(I,J))
  111.             FI = FLOP (FI+ORTR(J)*AI(I,J)+ORTI(J)*AR(I,J))
  112. 140       CONTINUE
  113.           FR = FLOP (FR/H)
  114.           FI = FLOP (FI/H)
  115.           DO 150 J = M, IGH
  116.             AR(I,J) = FLOP (AR(I,J)-FR*ORTR(J)-FI*ORTI(J))
  117.             AI(I,J) = FLOP (AI(I,J)+FR*ORTI(J)-FI*ORTR(J))
  118. 150       CONTINUE
  119. 160     CONTINUE
  120.         ORTR(M) = FLOP (SCALE*ORTR(M))
  121.         ORTI(M) = FLOP (SCALE*ORTI(M))
  122.         AR(M,M-1) = FLOP (-G*AR(M,M-1))
  123.         AI(M,M-1) = FLOP (-G*AI(M,M-1))
  124. 180   CONTINUE
  125. C
  126. 200   CONTINUE
  127.       RETURN
  128.       END
  129.