1 *
2 * \$Id\$
3 *
4 * \$Log\$
5 * Revision 1.1.1.1  1996/04/01 15:02:37  mclareni
6 * Mathlib gen
7 *
8 *
9 #include "gen/pilot.h"
10       SUBROUTINE INVIT(NM,N,A,WR,WI,SELECT,MM,M,Z,IERR,RM3,RV1,RV2)
11       INTEGER I,J,K,M,N,S,II,IP,MM,MP,NM,UK,IP1,ITS,KM1,IERR
12       REAL A(NM,N),WR(N),WI(N),Z(NM,MM),RM3(N,*),RV1(N),RV2(N)
13       REAL W,X,Y,EPS3,NORM,NORMV,GROWTO,ILAMBD,MACHEP,RLAMBD,UKROOT
14       LOGICAL SELECT(N)
15       COMPLEX Z3
16       REAL T3(2)
17       EQUIVALENCE (Z3,T3(1))
18 #if defined(CERNLIB_CDC)
19       MACHEP=2.**(-47)
20 #endif
21 #if !defined(CERNLIB_CDC)
22       MACHEP=2.**(-23)
23 #endif
24       IERR = 0
25       UK = 0
26       S = 1
27       IP = 0
28       DO 980 K = 1, N
29          IF (WI(K) .EQ. 0.0 .OR. IP .LT. 0) GO TO 100
30          IP = 1
31          IF (SELECT(K) .AND. SELECT(K+1)) SELECT(K+1) = .FALSE.
32   100    IF (.NOT. SELECT(K)) GO TO 960
33          IF (WI(K) .NE. 0.0) S = S + 1
34          IF (S .GT. MM) GO TO 1000
35          IF (UK .GE. K) GO TO 200
36          DO 120 UK = K, N
37             IF (UK .EQ. N) GO TO 140
38             IF (A(UK+1,UK) .EQ. 0.0) GO TO 140
39   120    CONTINUE
40   140    NORM = 0.0
41          MP = 1
42          DO 180 I = 1, UK
43             X = 0.0
44             DO 160 J = MP, UK
45   160       X = X + ABS(A(I,J))
46             IF (X .GT. NORM) NORM = X
47             MP = I
48   180    CONTINUE
49          IF (NORM .EQ. 0.0) NORM = 1.0
50          EPS3 = MACHEP * NORM
51          UKROOT = SQRT(REAL(UK))
52          GROWTO = 1.0E-1 / UKROOT
53   200    RLAMBD = WR(K)
54          ILAMBD = WI(K)
55          IF (K .EQ. 1) GO TO 280
56          KM1 = K - 1
57          GO TO 240
58   220    RLAMBD = RLAMBD + EPS3
59   240    DO 260 II = 1, KM1
60             I = K - II
61             IF (SELECT(I) .AND. ABS(WR(I)-RLAMBD) .LT. EPS3 .AND.
62      X         ABS(WI(I)-ILAMBD) .LT. EPS3) GO TO 220
63   260    CONTINUE
64          WR(K) = RLAMBD
65          IP1 = K + IP
66          WR(IP1) = RLAMBD
67   280    MP = 1
68          DO 320 I = 1, UK
69             DO 300 J = MP, UK
70   300       RM3(J,I) = A(I,J)
71             RM3(I,I) = RM3(I,I) - RLAMBD
72             MP = I
73             RV1(I) = EPS3
74   320    CONTINUE
75          ITS = 0
76          IF (ILAMBD .NE. 0.0) GO TO 520
77          IF (UK .EQ. 1) GO TO 420
78          DO 400 I = 2, UK
79             MP = I - 1
80             IF (ABS(RM3(MP,I)) .LE. ABS(RM3(MP,MP))) GO TO 360
81             DO 340 J = MP, UK
82                Y = RM3(J,I)
83                RM3(J,I) = RM3(J,MP)
84                RM3(J,MP) = Y
85   340       CONTINUE
86   360       IF (RM3(MP,MP) .EQ. 0.0) RM3(MP,MP) = EPS3
87             X = RM3(MP,I) / RM3(MP,MP)
88             IF (X .EQ. 0.0) GO TO 400
89             DO 380 J = I, UK
90   380       RM3(J,I) = RM3(J,I) - X * RM3(J,MP)
91   400    CONTINUE
92   420    IF (RM3(UK,UK) .EQ. 0.0) RM3(UK,UK) = EPS3
93   440    DO 500 II = 1, UK
94             I = UK + 1 - II
95             Y = RV1(I)
96             IF (I .EQ. UK) GO TO 480
97             IP1 = I + 1
98             DO 460 J = IP1, UK
99   460       Y = Y - RM3(J,I) * RV1(J)
100   480       RV1(I) = Y / RM3(I,I)
101   500    CONTINUE
102          GO TO 740
103   520    RM3(1,3) = -ILAMBD
104          DO 540 I = 2, UK
105   540    RM3(1,I+2) = 0.0
106          DO 640 I = 2, UK
107             MP = I - 1
108             W = RM3(MP,I)
109             X = RM3(MP,MP) * RM3(MP,MP) + RM3(MP,I+1) * RM3(MP,I+1)
110             IF (W * W .LE. X) GO TO 580
111             X = RM3(MP,MP) / W
112             Y = RM3(MP,I+1) / W
113             RM3(MP,MP) = W
114             RM3(MP,I+1) = 0.0
115             DO 560 J = I, UK
116                W = RM3(J,I)
117                RM3(J,I) = RM3(J,MP) - X * W
118                RM3(J,MP) = W
119                RM3(I,J+2) = RM3(MP,J+2) - Y * W
120                RM3(MP,J+2) = 0.0
121   560       CONTINUE
122             RM3(MP,I+2) = -ILAMBD
123             RM3(I,I) = RM3(I,I) - Y * ILAMBD
124             RM3(I,I+2) = RM3(I,I+2) + X * ILAMBD
125             GO TO 640
126   580       IF (X .NE. 0.0) GO TO 600
127             RM3(MP,MP) = EPS3
128             RM3(MP,I+1) = 0.0
129             X = EPS3 * EPS3
130   600       W = W / X
131             X = RM3(MP,MP) * W
132             Y = -RM3(MP,I+1) * W
133             DO 620 J = I, UK
134                RM3(J,I) = RM3(J,I) - X * RM3(J,MP) + Y * RM3(MP,J+2)
135                RM3(I,J+2) = -X * RM3(MP,J+2) - Y * RM3(J,MP)
136   620       CONTINUE
137             RM3(I,I+2) = RM3(I,I+2) - ILAMBD
138   640    CONTINUE
139          IF (RM3(UK,UK) .EQ. 0.0 .AND.
140      X      RM3(UK,UK+2) .EQ. 0.0) RM3(UK,UK) = EPS3
141   660    DO 720 II = 1, UK
142             I = UK + 1 - II
143             X = RV1(I)
144             Y = 0.0
145             IF (I .EQ. UK) GO TO 700
146             IP1 = I + 1
147             DO 680 J = IP1, UK
148                X = X - RM3(J,I) * RV1(J) + RM3(I,J+2) * RV2(J)
149                Y = Y - RM3(J,I) * RV2(J) - RM3(I,J+2) * RV1(J)
150   680       CONTINUE
151   700       Z3 = CMPLX(X,Y) / CMPLX(RM3(I,I),RM3(I,I+2))
152             RV1(I) = T3(1)
153             RV2(I) = T3(2)
154   720    CONTINUE
155   740    ITS = ITS + 1
156          NORM = 0.0
157          NORMV = 0.0
158          DO 780 I = 1, UK
159             IF (ILAMBD .EQ. 0.0) X = ABS(RV1(I))
160             IF (ILAMBD .NE. 0.0) X = ABS(CMPLX(RV1(I),RV2(I)))
161             IF (NORMV .GE. X) GO TO 760
162             NORMV = X
163             J = I
164   760       NORM = NORM + X
165   780    CONTINUE
166          IF (NORM .LT. GROWTO) GO TO 840
167          X = RV1(J)
168          IF (ILAMBD .EQ. 0.0) X = 1.0 / X
169          IF (ILAMBD .NE. 0.0) Y = RV2(J)
170          DO 820 I = 1, UK
171             IF (ILAMBD .NE. 0.0) GO TO 800
172             Z(I,S) = RV1(I) * X
173             GO TO 820
174   800       Z3 = CMPLX(RV1(I),RV2(I)) / CMPLX(X,Y)
175             Z(I,S-1) = T3(1)
176             Z(I,S) = T3(2)
177   820    CONTINUE
178          IF (UK .EQ. N) GO TO 940
179          J = UK + 1
180          GO TO 900
181   840    IF (ITS .GE. UK) GO TO 880
182          X = UKROOT
183          Y = EPS3 / (X + 1.0)
184          RV1(1) = EPS3
185          DO 860 I = 2, UK
186   860    RV1(I) = Y
187          J = UK - ITS + 1
188          RV1(J) = RV1(J) - EPS3 * X
189          IF (ILAMBD .EQ. 0.0) GO TO 440
190          GO TO 660
191   880    J = 1
192          IERR = -K
193   900    DO 920 I = J, N
194             Z(I,S) = 0.0
195             IF (ILAMBD .NE. 0.0) Z(I,S-1) = 0.0
196   920    CONTINUE
197   940    S = S + 1
198   960    IF (IP .EQ. (-1)) IP = 0
199          IF (IP .EQ. 1) IP = -1
200   980 CONTINUE
201       GO TO 1001
202  1000 IF (IERR .NE. 0) IERR = IERR - N
203       IF (IERR .EQ. 0) IERR = -(2 * N + 1)
204  1001 M = S - 1 - ABS(IP)
205       RETURN
206       END