subroutine qelg(n,epstab,result,abserr,res3la,nres) c***begin prologue qelg c***refer to qagie,qagoe,qagpe,qagse c***routines called r1mach 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 real 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 - real 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 - real c resulting approximation to the integral c c abserr - real c estimate of the absolute error computed from c result and the 3 previous results c c res3la - real 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 qelg c real abserr,delta1,delta2,delta3,r1mach, * 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 c e1 computation of a new element in c e2 the epsilon table is based 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 qelg epmach = r1mach(4) oflow = r1mach(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 = abs(e1) delta2 = e2-e1 err2 = abs(delta2) tol2 = amax1(abs(e2),e1abs)*epmach delta3 = e1-e0 err3 = abs(delta3) tol3 = amax1(e1abs,abs(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 = abs(delta1) tol1 = amax1(e1abs,abs(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.1e+01/delta1+0.1e+01/delta2-0.1e+01/delta3 epsinf = abs(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.1e-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.1e+01/ss epstab(k1) = res k1 = k1-2 error = err2+abs(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 = abs(result-res3la(3))+abs(result-res3la(2)) * +abs(result-res3la(1)) res3la(1) = res3la(2) res3la(2) = res3la(3) res3la(3) = result 100 abserr = amax1(abserr,0.5e+01*epmach*abs(result)) return end