GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dqelg.f
Go to the documentation of this file.
1  SUBROUTINE dqelg(N,EPSTAB,RESULT,ABSERR,RES3LA,NRES)
2 C***BEGIN PROLOGUE DQELG
3 C***REFER TO DQAGIE,DQAGOE,DQAGPE,DQAGSE
4 C***ROUTINES CALLED D1MACH
5 C***REVISION DATE 830518 (YYMMDD)
6 C***KEYWORDS EPSILON ALGORITHM, CONVERGENCE ACCELERATION,
7 C EXTRAPOLATION
8 C***AUTHOR PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
9 C DE DONCKER,ELISE,APPL. MATH & PROGR. DIV. - K.U.LEUVEN
10 C***PURPOSE THE ROUTINE DETERMINES THE LIMIT OF A GIVEN SEQUENCE OF
11 C APPROXIMATIONS, BY MEANS OF THE EPSILON ALGORITHM OF
12 C P.WYNN. AN ESTIMATE OF THE ABSOLUTE ERROR IS ALSO GIVEN.
13 C THE CONDENSED EPSILON TABLE IS COMPUTED. ONLY THOSE
14 C ELEMENTS NEEDED FOR THE COMPUTATION OF THE NEXT DIAGONAL
15 C ARE PRESERVED.
16 C***DESCRIPTION
17 C
18 C EPSILON ALGORITHM
19 C STANDARD FORTRAN SUBROUTINE
20 C DOUBLE PRECISION VERSION
21 C
22 C PARAMETERS
23 C N - INTEGER
24 C EPSTAB(N) CONTAINS THE NEW ELEMENT IN THE
25 C FIRST COLUMN OF THE EPSILON TABLE.
26 C
27 C EPSTAB - DOUBLE PRECISION
28 C VECTOR OF DIMENSION 52 CONTAINING THE ELEMENTS
29 C OF THE TWO LOWER DIAGONALS OF THE TRIANGULAR
30 C EPSILON TABLE. THE ELEMENTS ARE NUMBERED
31 C STARTING AT THE RIGHT-HAND CORNER OF THE
32 C TRIANGLE.
33 C
34 C RESULT - DOUBLE PRECISION
35 C RESULTING APPROXIMATION TO THE INTEGRAL
36 C
37 C ABSERR - DOUBLE PRECISION
38 C ESTIMATE OF THE ABSOLUTE ERROR COMPUTED FROM
39 C RESULT AND THE 3 PREVIOUS RESULTS
40 C
41 C RES3LA - DOUBLE PRECISION
42 C VECTOR OF DIMENSION 3 CONTAINING THE LAST 3
43 C RESULTS
44 C
45 C NRES - INTEGER
46 C NUMBER OF CALLS TO THE ROUTINE
47 C (SHOULD BE ZERO AT FIRST CALL)
48 C
49 C***END PROLOGUE DQELG
50 C
51  DOUBLE PRECISION ABSERR,DABS,DELTA1,DELTA2,DELTA3,DMAX1,D1MACH,
52  * EPMACH,EPSINF,EPSTAB,ERROR,ERR1,ERR2,ERR3,E0,E1,E1ABS,E2,E3,
53  * OFLOW,RES,RESULT,RES3LA,SS,TOL1,TOL2,TOL3
54  INTEGER I,IB,IB2,IE,INDX,K1,K2,K3,LIMEXP,N,NEWELM,NRES,NUM
55  dimension epstab(52),res3la(3)
56 C
57 C LIST OF MAJOR VARIABLES
58 C -----------------------
59 C
60 C E0 - THE 4 ELEMENTS ON WHICH THE COMPUTATION OF A NEW
61 C E1 ELEMENT IN THE EPSILON TABLE IS BASED
62 C E2
63 C E3 E0
64 C E3 E1 NEW
65 C E2
66 C NEWELM - NUMBER OF ELEMENTS TO BE COMPUTED IN THE NEW
67 C DIAGONAL
68 C ERROR - ERROR = ABS(E1-E0)+ABS(E2-E1)+ABS(NEW-E2)
69 C RESULT - THE ELEMENT IN THE NEW DIAGONAL WITH LEAST VALUE
70 C OF ERROR
71 C
72 C MACHINE DEPENDENT CONSTANTS
73 C ---------------------------
74 C
75 C EPMACH IS THE LARGEST RELATIVE SPACING.
76 C OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
77 C LIMEXP IS THE MAXIMUM NUMBER OF ELEMENTS THE EPSILON
78 C TABLE CAN CONTAIN. IF THIS NUMBER IS REACHED, THE UPPER
79 C DIAGONAL OF THE EPSILON TABLE IS DELETED.
80 C
81 C***FIRST EXECUTABLE STATEMENT DQELG
82  epmach = d1mach(4)
83  oflow = d1mach(2)
84  nres = nres+1
85  abserr = oflow
86  result = epstab(n)
87  IF(n.LT.3) GO TO 100
88  limexp = 50
89  epstab(n+2) = epstab(n)
90  newelm = (n-1)/2
91  epstab(n) = oflow
92  num = n
93  k1 = n
94  DO 40 i = 1,newelm
95  k2 = k1-1
96  k3 = k1-2
97  res = epstab(k1+2)
98  e0 = epstab(k3)
99  e1 = epstab(k2)
100  e2 = res
101  e1abs = dabs(e1)
102  delta2 = e2-e1
103  err2 = dabs(delta2)
104  tol2 = dmax1(dabs(e2),e1abs)*epmach
105  delta3 = e1-e0
106  err3 = dabs(delta3)
107  tol3 = dmax1(e1abs,dabs(e0))*epmach
108  IF(err2.GT.tol2.OR.err3.GT.tol3) GO TO 10
109 C
110 C IF E0, E1 AND E2 ARE EQUAL TO WITHIN MACHINE
111 C ACCURACY, CONVERGENCE IS ASSUMED.
112 C RESULT = E2
113 C ABSERR = ABS(E1-E0)+ABS(E2-E1)
114 C
115  result = res
116  abserr = err2+err3
117 C ***JUMP OUT OF DO-LOOP
118  GO TO 100
119  10 e3 = epstab(k1)
120  epstab(k1) = e1
121  delta1 = e1-e3
122  err1 = dabs(delta1)
123  tol1 = dmax1(e1abs,dabs(e3))*epmach
124 C
125 C IF TWO ELEMENTS ARE VERY CLOSE TO EACH OTHER, OMIT
126 C A PART OF THE TABLE BY ADJUSTING THE VALUE OF N
127 C
128  IF(err1.LE.tol1.OR.err2.LE.tol2.OR.err3.LE.tol3) GO TO 20
129  ss = 0.1d+01/delta1+0.1d+01/delta2-0.1d+01/delta3
130  epsinf = dabs(ss*e1)
131 C
132 C TEST TO DETECT IRREGULAR BEHAVIOUR IN THE TABLE, AND
133 C EVENTUALLY OMIT A PART OF THE TABLE ADJUSTING THE VALUE
134 C OF N.
135 C
136  IF(epsinf.GT.0.1d-03) GO TO 30
137  20 n = i+i-1
138 C ***JUMP OUT OF DO-LOOP
139  GO TO 50
140 C
141 C COMPUTE A NEW ELEMENT AND EVENTUALLY ADJUST
142 C THE VALUE OF RESULT.
143 C
144  30 res = e1+0.1d+01/ss
145  epstab(k1) = res
146  k1 = k1-2
147  error = err2+dabs(res-e2)+err3
148  IF(error.GT.abserr) GO TO 40
149  abserr = error
150  result = res
151  40 CONTINUE
152 C
153 C SHIFT THE TABLE.
154 C
155  50 IF(n.EQ.limexp) n = 2*(limexp/2)-1
156  ib = 1
157  IF((num/2)*2.EQ.num) ib = 2
158  ie = newelm+1
159  DO 60 i=1,ie
160  ib2 = ib+2
161  epstab(ib) = epstab(ib2)
162  ib = ib2
163  60 CONTINUE
164  IF(num.EQ.n) GO TO 80
165  indx = num-n+1
166  DO 70 i = 1,n
167  epstab(i)= epstab(indx)
168  indx = indx+1
169  70 CONTINUE
170  80 IF(nres.GE.4) GO TO 90
171  res3la(nres) = result
172  abserr = oflow
173  GO TO 100
174 C
175 C COMPUTE ERROR ESTIMATE
176 C
177  90 abserr = dabs(result-res3la(3))+dabs(result-res3la(2))
178  * +dabs(result-res3la(1))
179  res3la(1) = res3la(2)
180  res3la(2) = res3la(3)
181  res3la(3) = result
182  100 abserr = dmax1(abserr,0.5d+01*epmach*dabs(result))
183  RETURN
184  END
subroutine dqelg(N, EPSTAB, RESULT, ABSERR, RES3LA, NRES)
Definition: dqelg.f:2