src/numerical/slatec/fortran/dqelg.f

   1 *DECK DQELG
   2       SUBROUTINE DQELG (N, EPSTAB, RESULT, ABSERR, RES3LA, NRES)
   3 C***BEGIN PROLOGUE  DQELG
   4 C***SUBSIDIARY
   5 C***PURPOSE  The routine determines the limit of a given sequence of
   6 C            approximations, by means of the Epsilon algorithm of
   7 C            P.Wynn. An estimate of the absolute error is also given.
   8 C            The condensed Epsilon table is computed. Only those
   9 C            elements needed for the computation of the next diagonal
  10 C            are preserved.
  11 C***LIBRARY   SLATEC
  12 C***TYPE      DOUBLE PRECISION (QELG-S, DQELG-D)
  13 C***KEYWORDS  CONVERGENCE ACCELERATION, EPSILON ALGORITHM, EXTRAPOLATION
  14 C***AUTHOR  Piessens, Robert
  15 C             Applied Mathematics and Programming Division
  16 C             K. U. Leuven
  17 C           de Doncker, Elise
  18 C             Applied Mathematics and Programming Division
  19 C             K. U. Leuven
  20 C***DESCRIPTION
  21 C
  22 C           Epsilon algorithm
  23 C           Standard fortran subroutine
  24 C           Double precision version
  25 C
  26 C           PARAMETERS
  27 C              N      - Integer
  28 C                       EPSTAB(N) contains the new element in the
  29 C                       first column of the epsilon table.
  30 C
  31 C              EPSTAB - Double precision
  32 C                       Vector of dimension 52 containing the elements
  33 C                       of the two lower diagonals of the triangular
  34 C                       epsilon table. The elements are numbered
  35 C                       starting at the right-hand corner of the
  36 C                       triangle.
  37 C
  38 C              RESULT - Double precision
  39 C                       Resulting approximation to the integral
  40 C
  41 C              ABSERR - Double precision
  42 C                       Estimate of the absolute error computed from
  43 C                       RESULT and the 3 previous results
  44 C
  45 C              RES3LA - Double precision
  46 C                       Vector of dimension 3 containing the last 3
  47 C                       results
  48 C
  49 C              NRES   - Integer
  50 C                       Number of calls to the routine
  51 C                       (should be zero at first call)
  52 C
  53 C***SEE ALSO  DQAGIE, DQAGOE, DQAGPE, DQAGSE
  54 C***ROUTINES CALLED  D1MACH
  55 C***REVISION HISTORY  (YYMMDD)
  56 C   800101  DATE WRITTEN
  57 C   890531  Changed all specific intrinsics to generic.  (WRB)
  58 C   890531  REVISION DATE from Version 3.2
  59 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  60 C   900328  Added TYPE section.  (WRB)
  61 C***END PROLOGUE  DQELG
  62 C
  63       DOUBLE PRECISION ABSERR,DELTA1,DELTA2,DELTA3,D1MACH,
  64      1  EPMACH,EPSINF,EPSTAB,ERROR,ERR1,ERR2,ERR3,E0,E1,E1ABS,E2,E3,
  65      2  OFLOW,RES,RESULT,RES3LA,SS,TOL1,TOL2,TOL3
  66       INTEGER I,IB,IB2,IE,INDX,K1,K2,K3,LIMEXP,N,NEWELM,NRES,NUM
  67       DIMENSION EPSTAB(52),RES3LA(3)
  68 C
  69 C           LIST OF MAJOR VARIABLES
  70 C           -----------------------
  71 C
  72 C           E0     - THE 4 ELEMENTS ON WHICH THE COMPUTATION OF A NEW
  73 C           E1       ELEMENT IN THE EPSILON TABLE IS BASED
  74 C           E2
  75 C           E3                 E0
  76 C                        E3    E1    NEW
  77 C                              E2
  78 C           NEWELM - NUMBER OF ELEMENTS TO BE COMPUTED IN THE NEW
  79 C                    DIAGONAL
  80 C           ERROR  - ERROR = ABS(E1-E0)+ABS(E2-E1)+ABS(NEW-E2)
  81 C           RESULT - THE ELEMENT IN THE NEW DIAGONAL WITH LEAST VALUE
  82 C                    OF ERROR
  83 C
  84 C           MACHINE DEPENDENT CONSTANTS
  85 C           ---------------------------
  86 C
  87 C           EPMACH IS THE LARGEST RELATIVE SPACING.
  88 C           OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
  89 C           LIMEXP IS THE MAXIMUM NUMBER OF ELEMENTS THE EPSILON
  90 C           TABLE CAN CONTAIN. IF THIS NUMBER IS REACHED, THE UPPER
  91 C           DIAGONAL OF THE EPSILON TABLE IS DELETED.
  92 C
  93 C***FIRST EXECUTABLE STATEMENT  DQELG
  94       EPMACH = D1MACH(4)
  95       OFLOW = D1MACH(2)
  96       NRES = NRES+1
  97       ABSERR = OFLOW
  98       RESULT = EPSTAB(N)
  99       IF(N.LT.3) GO TO 100
 100       LIMEXP = 50
 101       EPSTAB(N+2) = EPSTAB(N)
 102       NEWELM = (N-1)/2
 103       EPSTAB(N) = OFLOW
 104       NUM = N
 105       K1 = N
 106       DO 40 I = 1,NEWELM
 107         K2 = K1-1
 108         K3 = K1-2
 109         RES = EPSTAB(K1+2)
 110         E0 = EPSTAB(K3)
 111         E1 = EPSTAB(K2)
 112         E2 = RES
 113         E1ABS = ABS(E1)
 114         DELTA2 = E2-E1
 115         ERR2 = ABS(DELTA2)
 116         TOL2 = MAX(ABS(E2),E1ABS)*EPMACH
 117         DELTA3 = E1-E0
 118         ERR3 = ABS(DELTA3)
 119         TOL3 = MAX(E1ABS,ABS(E0))*EPMACH
 120         IF(ERR2.GT.TOL2.OR.ERR3.GT.TOL3) GO TO 10
 121 C
 122 C           IF E0, E1 AND E2 ARE EQUAL TO WITHIN MACHINE
 123 C           ACCURACY, CONVERGENCE IS ASSUMED.
 124 C           RESULT = E2
 125 C           ABSERR = ABS(E1-E0)+ABS(E2-E1)
 126 C
 127         RESULT = RES
 128         ABSERR = ERR2+ERR3
 129 C ***JUMP OUT OF DO-LOOP
 130         GO TO 100
 131    10   E3 = EPSTAB(K1)
 132         EPSTAB(K1) = E1
 133         DELTA1 = E1-E3
 134         ERR1 = ABS(DELTA1)
 135         TOL1 = MAX(E1ABS,ABS(E3))*EPMACH
 136 C
 137 C           IF TWO ELEMENTS ARE VERY CLOSE TO EACH OTHER, OMIT
 138 C           A PART OF THE TABLE BY ADJUSTING THE VALUE OF N
 139 C
 140         IF(ERR1.LE.TOL1.OR.ERR2.LE.TOL2.OR.ERR3.LE.TOL3) GO TO 20
 141         SS = 0.1D+01/DELTA1+0.1D+01/DELTA2-0.1D+01/DELTA3
 142         EPSINF = ABS(SS*E1)
 143 C
 144 C           TEST TO DETECT IRREGULAR BEHAVIOUR IN THE TABLE, AND
 145 C           EVENTUALLY OMIT A PART OF THE TABLE ADJUSTING THE VALUE
 146 C           OF N.
 147 C
 148         IF(EPSINF.GT.0.1D-03) GO TO 30
 149    20   N = I+I-1
 150 C ***JUMP OUT OF DO-LOOP
 151         GO TO 50
 152 C
 153 C           COMPUTE A NEW ELEMENT AND EVENTUALLY ADJUST
 154 C           THE VALUE OF RESULT.
 155 C
 156    30   RES = E1+0.1D+01/SS
 157         EPSTAB(K1) = RES
 158         K1 = K1-2
 159         ERROR = ERR2+ABS(RES-E2)+ERR3
 160         IF(ERROR.GT.ABSERR) GO TO 40
 161         ABSERR = ERROR
 162         RESULT = RES
 163    40 CONTINUE
 164 C
 165 C           SHIFT THE TABLE.
 166 C
 167    50 IF(N.EQ.LIMEXP) N = 2*(LIMEXP/2)-1
 168       IB = 1
 169       IF((NUM/2)*2.EQ.NUM) IB = 2
 170       IE = NEWELM+1
 171       DO 60 I=1,IE
 172         IB2 = IB+2
 173         EPSTAB(IB) = EPSTAB(IB2)
 174         IB = IB2
 175    60 CONTINUE
 176       IF(NUM.EQ.N) GO TO 80
 177       INDX = NUM-N+1
 178       DO 70 I = 1,N
 179         EPSTAB(I)= EPSTAB(INDX)
 180         INDX = INDX+1
 181    70 CONTINUE
 182    80 IF(NRES.GE.4) GO TO 90
 183       RES3LA(NRES) = RESULT
 184       ABSERR = OFLOW
 185       GO TO 100
 186 C
 187 C           COMPUTE ERROR ESTIMATE
 188 C
 189    90 ABSERR = ABS(RESULT-RES3LA(3))+ABS(RESULT-RES3LA(2))
 190      1  +ABS(RESULT-RES3LA(1))
 191       RES3LA(1) = RES3LA(2)
 192       RES3LA(2) = RES3LA(3)
 193       RES3LA(3) = RESULT
 194   100 ABSERR = MAX(ABSERR,0.5D+01*EPMACH*ABS(RESULT))
 195       RETURN
 196       END