1      SUBROUTINE CORTH(NM,N,LOW,IGH,AR,AI,ORTR,ORTI)
2C***BEGIN PROLOGUE  CORTH
3C***DATE WRITTEN   760101   (YYMMDD)
4C***REVISION DATE  830518   (YYMMDD)
5C***CATEGORY NO.  D4C1B2
6C***KEYWORDS  EIGENVALUES,EIGENVECTORS,EISPACK
7C***AUTHOR  SMITH, B. T., ET AL.
8C***PURPOSE  Reduces complex general matrix to complex upper Hessenberg
9C            using unitary similarity transformations.
10C***DESCRIPTION
11C
12C     This subroutine is a translation of a complex analogue of
13C     the ALGOL procedure ORTHES, NUM. MATH. 12, 349-368(1968)
14C     by Martin and Wilkinson.
15C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
16C
17C     Given a COMPLEX GENERAL matrix, this subroutine
18C     reduces a submatrix situated in rows and columns
19C     LOW through IGH to upper Hessenberg form by
20C     unitary similarity transformations.
21C
22C     On INPUT
23C
24C        NM must be set to the row dimension of two-dimensional
25C          array parameters as declared in the calling program
26C          dimension statement.
27C
28C        N is the order of the matrix.
29C
30C        LOW and IGH are integers determined by the balancing
31C          subroutine  CBAL.  If  CBAL  has not been used,
32C          set LOW=1, IGH=N.
33C
34C        AR and AI contain the real and imaginary parts,
35C          respectively, of the complex input matrix.
36C
37C     On OUTPUT
38C
39C        AR and AI contain the real and imaginary parts,
40C          respectively, of the Hessenberg matrix.  Information
41C          about the unitary transformations used in the reduction
42C          is stored in the remaining triangles under the
43C          Hessenberg matrix.
44C
45C        ORTR and ORTI contain further information about the
46C          transformations.  Only elements LOW through IGH are used.
47C
48C     Calls PYTHAG(A,B) for sqrt(A**2 + B**2).
49C
50C     Questions and comments should be directed to B. S. Garbow,
51C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
52C     ------------------------------------------------------------------
53C***REFERENCES  B. T. SMITH, J. M. BOYLE, J. J. DONGARRA, B. S. GARBOW,
54C                 Y. IKEBE, V. C. KLEMA, C. B. MOLER, *MATRIX EIGEN-
55C                 SYSTEM ROUTINES - EISPACK GUIDE*, SPRINGER-VERLAG,
56C                 1976.
57C***ROUTINES CALLED  PYTHAG
58C***END PROLOGUE  CORTH
59C
60      INTEGER I,J,M,N,II,JJ,LA,MP,NM,IGH,KP1,LOW
61      REAL AR(NM,N),AI(NM,N),ORTR(IGH),ORTI(IGH)
62      REAL F,G,H,FI,FR,SCALE
63      REAL PYTHAG
64C
65C***FIRST EXECUTABLE STATEMENT  CORTH
66      LA = IGH - 1
67      KP1 = LOW + 1
68      IF (LA .LT. KP1) GO TO 200
69C
70      DO 180 M = KP1, LA
71         H = 0.0E0
72         ORTR(M) = 0.0E0
73         ORTI(M) = 0.0E0
74         SCALE = 0.0E0
75C     .......... SCALE COLUMN (ALGOL TOL THEN NOT NEEDED) ..........
76         DO 90 I = M, IGH
77   90    SCALE = SCALE + ABS(AR(I,M-1)) + ABS(AI(I,M-1))
78C
79         IF (SCALE .EQ. 0.0E0) GO TO 180
80         MP = M + IGH
81C     .......... FOR I=IGH STEP -1 UNTIL M DO -- ..........
82         DO 100 II = M, IGH
83            I = MP - II
84            ORTR(I) = AR(I,M-1) / SCALE
85            ORTI(I) = AI(I,M-1) / SCALE
86            H = H + ORTR(I) * ORTR(I) + ORTI(I) * ORTI(I)
87  100    CONTINUE
88C
89         G = SQRT(H)
90         F = PYTHAG(ORTR(M),ORTI(M))
91         IF (F .EQ. 0.0E0) GO TO 103
92         H = H + F * G
93         G = G / F
94         ORTR(M) = (1.0E0 + G) * ORTR(M)
95         ORTI(M) = (1.0E0 + G) * ORTI(M)
96         GO TO 105
97C
98  103    ORTR(M) = G
99         AR(M,M-1) = SCALE
100C     .......... FORM (I-(U*UT)/H) * A ..........
101  105    DO 130 J = M, N
102            FR = 0.0E0
103            FI = 0.0E0
104C     .......... FOR I=IGH STEP -1 UNTIL M DO -- ..........
105            DO 110 II = M, IGH
106               I = MP - II
107               FR = FR + ORTR(I) * AR(I,J) + ORTI(I) * AI(I,J)
108               FI = FI + ORTR(I) * AI(I,J) - ORTI(I) * AR(I,J)
109  110       CONTINUE
110C
111            FR = FR / H
112            FI = FI / H
113C
114            DO 120 I = M, IGH
115               AR(I,J) = AR(I,J) - FR * ORTR(I) + FI * ORTI(I)
116               AI(I,J) = AI(I,J) - FR * ORTI(I) - FI * ORTR(I)
117  120       CONTINUE
118C
119  130    CONTINUE
120C     .......... FORM (I-(U*UT)/H)*A*(I-(U*UT)/H) ..........
121         DO 160 I = 1, IGH
122            FR = 0.0E0
123            FI = 0.0E0
124C     .......... FOR J=IGH STEP -1 UNTIL M DO -- ..........
125            DO 140 JJ = M, IGH
126               J = MP - JJ
127               FR = FR + ORTR(J) * AR(I,J) - ORTI(J) * AI(I,J)
128               FI = FI + ORTR(J) * AI(I,J) + ORTI(J) * AR(I,J)
129  140       CONTINUE
130C
131            FR = FR / H
132            FI = FI / H
133C
134            DO 150 J = M, IGH
135               AR(I,J) = AR(I,J) - FR * ORTR(J) - FI * ORTI(J)
136               AI(I,J) = AI(I,J) + FR * ORTI(J) - FI * ORTR(J)
137  150       CONTINUE
138C
139  160    CONTINUE
140C
141         ORTR(M) = SCALE * ORTR(M)
142         ORTI(M) = SCALE * ORTI(M)
143         AR(M,M-1) = -G * AR(M,M-1)
144         AI(M,M-1) = -G * AI(M,M-1)
145  180 CONTINUE
146C
147  200 RETURN
148      END
149