MINICERN/mathlib/gen/f/ratqr.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 RATQR(N,EPS1,D,E,E2,M,W,IND,BD,TYPE,IDEF,IERR)
  11       INTEGER I,J,K,M,N,II,JJ,K1,IDEF,IERR,JDEF
  12       REAL D(N),E(N),E2(N),W(N),BD(N)
  13       REAL F,P,Q,R,S,EP,QP,ERR,TOT,EPS1,DELTA,MACHEP
  14       INTEGER IND(M)
  15       LOGICAL TYPE
  16 #if defined(CERNLIB_CDC)
  17       MACHEP=2.**(-47)
  18 #endif
  19 #if !defined(CERNLIB_CDC)
  20       MACHEP=2.**(-23)
  21 #endif
  22       IERR = 0
  23       JDEF = IDEF
  24       DO 20 I = 1, N
  25    20 W(I) = D(I)
  26       IF (TYPE) GO TO 40
  27       J = 1
  28       GO TO 400
  29    40 ERR = 0.0
  30       S = 0.0
  31       TOT = W(1)
  32       Q = 0.0
  33       J = 0
  34       DO 100 I = 1, N
  35          P = Q
  36          IF (I .EQ. 1) GO TO 60
  37          IF (P .GT. MACHEP * (ABS(D(I)) + ABS(D(I-1)))) GO TO 80
  38    60    E2(I) = 0.0
  39          J = J + 1
  40    80    BD(I) = E2(I)
  41          IND(I) = J
  42          Q = 0.0
  43          IF (I .NE. N) Q = ABS(E(I+1))
  44          TOT = MIN(W(I)-P-Q,TOT)
  45   100 CONTINUE
  46       IF (JDEF .EQ. 1 .AND. TOT .LT. 0.0) GO TO 140
  47       DO 110 I = 1, N
  48   110 W(I) = W(I) - TOT
  49       GO TO 160
  50   140 TOT = 0.0
  51   160 DO 360 K = 1, M
  52   180    TOT = TOT + S
  53          DELTA = W(N) - S
  54          I = N
  55          F = ABS(MACHEP*TOT)
  56          IF (EPS1 .LT. F) EPS1 = F
  57          IF (DELTA .GT. EPS1) GO TO 190
  58          IF (DELTA .LT. (-EPS1)) GO TO 1000
  59          GO TO 300
  60   190    IF (K .EQ. N) GO TO 210
  61          K1 = K + 1
  62          DO 200 J = K1, N
  63             IF (BD(J) .LE. (MACHEP*(W(J)+W(J-1))) ** 2) BD(J) = 0.0
  64   200    CONTINUE
  65   210    F = BD(N) / DELTA
  66          QP = DELTA + F
  67          P = 1.0
  68          IF (K .EQ. N) GO TO 260
  69          K1 = N - K
  70          DO 240 II = 1, K1
  71             I = N - II
  72             Q = W(I) - S - F
  73             R = Q / QP
  74             P = P * R + 1.0
  75             EP = F * R
  76             W(I+1) = QP + EP
  77             DELTA = Q - EP
  78             IF (DELTA .GT. EPS1) GO TO 220
  79             IF (DELTA .LT. (-EPS1)) GO TO 1000
  80             GO TO 300
  81   220       F = BD(I) / Q
  82             QP = DELTA + F
  83             BD(I+1) = QP * EP
  84   240    CONTINUE
  85   260    W(K) = QP
  86          S = QP / P
  87          IF (TOT + S .GT. TOT) GO TO 180
  88          IERR = 5 * N + K
  89          S = 0.0
  90          DELTA = QP
  91          DO 280 J = K, N
  92             IF (W(J) .GT. DELTA) GO TO 280
  93             I = J
  94             DELTA = W(J)
  95   280    CONTINUE
  96   300    IF (I .LT. N) BD(I+1) = BD(I) * F / QP
  97          II = IND(I)
  98          IF (I .EQ. K) GO TO 340
  99          K1 = I - K
 100          DO 320 JJ = 1, K1
 101             J = I - JJ
 102             W(J+1) = W(J) - S
 103             BD(J+1) = BD(J)
 104             IND(J+1) = IND(J)
 105   320    CONTINUE
 106   340    W(K) = TOT
 107          ERR = ERR + ABS(DELTA)
 108          BD(K) = ERR
 109          IND(K) = II
 110   360 CONTINUE
 111       IF (TYPE) GO TO 1001
 112       F = BD(1)
 113       E2(1) = 2.0
 114       BD(1) = F
 115       J = 2
 116   400 DO 500 I = 1, N
 117   500 W(I) = -W(I)
 118       JDEF = -JDEF
 119       GO TO (40,1001), J
 120  1000 IERR = 6 * N + 1
 121  1001 RETURN
 122       END