GNU Octave  3.8.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
qelg.f
Go to the documentation of this file.
1  subroutine qelg(n,epstab,result,abserr,res3la,nres)
2 c***begin prologue qelg
3 c***refer to qagie,qagoe,qagpe,qagse
4 c***routines called r1mach
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 real 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 - real
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 - real
35 c resulting approximation to the integral
36 c
37 c abserr - real
38 c estimate of the absolute error computed from
39 c result and the 3 previous results
40 c
41 c res3la - real
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 qelg
50 c
51  real abserr,delta1,delta2,delta3,r1mach,
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
61 c e1 computation of a new element in
62 c e2 the epsilon table is based
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 qelg
82  epmach = r1mach(4)
83  oflow = r1mach(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 = abs(e1)
102  delta2 = e2-e1
103  err2 = abs(delta2)
104  tol2 = amax1(abs(e2),e1abs)*epmach
105  delta3 = e1-e0
106  err3 = abs(delta3)
107  tol3 = amax1(e1abs,abs(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 = abs(delta1)
123  tol1 = amax1(e1abs,abs(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.1e+01/delta1+0.1e+01/delta2-0.1e+01/delta3
130  epsinf = abs(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.1e-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.1e+01/ss
145  epstab(k1) = res
146  k1 = k1-2
147  error = err2+abs(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 = abs(result-res3la(3))+abs(result-res3la(2))
178  * +abs(result-res3la(1))
179  res3la(1) = res3la(2)
180  res3la(2) = res3la(3)
181  res3la(3) = result
182  100 abserr = amax1(abserr,0.5e+01*epmach*abs(result))
183  return
184  end