1*DECK DQELG
2      SUBROUTINE DQELG (N, EPSTAB, RESULT, ABSERR, RES3LA, NRES)
3C***BEGIN PROLOGUE  DQELG
4C***SUBSIDIARY
5C***PURPOSE  The routine determines the limit of a given sequence of
6C            approximations, by means of the Epsilon algorithm of
7C            P.Wynn. An estimate of the absolute error is also given.
8C            The condensed Epsilon table is computed. Only those
9C            elements needed for the computation of the next diagonal
10C            are preserved.
11C***LIBRARY   SLATEC
12C***TYPE      DOUBLE PRECISION (QELG-S, DQELG-D)
13C***KEYWORDS  CONVERGENCE ACCELERATION, EPSILON ALGORITHM, EXTRAPOLATION
14C***AUTHOR  Piessens, Robert
15C             Applied Mathematics and Programming Division
16C             K. U. Leuven
17C           de Doncker, Elise
18C             Applied Mathematics and Programming Division
19C             K. U. Leuven
20C***DESCRIPTION
21C
22C           Epsilon algorithm
23C           Standard fortran subroutine
24C           Double precision version
25C
26C           PARAMETERS
27C              N      - Integer
28C                       EPSTAB(N) contains the new element in the
29C                       first column of the epsilon table.
30C
31C              EPSTAB - Double precision
32C                       Vector of dimension 52 containing the elements
33C                       of the two lower diagonals of the triangular
34C                       epsilon table. The elements are numbered
35C                       starting at the right-hand corner of the
36C                       triangle.
37C
38C              RESULT - Double precision
39C                       Resulting approximation to the integral
40C
41C              ABSERR - Double precision
42C                       Estimate of the absolute error computed from
43C                       RESULT and the 3 previous results
44C
45C              RES3LA - Double precision
46C                       Vector of dimension 3 containing the last 3
47C                       results
48C
49C              NRES   - Integer
50C                       Number of calls to the routine
51C                       (should be zero at first call)
52C
53C***SEE ALSO  DQAGIE, DQAGOE, DQAGPE, DQAGSE
54C***ROUTINES CALLED  D1MACH
55C***REVISION HISTORY  (YYMMDD)
56C   800101  DATE WRITTEN
57C   890531  Changed all specific intrinsics to generic.  (WRB)
58C   890531  REVISION DATE from Version 3.2
59C   891214  Prologue converted to Version 4.0 format.  (BAB)
60C   900328  Added TYPE section.  (WRB)
61C***END PROLOGUE  DQELG
62C
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)
68C
69C           LIST OF MAJOR VARIABLES
70C           -----------------------
71C
72C           E0     - THE 4 ELEMENTS ON WHICH THE COMPUTATION OF A NEW
73C           E1       ELEMENT IN THE EPSILON TABLE IS BASED
74C           E2
75C           E3                 E0
76C                        E3    E1    NEW
77C                              E2
78C           NEWELM - NUMBER OF ELEMENTS TO BE COMPUTED IN THE NEW
79C                    DIAGONAL
80C           ERROR  - ERROR = ABS(E1-E0)+ABS(E2-E1)+ABS(NEW-E2)
81C           RESULT - THE ELEMENT IN THE NEW DIAGONAL WITH LEAST VALUE
82C                    OF ERROR
83C
84C           MACHINE DEPENDENT CONSTANTS
85C           ---------------------------
86C
87C           EPMACH IS THE LARGEST RELATIVE SPACING.
88C           OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
89C           LIMEXP IS THE MAXIMUM NUMBER OF ELEMENTS THE EPSILON
90C           TABLE CAN CONTAIN. IF THIS NUMBER IS REACHED, THE UPPER
91C           DIAGONAL OF THE EPSILON TABLE IS DELETED.
92C
93C***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
121C
122C           IF E0, E1 AND E2 ARE EQUAL TO WITHIN MACHINE
123C           ACCURACY, CONVERGENCE IS ASSUMED.
124C           RESULT = E2
125C           ABSERR = ABS(E1-E0)+ABS(E2-E1)
126C
127        RESULT = RES
128        ABSERR = ERR2+ERR3
129C ***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
136C
137C           IF TWO ELEMENTS ARE VERY CLOSE TO EACH OTHER, OMIT
138C           A PART OF THE TABLE BY ADJUSTING THE VALUE OF N
139C
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)
143C
144C           TEST TO DETECT IRREGULAR BEHAVIOUR IN THE TABLE, AND
145C           EVENTUALLY OMIT A PART OF THE TABLE ADJUSTING THE VALUE
146C           OF N.
147C
148        IF(EPSINF.GT.0.1D-03) GO TO 30
149   20   N = I+I-1
150C ***JUMP OUT OF DO-LOOP
151        GO TO 50
152C
153C           COMPUTE A NEW ELEMENT AND EVENTUALLY ADJUST
154C           THE VALUE OF RESULT.
155C
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
164C
165C           SHIFT THE TABLE.
166C
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
186C
187C           COMPUTE ERROR ESTIMATE
188C
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
197