[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / f / invit.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1996/04/01 15:02:37  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
      SUBROUTINE INVIT(NM,N,A,WR,WI,SELECT,MM,M,Z,IERR,RM3,RV1,RV2)
      INTEGER I,J,K,M,N,S,II,IP,MM,MP,NM,UK,IP1,ITS,KM1,IERR
      REAL A(NM,N),WR(N),WI(N),Z(NM,MM),RM3(N,*),RV1(N),RV2(N)
      REAL W,X,Y,EPS3,NORM,NORMV,GROWTO,ILAMBD,MACHEP,RLAMBD,UKROOT
      LOGICAL SELECT(N)
      COMPLEX Z3
      REAL T3(2)
      EQUIVALENCE (Z3,T3(1))
#if defined(CERNLIB_CDC)
      MACHEP=2.**(-47)
#endif
#if !defined(CERNLIB_CDC)
      MACHEP=2.**(-23)
#endif
      IERR = 0
      UK = 0
      S = 1
      IP = 0
      DO 980 K = 1, N
         IF (WI(K) .EQ. 0.0 .OR. IP .LT. 0) GO TO 100
         IP = 1
         IF (SELECT(K) .AND. SELECT(K+1)) SELECT(K+1) = .FALSE.
  100    IF (.NOT. SELECT(K)) GO TO 960
         IF (WI(K) .NE. 0.0) S = S + 1
         IF (S .GT. MM) GO TO 1000
         IF (UK .GE. K) GO TO 200
         DO 120 UK = K, N
            IF (UK .EQ. N) GO TO 140
            IF (A(UK+1,UK) .EQ. 0.0) GO TO 140
  120    CONTINUE
  140    NORM = 0.0
         MP = 1
         DO 180 I = 1, UK
            X = 0.0
            DO 160 J = MP, UK
  160       X = X + ABS(A(I,J))
            IF (X .GT. NORM) NORM = X
            MP = I
  180    CONTINUE
         IF (NORM .EQ. 0.0) NORM = 1.0
         EPS3 = MACHEP * NORM
         UKROOT = SQRT(REAL(UK))
         GROWTO = 1.0E-1 / UKROOT
  200    RLAMBD = WR(K)
         ILAMBD = WI(K)
         IF (K .EQ. 1) GO TO 280
         KM1 = K - 1
         GO TO 240
  220    RLAMBD = RLAMBD + EPS3
  240    DO 260 II = 1, KM1
            I = K - II
            IF (SELECT(I) .AND. ABS(WR(I)-RLAMBD) .LT. EPS3 .AND.
     X         ABS(WI(I)-ILAMBD) .LT. EPS3) GO TO 220
  260    CONTINUE
         WR(K) = RLAMBD
         IP1 = K + IP
         WR(IP1) = RLAMBD
  280    MP = 1
         DO 320 I = 1, UK
            DO 300 J = MP, UK
  300       RM3(J,I) = A(I,J)
            RM3(I,I) = RM3(I,I) - RLAMBD
            MP = I
            RV1(I) = EPS3
  320    CONTINUE
         ITS = 0
         IF (ILAMBD .NE. 0.0) GO TO 520
         IF (UK .EQ. 1) GO TO 420
         DO 400 I = 2, UK
            MP = I - 1
            IF (ABS(RM3(MP,I)) .LE. ABS(RM3(MP,MP))) GO TO 360
            DO 340 J = MP, UK
               Y = RM3(J,I)
               RM3(J,I) = RM3(J,MP)
               RM3(J,MP) = Y
  340       CONTINUE
  360       IF (RM3(MP,MP) .EQ. 0.0) RM3(MP,MP) = EPS3
            X = RM3(MP,I) / RM3(MP,MP)
            IF (X .EQ. 0.0) GO TO 400
            DO 380 J = I, UK
  380       RM3(J,I) = RM3(J,I) - X * RM3(J,MP)
  400    CONTINUE
  420    IF (RM3(UK,UK) .EQ. 0.0) RM3(UK,UK) = EPS3
  440    DO 500 II = 1, UK
            I = UK + 1 - II
            Y = RV1(I)
            IF (I .EQ. UK) GO TO 480
            IP1 = I + 1
            DO 460 J = IP1, UK
  460       Y = Y - RM3(J,I) * RV1(J)
  480       RV1(I) = Y / RM3(I,I)
  500    CONTINUE
         GO TO 740
  520    RM3(1,3) = -ILAMBD
         DO 540 I = 2, UK
  540    RM3(1,I+2) = 0.0
         DO 640 I = 2, UK
            MP = I - 1
            W = RM3(MP,I)
            X = RM3(MP,MP) * RM3(MP,MP) + RM3(MP,I+1) * RM3(MP,I+1)
            IF (W * W .LE. X) GO TO 580
            X = RM3(MP,MP) / W
            Y = RM3(MP,I+1) / W
            RM3(MP,MP) = W
            RM3(MP,I+1) = 0.0
            DO 560 J = I, UK
               W = RM3(J,I)
               RM3(J,I) = RM3(J,MP) - X * W
               RM3(J,MP) = W
               RM3(I,J+2) = RM3(MP,J+2) - Y * W
               RM3(MP,J+2) = 0.0
  560       CONTINUE
            RM3(MP,I+2) = -ILAMBD
            RM3(I,I) = RM3(I,I) - Y * ILAMBD
            RM3(I,I+2) = RM3(I,I+2) + X * ILAMBD
            GO TO 640
  580       IF (X .NE. 0.0) GO TO 600
            RM3(MP,MP) = EPS3
            RM3(MP,I+1) = 0.0
            X = EPS3 * EPS3
  600       W = W / X
            X = RM3(MP,MP) * W
            Y = -RM3(MP,I+1) * W
            DO 620 J = I, UK
               RM3(J,I) = RM3(J,I) - X * RM3(J,MP) + Y * RM3(MP,J+2)
               RM3(I,J+2) = -X * RM3(MP,J+2) - Y * RM3(J,MP)
  620       CONTINUE
            RM3(I,I+2) = RM3(I,I+2) - ILAMBD
  640    CONTINUE
         IF (RM3(UK,UK) .EQ. 0.0 .AND.
     X      RM3(UK,UK+2) .EQ. 0.0) RM3(UK,UK) = EPS3
  660    DO 720 II = 1, UK
            I = UK + 1 - II
            X = RV1(I)
            Y = 0.0
            IF (I .EQ. UK) GO TO 700
            IP1 = I + 1
            DO 680 J = IP1, UK
               X = X - RM3(J,I) * RV1(J) + RM3(I,J+2) * RV2(J)
               Y = Y - RM3(J,I) * RV2(J) - RM3(I,J+2) * RV1(J)
  680       CONTINUE
  700       Z3 = CMPLX(X,Y) / CMPLX(RM3(I,I),RM3(I,I+2))
            RV1(I) = T3(1)
            RV2(I) = T3(2)
  720    CONTINUE
  740    ITS = ITS + 1
         NORM = 0.0
         NORMV = 0.0
         DO 780 I = 1, UK
            IF (ILAMBD .EQ. 0.0) X = ABS(RV1(I))
            IF (ILAMBD .NE. 0.0) X = ABS(CMPLX(RV1(I),RV2(I)))
            IF (NORMV .GE. X) GO TO 760
            NORMV = X
            J = I
  760       NORM = NORM + X
  780    CONTINUE
         IF (NORM .LT. GROWTO) GO TO 840
         X = RV1(J)
         IF (ILAMBD .EQ. 0.0) X = 1.0 / X
         IF (ILAMBD .NE. 0.0) Y = RV2(J)
         DO 820 I = 1, UK
            IF (ILAMBD .NE. 0.0) GO TO 800
            Z(I,S) = RV1(I) * X
            GO TO 820
  800       Z3 = CMPLX(RV1(I),RV2(I)) / CMPLX(X,Y)
            Z(I,S-1) = T3(1)
            Z(I,S) = T3(2)
  820    CONTINUE
         IF (UK .EQ. N) GO TO 940
         J = UK + 1
         GO TO 900
  840    IF (ITS .GE. UK) GO TO 880
         X = UKROOT
         Y = EPS3 / (X + 1.0)
         RV1(1) = EPS3
         DO 860 I = 2, UK
  860    RV1(I) = Y
         J = UK - ITS + 1
         RV1(J) = RV1(J) - EPS3 * X
         IF (ILAMBD .EQ. 0.0) GO TO 440
         GO TO 660
  880    J = 1
         IERR = -K
  900    DO 920 I = J, N
            Z(I,S) = 0.0
            IF (ILAMBD .NE. 0.0) Z(I,S-1) = 0.0
  920    CONTINUE
  940    S = S + 1
  960    IF (IP .EQ. (-1)) IP = 0
         IF (IP .EQ. 1) IP = -1
  980 CONTINUE
      GO TO 1001
 1000 IF (IERR .NE. 0) IERR = IERR - N
      IF (IERR .EQ. 0) IERR = -(2 * N + 1)
 1001 M = S - 1 - ABS(IP)
      RETURN
      END
Commit	Line	Data
fe4da5cc	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