1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
|
subroutine dqelg(n,epstab,result,abserr,res3la,nres)
c***begin prologue dqelg
c***refer to dqagie,dqagoe,dqagpe,dqagse
c***routines called d1mach
c***revision date 830518 (yymmdd)
c***keywords epsilon algorithm, convergence acceleration,
c extrapolation
c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
c de doncker,elise,appl. math & progr. div. - k.u.leuven
c***purpose the routine determines the limit of a given sequence of
c approximations, by means of the epsilon algorithm of
c p.wynn. an estimate of the absolute error is also given.
c the condensed epsilon table is computed. only those
c elements needed for the computation of the next diagonal
c are preserved.
c***description
c
c epsilon algorithm
c standard fortran subroutine
c double precision version
c
c parameters
c n - integer
c epstab(n) contains the new element in the
c first column of the epsilon table.
c
c epstab - double precision
c vector of dimension 52 containing the elements
c of the two lower diagonals of the triangular
c epsilon table. the elements are numbered
c starting at the right-hand corner of the
c triangle.
c
c result - double precision
c resulting approximation to the integral
c
c abserr - double precision
c estimate of the absolute error computed from
c result and the 3 previous results
c
c res3la - double precision
c vector of dimension 3 containing the last 3
c results
c
c nres - integer
c number of calls to the routine
c (should be zero at first call)
c
c***end prologue dqelg
c
double precision abserr,dabs,delta1,delta2,delta3,dmax1,d1mach,
* epmach,epsinf,epstab,error,err1,err2,err3,e0,e1,e1abs,e2,e3,
* oflow,res,result,res3la,ss,tol1,tol2,tol3
integer i,ib,ib2,ie,indx,k1,k2,k3,limexp,n,newelm,nres,num
dimension epstab(52),res3la(3)
c
c list of major variables
c -----------------------
c
c e0 - the 4 elements on which the computation of a new
c e1 element in the epsilon table is based
c e2
c e3 e0
c e3 e1 new
c e2
c newelm - number of elements to be computed in the new
c diagonal
c error - error = abs(e1-e0)+abs(e2-e1)+abs(new-e2)
c result - the element in the new diagonal with least value
c of error
c
c machine dependent constants
c ---------------------------
c
c epmach is the largest relative spacing.
c oflow is the largest positive magnitude.
c limexp is the maximum number of elements the epsilon
c table can contain. if this number is reached, the upper
c diagonal of the epsilon table is deleted.
c
c***first executable statement dqelg
epmach = d1mach(4)
oflow = d1mach(2)
nres = nres+1
abserr = oflow
result = epstab(n)
if(n.lt.3) go to 100
limexp = 50
epstab(n+2) = epstab(n)
newelm = (n-1)/2
epstab(n) = oflow
num = n
k1 = n
do 40 i = 1,newelm
k2 = k1-1
k3 = k1-2
res = epstab(k1+2)
e0 = epstab(k3)
e1 = epstab(k2)
e2 = res
e1abs = dabs(e1)
delta2 = e2-e1
err2 = dabs(delta2)
tol2 = dmax1(dabs(e2),e1abs)*epmach
delta3 = e1-e0
err3 = dabs(delta3)
tol3 = dmax1(e1abs,dabs(e0))*epmach
if(err2.gt.tol2.or.err3.gt.tol3) go to 10
c
c if e0, e1 and e2 are equal to within machine
c accuracy, convergence is assumed.
c result = e2
c abserr = abs(e1-e0)+abs(e2-e1)
c
result = res
abserr = err2+err3
c ***jump out of do-loop
go to 100
10 e3 = epstab(k1)
epstab(k1) = e1
delta1 = e1-e3
err1 = dabs(delta1)
tol1 = dmax1(e1abs,dabs(e3))*epmach
c
c if two elements are very close to each other, omit
c a part of the table by adjusting the value of n
c
if(err1.le.tol1.or.err2.le.tol2.or.err3.le.tol3) go to 20
ss = 0.1d+01/delta1+0.1d+01/delta2-0.1d+01/delta3
epsinf = dabs(ss*e1)
c
c test to detect irregular behaviour in the table, and
c eventually omit a part of the table adjusting the value
c of n.
c
if(epsinf.gt.0.1d-03) go to 30
20 n = i+i-1
c ***jump out of do-loop
go to 50
c
c compute a new element and eventually adjust
c the value of result.
c
30 res = e1+0.1d+01/ss
epstab(k1) = res
k1 = k1-2
error = err2+dabs(res-e2)+err3
if(error.gt.abserr) go to 40
abserr = error
result = res
40 continue
c
c shift the table.
c
50 if(n.eq.limexp) n = 2*(limexp/2)-1
ib = 1
if((num/2)*2.eq.num) ib = 2
ie = newelm+1
do 60 i=1,ie
ib2 = ib+2
epstab(ib) = epstab(ib2)
ib = ib2
60 continue
if(num.eq.n) go to 80
indx = num-n+1
do 70 i = 1,n
epstab(i)= epstab(indx)
indx = indx+1
70 continue
80 if(nres.ge.4) go to 90
res3la(nres) = result
abserr = oflow
go to 100
c
c compute error estimate
c
90 abserr = dabs(result-res3la(3))+dabs(result-res3la(2))
* +dabs(result-res3la(1))
res3la(1) = res3la(2)
res3la(2) = res3la(3)
res3la(3) = result
100 abserr = dmax1(abserr,0.5d+01*epmach*dabs(result))
return
end
|
