GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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
std::string dimension(void) const
void error(const char *fmt,...)
Definition: error.cc:476
subroutine qelg(n, epstab, result, abserr, res3la, nres)
Definition: qelg.f:1
T abs(T x)
Definition: pr-output.cc:3062