MINICERN/mathlib/gen/f/tsturm.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:38  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       SUBROUTINE TSTURM(NM,N,EPS1,D,E,E2,LB,UB,MM,M,W,Z,
  11      X                  IERR,RV1,RV2,RV3,RV4,RV5,RV6)
  12       INTEGER I,J,K,M,N,P,Q,R,S,II,IP,JJ,MM,M1,M2,NM,ITS,
  13      X        IERR,GROUP,ISTURM
  14       REAL D(N),E(N),E2(N),W(MM),Z(NM,MM),
  15      X       RV1(N),RV2(N),RV3(N),RV4(N),RV5(N),RV6(N)
  16       REAL U,V,LB,T1,T2,UB,UK,XU,X0,X1,EPS1,EPS2,EPS3,EPS4,
  17      X       NORM,MACHEP
  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       T1 = LB
  26       T2 = UB
  27       DO 40 I = 1, N
  28          IF (I .EQ. 1) GO TO 20
  29          IF (ABS(E(I)) .GT. MACHEP * (ABS(D(I)) + ABS(D(I-1))))
  30      X      GO TO 40
  31    20    E2(I) = 0.0
  32    40 CONTINUE
  33       P = 1
  34       Q = N
  35       X1 = UB
  36       ISTURM = 1
  37       GO TO 320
  38    60 M = S
  39       X1 = LB
  40       ISTURM = 2
  41       GO TO 320
  42    80 M = M - S
  43       IF (M .GT. MM) GO TO 980
  44       Q = 0
  45       R = 0
  46   100 IF (R .EQ. M) GO TO 1001
  47       P = Q + 1
  48       XU = D(P)
  49       X0 = D(P)
  50       U = 0.0
  51       DO 120 Q = P, N
  52          X1 = U
  53          U = 0.0
  54          V = 0.0
  55          IF (Q .EQ. N) GO TO 110
  56          U = ABS(E(Q+1))
  57          V = E2(Q+1)
  58   110    XU = MIN(D(Q)-(X1+U),XU)
  59          X0 = MAX(D(Q)+(X1+U),X0)
  60          IF (V .EQ. 0.0) GO TO 140
  61   120 CONTINUE
  62   140 X1 = MAX(ABS(XU),ABS(X0)) * MACHEP
  63       IF (EPS1 .LE. 0.0) EPS1 = -X1
  64       IF (P .NE. Q) GO TO 180
  65       IF (T1 .GT. D(P) .OR. D(P) .GE. T2) GO TO 940
  66       R = R + 1
  67       DO 160 I = 1, N
  68   160 Z(I,R) = 0.0
  69       W(R) = D(P)
  70       Z(P,R) = 1.0
  71       GO TO 940
  72   180 X1 = X1 * (Q-P+1)
  73       LB = MAX(T1,XU-X1)
  74       UB = MIN(T2,X0+X1)
  75       X1 = LB
  76       ISTURM = 3
  77       GO TO 320
  78   200 M1 = S + 1
  79       X1 = UB
  80       ISTURM = 4
  81       GO TO 320
  82   220 M2 = S
  83       IF (M1 .GT. M2) GO TO 940
  84       X0 = UB
  85       ISTURM = 5
  86       DO 240 I = M1, M2
  87          RV5(I) = UB
  88          RV4(I) = LB
  89   240 CONTINUE
  90       K = M2
  91   250    XU = LB
  92          DO 260 II = M1, K
  93             I = M1 + K - II
  94             IF (XU .GE. RV4(I)) GO TO 260
  95             XU = RV4(I)
  96             GO TO 280
  97   260    CONTINUE
  98   280    IF (X0 .GT. RV5(K)) X0 = RV5(K)
  99   300    X1 = (XU + X0) * 0.5
 100          IF ((X0 - XU) .LE. (2.0 * MACHEP *
 101      X      (ABS(XU) + ABS(X0)) + ABS(EPS1))) GO TO 420
 102   320    S = P - 1
 103          U = 1.0
 104          DO 340 I = P, Q
 105             IF (U .NE. 0.0) GO TO 325
 106             V = ABS(E(I)) / MACHEP
 107             GO TO 330
 108   325       V = E2(I) / U
 109   330       U = D(I) - X1 - V
 110             IF (U .LT. 0.0) S = S + 1
 111   340    CONTINUE
 112          GO TO (60,80,200,220,360), ISTURM
 113   360    IF (S .GE. K) GO TO 400
 114          XU = X1
 115          IF (S .GE. M1) GO TO 380
 116          RV4(M1) = X1
 117          GO TO 300
 118   380    RV4(S+1) = X1
 119          IF (RV5(S) .GT. X1) RV5(S) = X1
 120          GO TO 300
 121   400    X0 = X1
 122          GO TO 300
 123   420    RV5(K) = X1
 124       K = K - 1
 125       IF (K .GE. M1) GO TO 250
 126       NORM = ABS(D(P))
 127       IP = P + 1
 128       DO 500 I = IP, Q
 129   500 NORM = NORM + ABS(D(I)) + ABS(E(I))
 130       EPS2 = 1.0E-3 * NORM
 131       EPS3 = MACHEP * NORM
 132       UK = Q-P+1
 133       EPS4 = UK * EPS3
 134       UK = EPS4 / SQRT(UK)
 135       GROUP = 0
 136       S = P
 137       DO 920 K = M1, M2
 138          R = R + 1
 139          ITS = 1
 140          W(R) = RV5(K)
 141          X1 = RV5(K)
 142          IF (K .EQ. M1) GO TO 520
 143          IF (X1 - X0 .GE. EPS2) GROUP = -1
 144          GROUP = GROUP + 1
 145          IF (X1 .LE. X0) X1 = X0 + EPS3
 146   520    V = 0.0
 147          DO 580 I = P, Q
 148             RV6(I) = UK
 149             IF (I .EQ. P) GO TO 560
 150             IF (ABS(E(I)) .LT. ABS(U)) GO TO 540
 151             XU = U / E(I)
 152             RV4(I) = XU
 153             RV1(I-1) = E(I)
 154             RV2(I-1) = D(I) - X1
 155             RV3(I-1) = 0.0
 156             IF (I .NE. Q) RV3(I-1) = E(I+1)
 157             U = V - XU * RV2(I-1)
 158             V = -XU * RV3(I-1)
 159             GO TO 580
 160   540       XU = E(I) / U
 161             RV4(I) = XU
 162             RV1(I-1) = U
 163             RV2(I-1) = V
 164             RV3(I-1) = 0.0
 165   560       U = D(I) - X1 - XU * V
 166             IF (I .NE. Q) V = E(I+1)
 167   580    CONTINUE
 168          IF (U .EQ. 0.0) U = EPS3
 169          RV1(Q) = U
 170          RV2(Q) = 0.0
 171          RV3(Q) = 0.0
 172   600    DO 620 II = P, Q
 173             I = P + Q - II
 174             RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
 175             V = U
 176             U = RV6(I)
 177   620    CONTINUE
 178          IF (GROUP .EQ. 0) GO TO 700
 179          DO 680 JJ = 1, GROUP
 180             J = R - GROUP - 1 + JJ
 181             XU = 0.0
 182             DO 640 I = P, Q
 183   640       XU = XU + RV6(I) * Z(I,J)
 184             DO 660 I = P, Q
 185   660       RV6(I) = RV6(I) - XU * Z(I,J)
 186   680    CONTINUE
 187   700    NORM = 0.0
 188          DO 720 I = P, Q
 189   720    NORM = NORM + ABS(RV6(I))
 190          IF (NORM .GE. 1.0) GO TO 840
 191          IF (ITS .EQ. 5) GO TO 960
 192          IF (NORM .NE. 0.0) GO TO 740
 193          RV6(S) = EPS4
 194          S = S + 1
 195          IF (S .GT. Q) S = P
 196          GO TO 780
 197   740    XU = EPS4 / NORM
 198          DO 760 I = P, Q
 199   760    RV6(I) = RV6(I) * XU
 200   780    DO 820 I = IP, Q
 201             U = RV6(I)
 202             IF (RV1(I-1) .NE. E(I)) GO TO 800
 203             U = RV6(I-1)
 204             RV6(I-1) = RV6(I)
 205   800       RV6(I) = U - RV4(I) * RV6(I-1)
 206   820    CONTINUE
 207          ITS = ITS + 1
 208          GO TO 600
 209   840    U = 0.0
 210          DO 860 I = P, Q
 211   860    U = U + RV6(I)**2
 212          XU = 1.0 / SQRT(U)
 213          DO 880 I = 1, N
 214   880    Z(I,R) = 0.0
 215          DO 900 I = P, Q
 216   900    Z(I,R) = RV6(I) * XU
 217          X0 = X1
 218   920 CONTINUE
 219   940 IF (Q .LT. N) GO TO 100
 220       GO TO 1001
 221   960 IERR = 4 * N + R
 222       GO TO 1001
 223   980 IERR = 3 * N + 1
 224  1001 LB = T1
 225       UB = T2
 226       RETURN
 227       END