MINICERN/mathlib/gen/f/invit.F

   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