MINICERN/mathlib/gen/f/tinvit.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 TINVIT(NM,N,D,E,E2,M,W,IND,Z,
  11      X                  IERR,RV1,RV2,RV3,RV4,RV6)
  12       INTEGER I,J,M,N,P,Q,R,S,II,IP,JJ,NM,ITS,TAG,IERR,GROUP
  13       REAL D(N),E(N),E2(N),W(M),Z(NM,M),
  14      X       RV1(N),RV2(N),RV3(N),RV4(N),RV6(N)
  15       REAL U,V,UK,XU,X0,X1,EPS2,EPS3,EPS4,NORM,ORDER,MACHEP
  16       INTEGER IND(M)
  17 #if defined(CERNLIB_CDC)
  18       MACHEP=2.**(-47)
  19 #endif
  20 #if !defined(CERNLIB_CDC)
  21       MACHEP=2.**(-23)
  22 #endif
  23       IERR = 0
  24       TAG = 0
  25       ORDER = 1.0 - E2(1)
  26       Q = 0
  27   100 P = Q + 1
  28       DO 120 Q = P, N
  29          IF (Q .EQ. N) GO TO 140
  30          IF (E2(Q+1) .EQ. 0.0) GO TO 140
  31   120 CONTINUE
  32   140 TAG = TAG + 1
  33       S = 0
  34       DO 920 R = 1, M
  35          IF (IND(R) .NE. TAG) GO TO 920
  36          ITS = 1
  37          X1 = W(R)
  38          IF (S .NE. 0) GO TO 510
  39          XU = 1.0
  40          IF (P .NE. Q) GO TO 490
  41          RV6(P) = 1.0
  42          GO TO 870
  43   490    NORM = ABS(D(P))
  44          IP = P + 1
  45          DO 500 I = IP, Q
  46   500    NORM = NORM + ABS(D(I)) + ABS(E(I))
  47          EPS2 = 1.0E-3 * NORM
  48          EPS3 = MACHEP * NORM
  49          UK = Q-P+1
  50          EPS4 = UK * EPS3
  51          UK = EPS4 / SQRT(UK)
  52          S = P
  53   505    GROUP = 0
  54          GO TO 520
  55   510    IF (ABS(X1-X0) .GE. EPS2) GO TO 505
  56          GROUP = GROUP + 1
  57          IF (ORDER * (X1 - X0) .LE. 0.0) X1 = X0 + ORDER * EPS3
  58   520    V = 0.0
  59          DO 580 I = P, Q
  60             RV6(I) = UK
  61             IF (I .EQ. P) GO TO 560
  62             IF (ABS(E(I)) .LT. ABS(U)) GO TO 540
  63             XU = U / E(I)
  64             RV4(I) = XU
  65             RV1(I-1) = E(I)
  66             RV2(I-1) = D(I) - X1
  67             RV3(I-1) = 0.0
  68             IF (I .NE. Q) RV3(I-1) = E(I+1)
  69             U = V - XU * RV2(I-1)
  70             V = -XU * RV3(I-1)
  71             GO TO 580
  72   540       XU = E(I) / U
  73             RV4(I) = XU
  74             RV1(I-1) = U
  75             RV2(I-1) = V
  76             RV3(I-1) = 0.0
  77   560       U = D(I) - X1 - XU * V
  78             IF (I .NE. Q) V = E(I+1)
  79   580    CONTINUE
  80          IF (U .EQ. 0.0) U = EPS3
  81          RV1(Q) = U
  82          RV2(Q) = 0.0
  83          RV3(Q) = 0.0
  84   600    DO 620 II = P, Q
  85             I = P + Q - II
  86             RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
  87             V = U
  88             U = RV6(I)
  89   620    CONTINUE
  90          IF (GROUP .EQ. 0) GO TO 700
  91          J = R
  92          DO 680 JJ = 1, GROUP
  93   630       J = J - 1
  94             IF (IND(J) .NE. TAG) GO TO 630
  95             XU = 0.0
  96             DO 640 I = P, Q
  97   640       XU = XU + RV6(I) * Z(I,J)
  98             DO 660 I = P, Q
  99   660       RV6(I) = RV6(I) - XU * Z(I,J)
 100   680    CONTINUE
 101   700    NORM = 0.0
 102          DO 720 I = P, Q
 103   720    NORM = NORM + ABS(RV6(I))
 104          IF (NORM .GE. 1.0) GO TO 840
 105          IF (ITS .EQ. 5) GO TO 830
 106          IF (NORM .NE. 0.0) GO TO 740
 107          RV6(S) = EPS4
 108          S = S + 1
 109          IF (S .GT. Q) S = P
 110          GO TO 780
 111   740    XU = EPS4 / NORM
 112          DO 760 I = P, Q
 113   760    RV6(I) = RV6(I) * XU
 114   780    DO 820 I = IP, Q
 115             U = RV6(I)
 116             IF (RV1(I-1) .NE. E(I)) GO TO 800
 117             U = RV6(I-1)
 118             RV6(I-1) = RV6(I)
 119   800       RV6(I) = U - RV4(I) * RV6(I-1)
 120   820    CONTINUE
 121          ITS = ITS + 1
 122          GO TO 600
 123   830    IERR = -R
 124          XU = 0.0
 125          GO TO 870
 126   840    U = 0.0
 127          DO 860 I = P, Q
 128   860    U = U + RV6(I)**2
 129          XU = 1.0 / SQRT(U)
 130   870    DO 880 I = 1, N
 131   880    Z(I,R) = 0.0
 132          DO 900 I = P, Q
 133   900    Z(I,R) = RV6(I) * XU
 134          X0 = X1
 135   920 CONTINUE
 136       IF (Q .LT. N) GO TO 100
 137       RETURN
 138       END