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