1      SUBROUTINE CINVIT(NM,N,AR,AI,WR,WI,SELECT,MM,M,ZR,ZI,IERR,RM1,
2     1   RM2,RV1,RV2)
3C***BEGIN PROLOGUE  CINVIT
4C***DATE WRITTEN   760101   (YYMMDD)
5C***REVISION DATE  830518   (YYMMDD)
6C***CATEGORY NO.  D4C2B
7C***KEYWORDS  EIGENVALUES,EIGENVECTORS,EISPACK
8C***AUTHOR  SMITH, B. T., ET AL.
9C***PURPOSE  Computes eigenvectors of a complex upper Hessenberg
10C            associated with specified eigenvalues using inverse
11C            iteration.
12C***DESCRIPTION
13C
14C     This subroutine is a translation of the ALGOL procedure CXINVIT
15C     by Peters and Wilkinson.
16C     HANDBOOK FOR AUTO. COMP. VOL.II-LINEAR ALGEBRA, 418-439(1971).
17C
18C     This subroutine finds those eigenvectors of A COMPLEX UPPER
19C     Hessenberg matrix corresponding to specified eigenvalues,
20C     using inverse iteration.
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        AR and AI contain the real and imaginary parts,
31C          respectively, of the Hessenberg matrix.
32C
33C        WR and WI contain the real and imaginary parts, respectively,
34C          of the eigenvalues of the matrix.  The eigenvalues must be
35C          stored in a manner identical to that of subroutine  COMLR,
36C          which recognizes possible splitting of the matrix.
37C
38C        SELECT specifies the eigenvectors to be found.  The
39C          eigenvector corresponding to the J-th eigenvalue is
40C          specified by setting SELECT(J) to .TRUE.
41C
42C        MM should be set to an upper bound for the number of
43C          eigenvectors to be found.
44C
45C     On OUTPUT
46C
47C        AR, AI, WI, and SELECT are unaltered.
48C
49C        WR may have been altered since close eigenvalues are perturbed
50C          slightly in searching for independent eigenvectors.
51C
52C        M is the number of eigenvectors actually found.
53C
54C        ZR and ZI contain the real and imaginary parts, respectively,
55C          of the eigenvectors.  The eigenvectors are normalized
56C          so that the component of largest magnitude is 1.
57C          Any vector which fails the acceptance test is set to zero.
58C
59C        IERR is set to
60C          Zero       for normal return,
61C          -(2*N+1)   if more than MM eigenvectors have been specified,
62C          -K         if the iteration corresponding to the K-th
63C                     value fails,
64C          -(N+K)     if both error situations occur.
65C
66C        RM1, RM2, RV1, and RV2 are temporary storage arrays.
67C
68C     The ALGOL procedure GUESSVEC appears in CINVIT in line.
69C
70C     Calls PYTHAG(A,B) for sqrt(A**2 + B**2).
71C     Calls CDIV for complex division.
72C
73C     Questions and comments should be directed to B. S. Garbow,
74C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
75C     ------------------------------------------------------------------
76C***REFERENCES  B. T. SMITH, J. M. BOYLE, J. J. DONGARRA, B. S. GARBOW,
77C                 Y. IKEBE, V. C. KLEMA, C. B. MOLER, *MATRIX EIGEN-
78C                 SYSTEM ROUTINES - EISPACK GUIDE*, SPRINGER-VERLAG,
79C                 1976.
80C***ROUTINES CALLED  CDIV,PYTHAG
81C***END PROLOGUE  CINVIT
82C
83      INTEGER I,J,K,M,N,S,II,MM,MP,NM,UK,IP1,ITS,KM1,IERR
84      REAL AR(NM,N),AI(NM,N),WR(N),WI(N),ZR(NM,MM),ZI(NM,MM)
85      REAL RM1(N,N),RM2(N,N),RV1(N),RV2(N)
86      REAL X,Y,EPS3,NORM,NORMV,GROWTO,ILAMBD,RLAMBD,UKROOT
87      REAL PYTHAG
88      LOGICAL SELECT(N)
89C
90C***FIRST EXECUTABLE STATEMENT  CINVIT
91      IERR = 0
92      UK = 0
93      S = 1
94C
95      DO 980 K = 1, N
96         IF (.NOT. SELECT(K)) GO TO 980
97         IF (S .GT. MM) GO TO 1000
98         IF (UK .GE. K) GO TO 200
99C     .......... CHECK FOR POSSIBLE SPLITTING ..........
100         DO 120 UK = K, N
101            IF (UK .EQ. N) GO TO 140
102            IF (AR(UK+1,UK) .EQ. 0.0E0 .AND. AI(UK+1,UK) .EQ. 0.0E0)
103     1         GO TO 140
104  120    CONTINUE
105C     .......... COMPUTE INFINITY NORM OF LEADING UK BY UK
106C                (HESSENBERG) MATRIX ..........
107  140    NORM = 0.0E0
108         MP = 1
109C
110         DO 180 I = 1, UK
111            X = 0.0E0
112C
113            DO 160 J = MP, UK
114  160       X = X + PYTHAG(AR(I,J),AI(I,J))
115C
116            IF (X .GT. NORM) NORM = X
117            MP = I
118  180    CONTINUE
119C     .......... EPS3 REPLACES ZERO PIVOT IN DECOMPOSITION
120C                AND CLOSE ROOTS ARE MODIFIED BY EPS3 ..........
121         IF (NORM .EQ. 0.0E0) NORM = 1.0E0
122         EPS3 = NORM
123  190    EPS3 = 0.5E0*EPS3
124         IF (NORM + EPS3 .GT. NORM) GO TO 190
125         EPS3 = 2.0E0*EPS3
126C     .......... GROWTO IS THE CRITERION FOR GROWTH ..........
127         UKROOT = SQRT(FLOAT(UK))
128         GROWTO = 0.1E0 / UKROOT
129  200    RLAMBD = WR(K)
130         ILAMBD = WI(K)
131         IF (K .EQ. 1) GO TO 280
132         KM1 = K - 1
133         GO TO 240
134C     .......... PERTURB EIGENVALUE IF IT IS CLOSE
135C                TO ANY PREVIOUS EIGENVALUE ..........
136  220    RLAMBD = RLAMBD + EPS3
137C     .......... FOR I=K-1 STEP -1 UNTIL 1 DO -- ..........
138  240    DO 260 II = 1, KM1
139            I = K - II
140            IF (SELECT(I) .AND. ABS(WR(I)-RLAMBD) .LT. EPS3 .AND.
141     1         ABS(WI(I)-ILAMBD) .LT. EPS3) GO TO 220
142  260    CONTINUE
143C
144         WR(K) = RLAMBD
145C     .......... FORM UPPER HESSENBERG (AR,AI)-(RLAMBD,ILAMBD)*I
146C                AND INITIAL COMPLEX VECTOR ..........
147  280    MP = 1
148C
149         DO 320 I = 1, UK
150C
151            DO 300 J = MP, UK
152               RM1(I,J) = AR(I,J)
153               RM2(I,J) = AI(I,J)
154  300       CONTINUE
155C
156            RM1(I,I) = RM1(I,I) - RLAMBD
157            RM2(I,I) = RM2(I,I) - ILAMBD
158            MP = I
159            RV1(I) = EPS3
160  320    CONTINUE
161C     .......... TRIANGULAR DECOMPOSITION WITH INTERCHANGES,
162C                REPLACING ZERO PIVOTS BY EPS3 ..........
163         IF (UK .EQ. 1) GO TO 420
164C
165         DO 400 I = 2, UK
166            MP = I - 1
167            IF (PYTHAG(RM1(I,MP),RM2(I,MP)) .LE.
168     1         PYTHAG(RM1(MP,MP),RM2(MP,MP))) GO TO 360
169C
170            DO 340 J = MP, UK
171               Y = RM1(I,J)
172               RM1(I,J) = RM1(MP,J)
173               RM1(MP,J) = Y
174               Y = RM2(I,J)
175               RM2(I,J) = RM2(MP,J)
176               RM2(MP,J) = Y
177  340       CONTINUE
178C
179  360       IF (RM1(MP,MP) .EQ. 0.0E0 .AND. RM2(MP,MP) .EQ. 0.0E0)
180     1         RM1(MP,MP) = EPS3
181            CALL CDIV(RM1(I,MP),RM2(I,MP),RM1(MP,MP),RM2(MP,MP),X,Y)
182            IF (X .EQ. 0.0E0 .AND. Y .EQ. 0.0E0) GO TO 400
183C
184            DO 380 J = I, UK
185               RM1(I,J) = RM1(I,J) - X * RM1(MP,J) + Y * RM2(MP,J)
186               RM2(I,J) = RM2(I,J) - X * RM2(MP,J) - Y * RM1(MP,J)
187  380       CONTINUE
188C
189  400    CONTINUE
190C
191  420    IF (RM1(UK,UK) .EQ. 0.0E0 .AND. RM2(UK,UK) .EQ. 0.0E0)
192     1      RM1(UK,UK) = EPS3
193         ITS = 0
194C     .......... BACK SUBSTITUTION
195C                FOR I=UK STEP -1 UNTIL 1 DO -- ..........
196  660    DO 720 II = 1, UK
197            I = UK + 1 - II
198            X = RV1(I)
199            Y = 0.0E0
200            IF (I .EQ. UK) GO TO 700
201            IP1 = I + 1
202C
203            DO 680 J = IP1, UK
204               X = X - RM1(I,J) * RV1(J) + RM2(I,J) * RV2(J)
205               Y = Y - RM1(I,J) * RV2(J) - RM2(I,J) * RV1(J)
206  680       CONTINUE
207C
208  700       CALL CDIV(X,Y,RM1(I,I),RM2(I,I),RV1(I),RV2(I))
209  720    CONTINUE
210C     .......... ACCEPTANCE TEST FOR EIGENVECTOR
211C                AND NORMALIZATION ..........
212         ITS = ITS + 1
213         NORM = 0.0E0
214         NORMV = 0.0E0
215C
216         DO 780 I = 1, UK
217            X = PYTHAG(RV1(I),RV2(I))
218            IF (NORMV .GE. X) GO TO 760
219            NORMV = X
220            J = I
221  760       NORM = NORM + X
222  780    CONTINUE
223C
224         IF (NORM .LT. GROWTO) GO TO 840
225C     .......... ACCEPT VECTOR ..........
226         X = RV1(J)
227         Y = RV2(J)
228C
229         DO 820 I = 1, UK
230            CALL CDIV(RV1(I),RV2(I),X,Y,ZR(I,S),ZI(I,S))
231  820    CONTINUE
232C
233         IF (UK .EQ. N) GO TO 940
234         J = UK + 1
235         GO TO 900
236C     .......... IN-LINE PROCEDURE FOR CHOOSING
237C                A NEW STARTING VECTOR ..........
238  840    IF (ITS .GE. UK) GO TO 880
239         X = UKROOT
240         Y = EPS3 / (X + 1.0E0)
241         RV1(1) = EPS3
242C
243         DO 860 I = 2, UK
244  860    RV1(I) = Y
245C
246         J = UK - ITS + 1
247         RV1(J) = RV1(J) - EPS3 * X
248         GO TO 660
249C     .......... SET ERROR -- UNACCEPTED EIGENVECTOR ..........
250  880    J = 1
251         IERR = -K
252C     .......... SET REMAINING VECTOR COMPONENTS TO ZERO ..........
253  900    DO 920 I = J, N
254            ZR(I,S) = 0.0E0
255            ZI(I,S) = 0.0E0
256  920    CONTINUE
257C
258  940    S = S + 1
259  980 CONTINUE
260C
261      GO TO 1001
262C     .......... SET ERROR -- UNDERESTIMATE OF EIGENVECTOR
263C                SPACE REQUIRED ..........
264 1000 IF (IERR .NE. 0) IERR = IERR - N
265      IF (IERR .EQ. 0) IERR = -(2 * N + 1)
266 1001 M = S - 1
267      RETURN
268      END
269