| [f,a,b,npts2,points,epsabs,epsrel,limit,resultf,abserr,neval,ier,alist,blist,rlist,elist,pts,iord,level,ndin,last]=dqagpe(f,a,b,npts2,points,epsabs,epsrel,limit,resultf,abserr,neval,ier,alist,blist,rlist,elist,pts,iord,level,ndin,last); |
function [f,a,b,npts2,points,epsabs,epsrel,limit,resultf,abserr,neval,ier,alist,blist,rlist,elist,pts,iord,level,ndin,last]=dqagpe(f,a,b,npts2,points,epsabs,epsrel,limit,resultf,abserr,neval,ier,alist,blist,rlist,elist,pts,iord,level,ndin,last);
%***BEGIN PROLOGUE DQAGPE
%***PURPOSE Approximate a given definite integral I = Integral of F
% over (A,B), hopefully satisfying the accuracy claim:
% ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
% Break points of the integration interval, where local
% difficulties of the integrand may occur (e.g. singularities
% or discontinuities) are provided by the user.
%***LIBRARY SLATEC (QUADPACK)
%***CATEGORY H2A2A1
%***TYPE doubleprecision (QAGPE-S, DQAGPE-D)
%***KEYWORDS AUTOMATIC INTEGRATOR, EXTRAPOLATION, GENERAL-PURPOSE,
% GLOBALLY ADAPTIVE, QUADPACK, QUADRATURE,
% SINGULARITIES AT USER SPECIFIED POINTS
%***AUTHOR Piessens, Robert
% Applied Mathematics and Programming Division
% K. U. Leuven
% de Doncker, Elise
% Applied Mathematics and Programming Division
% K. U. Leuven
%***DESCRIPTION
%
% Computation of a definite integral
% Standard fortran subroutine
% doubleprecision version
%
% PARAMETERS
% ON ENTRY
% F - doubleprecision
% function subprogram defining the integrand
% function F(X). The actual name for F needs to be
% declared E X T E R N A L in the driver program.
%
% A - doubleprecision
% Lower limit of integration
%
% B - doubleprecision
% Upper limit of integration
%
% NPTS2 - Integer
% Number equal to two more than the number of
% user-supplied break points within the integration
% range, NPTS2.GE.2.
% If NPTS2.LT.2, the routine will end with IER = 6.
%
% POINTS - doubleprecision
% Vector of dimension NPTS2, the first (NPTS2-2)
% elements of which are the user provided break
% POINTS. If these POINTS do not constitute an
% ascending sequence there will be an automatic
% sorting.
%
% EPSABS - doubleprecision
% Absolute accuracy requested
% EPSREL - doubleprecision
% Relative accuracy requested
% If EPSABS.LE.0
% and EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
% the routine will end with IER = 6.
%
% LIMIT - Integer
% Gives an upper bound on the number of subintervals
% in the partition of (A,B), LIMIT.GE.NPTS2
% If LIMIT.LT.NPTS2, the routine will end with
% IER = 6.
%
% ON RETURN
% RESULTF - doubleprecision
% Approximation to the integral
%
% ABSERR - doubleprecision
% Estimate of the modulus of the absolute error,
% which should equal or exceed ABS(I-RESULT)
%
% NEVAL - Integer
% Number of integrand evaluations
%
% IER - Integer
% IER = 0 Normal and reliable termination of the
% routine. It is assumed that the requested
% accuracy has been achieved.
% IER.GT.0 Abnormal termination of the routine.
% The estimates for integral and error are
% less reliable. It is assumed that the
% requested accuracy has not been achieved.
% ERROR MESSAGES
% IER = 1 Maximum number of subdivisions allowed
% has been achieved. One can allow more
% subdivisions by increasing the value of
% LIMIT (and taking the according dimension
% adjustments into account). However, if
% this yields no improvement it is advised
% to analyze the integrand in order to
% determine the integration difficulties. If
% the position of a local difficulty can be
% determined (i.e. SINGULARITY,
% DISCONTINUITY within the interval), it
% should be supplied to the routine as an
% element of the vector points. If necessary
% an appropriate special-purpose integrator
% must be used, which is designed for
% handling the type of difficulty involved.
% = 2 The occurrence of roundoff error is
% detected, which prevents the requested
% tolerance from being achieved.
% The error may be under-estimated.
% = 3 Extremely bad integrand behaviour occurs
% At some points of the integration
% interval.
% = 4 The algorithm does not converge.
% Roundoff error is detected in the
% extrapolation table. It is presumed that
% the requested tolerance cannot be
% achieved, and that the returned result is
% the best which can be obtained.
% = 5 The integral is probably divergent, or
% slowly convergent. It must be noted that
% divergence can occur with any other value
% of IER.GT.0.
% = 6 The input is invalid because
% NPTS2.LT.2 or
% Break points are specified outside
% the integration range or
% (EPSABS.LE.0 and
% EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))
% or LIMIT.LT.NPTS2.
% RESULTF, ABSERR, NEVAL, LAST, RLIST(1),
% and ELIST(1) are set to zero. ALIST(1) and
% BLIST(1) are set to A and B respectively.
%
% ALIST - doubleprecision
% Vector of dimension at least LIMIT, the first
% LAST elements of which are the left end points
% of the subintervals in the partition of the given
% integration range (A,B)
%
% BLIST - doubleprecision
% Vector of dimension at least LIMIT, the first
% LAST elements of which are the right end points
% of the subintervals in the partition of the given
% integration range (A,B)
%
% RLIST - doubleprecision
% Vector of dimension at least LIMIT, the first
% LAST elements of which are the integral
% approximations on the subintervals
%
% ELIST - doubleprecision
% Vector of dimension at least LIMIT, the first
% LAST elements of which are the moduli of the
% absolute error estimates on the subintervals
%
% PTS - doubleprecision
% Vector of dimension at least NPTS2, containing the
% integration limits and the break points of the
% interval in ascending sequence.
%
% LEVEL - Integer
% Vector of dimension at least LIMIT, containing the
% subdivision levels of the subinterval, i.e. if
% (AA,BB) is a subinterval of (P1,P2) where P1 as
% well as P2 is a user-provided break point or
% integration limit, then (AA,BB) has level L if
% ABS(BB-AA) = ABS(P2-P1)*2**(-L).
%
% NDIN - Integer
% Vector of dimension at least NPTS2, after first
% integration over the intervals (PTS(I)),PTS(I+1),
% I = 0,1, ..., NPTS2-2, the error estimates over
% some of the intervals may have been increased
% artificially, in order to put their subdivision
% forward. If this happens for the subinterval
% numbered K, NDIN(K) is put to 1, otherwise
% NDIN(K) = 0.
%
% IORD - Integer
% Vector of dimension at least LIMIT, the first K
% elements of which are pointers to the
% error estimates over the subintervals,
% such that ELIST(IORD(1)), ..., ELIST(IORD(K))
% form a decreasing sequence, with K = LAST
% If LAST.LE.(LIMIT/2+2), and K = LIMIT+1-LAST
% otherwise
%
% LAST - Integer
% Number of subintervals actually produced in the
% subdivisions process
%
%***REFERENCES (NONE)
%***ROUTINES CALLED D1MACH, DQELG, DQK21, DQPSRT
%***REVISION HISTORY (YYMMDD)
% 800101 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890831 Modified array declarations. (WRB)
% 890831 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
%***end PROLOGUE DQAGPE
persistent a1 a2 abseps area area1 area12 area2 b1 b2 correc defab1 defab2 defabs dres epmach erlarg erlast errbnd errmax erro12 error1 error2 errsum ertest extrap gt i id ierro ind1 ind2 ip1 iroff1 iroff2 iroff3 j jlow jupbnd k ksgn ktmin levcur levmax maxerr nintmlv nintp1 noext npts nres nrmax numrl2 oflow res3la resa resabs reseps rlist2 signmlv temp uflow ;
if isempty(abseps), abseps=0; end;
if isempty(area), area=0; end;
if isempty(area1), area1=0; end;
if isempty(area12), area12=0; end;
if isempty(area2), area2=0; end;
if isempty(a1), a1=0; end;
if isempty(a2), a2=0; end;
if isempty(b1), b1=0; end;
if isempty(b2), b2=0; end;
if isempty(correc), correc=0; end;
if isempty(defabs), defabs=0; end;
if isempty(defab1), defab1=0; end;
if isempty(defab2), defab2=0; end;
if isempty(dres), dres=0; end;
if isempty(epmach), epmach=0; end;
if isempty(erlarg), erlarg=0; end;
if isempty(erlast), erlast=0; end;
if isempty(errbnd), errbnd=0; end;
if isempty(errmax), errmax=0; end;
if isempty(error1), error1=0; end;
if isempty(erro12), erro12=0; end;
if isempty(error2), error2=0; end;
if isempty(errsum), errsum=0; end;
if isempty(ertest), ertest=0; end;
if isempty(oflow), oflow=0; end;
if isempty(resa), resa=0; end;
if isempty(resabs), resabs=0; end;
if isempty(reseps), reseps=0; end;
if isempty(res3la), res3la=zeros(1,3); end;
if isempty(rlist2), rlist2=zeros(1,52); end;
if isempty(signmlv), signmlv=0; end;
if isempty(temp), temp=0; end;
if isempty(uflow), uflow=0; end;
if isempty(i), i=0; end;
if isempty(id), id=0; end;
if isempty(ierro), ierro=0; end;
if isempty(ind1), ind1=0; end;
if isempty(ind2), ind2=0; end;
if isempty(ip1), ip1=0; end;
if isempty(iroff1), iroff1=0; end;
if isempty(iroff2), iroff2=0; end;
if isempty(iroff3), iroff3=0; end;
if isempty(j), j=0; end;
if isempty(jlow), jlow=0; end;
if isempty(jupbnd), jupbnd=0; end;
if isempty(k), k=0; end;
if isempty(ksgn), ksgn=0; end;
if isempty(ktmin), ktmin=0; end;
if isempty(levcur), levcur=0; end;
if isempty(levmax), levmax=0; end;
if isempty(maxerr), maxerr=0; end;
if isempty(nintmlv), nintmlv=0; end;
if isempty(nintp1), nintp1=0; end;
if isempty(npts), npts=0; end;
if isempty(nres), nres=0; end;
if isempty(nrmax), nrmax=0; end;
if isempty(numrl2), numrl2=0; end;
if isempty(gt), gt=zeros(1,3); end;
if isempty(extrap), extrap=false; end;
if isempty(noext), noext=false; end;
%
%
alist_shape=size(alist);alist=reshape(alist,1,[]);
blist_shape=size(blist);blist=reshape(blist,1,[]);
elist_shape=size(elist);elist=reshape(elist,1,[]);
iord_shape=size(iord);iord=reshape(iord,1,[]);
level_shape=size(level);level=reshape(level,1,[]);
ndin_shape=size(ndin);ndin=reshape(ndin,1,[]);
points_shape=size(points);points=reshape(points,1,[]);
pts_shape=size(pts);pts=reshape(pts,1,[]);
rlist_shape=size(rlist);rlist=reshape(rlist,1,[]);
%
%
% THE DIMENSION OF RLIST2 IS DETERMINED BY THE VALUE OF
% LIMEXP IN SUBROUTINE EPSALG (RLIST2 SHOULD BE OF DIMENSION
% (LIMEXP+2) AT LEAST).
%
%
% LIST OF MAJOR VARIABLES
% -----------------------
%
% ALIST - LIST OF LEFT END POINTS OF ALL SUBINTERVALS
% CONSIDERED UP TO NOW
% BLIST - LIST OF RIGHT END POINTS OF ALL SUBINTERVALS
% CONSIDERED UP TO NOW
% RLIST(I) - APPROXIMATION TO THE INTEGRAL OVER
% (ALIST(I),BLIST(I))
% RLIST2 - ARRAY OF DIMENSION AT LEAST LIMEXP+2
% CONTAINING THE PART OF THE EPSILON TABLE WHICH
% IS STILL NEEDED FOR FURTHER COMPUTATIONS
% ELIST(I) - ERROR ESTIMATE APPLYING TO RLIST(I)
% MAXERR - POINTER TO THE INTERVAL WITH LARGEST ERROR
% ESTIMATE
% ERRMAX - ELIST(MAXERR)
% ERLAST - ERROR ON THE INTERVAL CURRENTLY SUBDIVIDED
% (BEFORE THAT SUBDIVISION HAS TAKEN PLACE)
% AREA - SUM OF THE INTEGRALS OVER THE SUBINTERVALS
% ERRSUM - SUM OF THE ERRORS OVER THE SUBINTERVALS
% ERRBND - REQUESTED ACCURACY MAX(EPSABS,EPSREL*
% ABS(RESULTF))
% *****1 - VARIABLE FOR THE LEFT SUBINTERVAL
% *****2 - VARIABLE FOR THE RIGHT SUBINTERVAL
% LAST - INDEX FOR SUBDIVISION
% NRES - NUMBER OF CALLS TO THE EXTRAPOLATION ROUTINE
% NUMRL2 - NUMBER OF ELEMENTS IN RLIST2. IF AN APPROPRIATE
% APPROXIMATION TO THE COMPOUNDED INTEGRAL HAS
% BEEN OBTAINED, IT IS PUT IN RLIST2(NUMRL2) AFTER
% NUMRL2 HAS BEEN INCREASED BY ONE.
% ERLARG - SUM OF THE ERRORS OVER THE INTERVALS LARGER
% THAN THE SMALLEST INTERVAL CONSIDERED UP TO NOW
% EXTRAP - LOGICAL VARIABLE DENOTING THAT THE ROUTINE
% IS ATTEMPTING TO PERFORM EXTRAPOLATION. I.E.
% BEFORE SUBDIVIDING THE SMALLEST INTERVAL WE
% TRY TO DECREASE THE VALUE OF ERLARG.
% NOEXT - LOGICAL VARIABLE DENOTING THAT EXTRAPOLATION IS
% NO LONGER ALLOWED (truemlv-VALUE)
%
% MACHINE DEPENDENT CONSTANTS
% ---------------------------
%
% EPMACH IS THE LARGEST RELATIVE SPACING.
% UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
% OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
%
%***FIRST EXECUTABLE STATEMENT DQAGPE
gt(:)=0;
[epmach ]=d1mach(4);
%
% TEST ON VALIDITY OF PARAMETERS
% -----------------------------
%
ier = 0;
neval = 0;
last = 0;
resultf = 0.0d+00;
abserr = 0.0d+00;
alist(1) = a;
blist(1) = b;
rlist(1) = 0.0d+00;
elist(1) = 0.0d+00;
iord(1) = 0;
level(1) = 0;
npts = fix(npts2 - 2);
if( npts2<2 || limit<=npts ||(epsabs<=0.0d+00 && epsrel<max(0.5d+02.*epmach,0.5d-28)) )
ier = 6;
end;
if( ier==6 )
alist_shape=zeros(alist_shape);alist_shape(:)=alist(1:numel(alist_shape));alist=alist_shape;
blist_shape=zeros(blist_shape);blist_shape(:)=blist(1:numel(blist_shape));blist=blist_shape;
elist_shape=zeros(elist_shape);elist_shape(:)=elist(1:numel(elist_shape));elist=elist_shape;
iord_shape=zeros(iord_shape);iord_shape(:)=iord(1:numel(iord_shape));iord=iord_shape;
level_shape=zeros(level_shape);level_shape(:)=level(1:numel(level_shape));level=level_shape;
ndin_shape=zeros(ndin_shape);ndin_shape(:)=ndin(1:numel(ndin_shape));ndin=ndin_shape;
points_shape=zeros(points_shape);points_shape(:)=points(1:numel(points_shape));points=points_shape;
pts_shape=zeros(pts_shape);pts_shape(:)=pts(1:numel(pts_shape));pts=pts_shape;
rlist_shape=zeros(rlist_shape);rlist_shape(:)=rlist(1:numel(rlist_shape));rlist=rlist_shape;
return;
end;
%
% IF ANY BREAK POINTS ARE PROVIDED, SORT THEM INTO AN
% ASCENDING SEQUENCE.
%
signmlv = 1.0d+00;
if( a>b )
signmlv = -1.0d+00;
end;
pts(1) = min(a,b);
if( npts~=0 )
for i = 1 : npts;
pts(i+1) = points(i);
end; i = fix(npts+1);
end;
pts(npts+2) = max(a,b);
nintmlv = npts + 1;
a1 = pts(1);
if( npts~=0 )
nintp1 = fix(nintmlv + 1);
for i = 1 : nintmlv;
ip1 = fix(i + 1);
for j = ip1 : nintp1;
if( pts(i)>pts(j) )
temp = pts(i);
pts(i) = pts(j);
pts(j) = temp;
end;
end; j = fix(nintp1+1);
end; i = fix(nintmlv+1);
if( pts(1)~=min(a,b) || pts(nintp1)~=max(a,b) )
ier = 6;
end;
if( ier==6 )
alist_shape=zeros(alist_shape);alist_shape(:)=alist(1:numel(alist_shape));alist=alist_shape;
blist_shape=zeros(blist_shape);blist_shape(:)=blist(1:numel(blist_shape));blist=blist_shape;
elist_shape=zeros(elist_shape);elist_shape(:)=elist(1:numel(elist_shape));elist=elist_shape;
iord_shape=zeros(iord_shape);iord_shape(:)=iord(1:numel(iord_shape));iord=iord_shape;
level_shape=zeros(level_shape);level_shape(:)=level(1:numel(level_shape));level=level_shape;
ndin_shape=zeros(ndin_shape);ndin_shape(:)=ndin(1:numel(ndin_shape));ndin=ndin_shape;
points_shape=zeros(points_shape);points_shape(:)=points(1:numel(points_shape));points=points_shape;
pts_shape=zeros(pts_shape);pts_shape(:)=pts(1:numel(pts_shape));pts=pts_shape;
rlist_shape=zeros(rlist_shape);rlist_shape(:)=rlist(1:numel(rlist_shape));rlist=rlist_shape;
return;
end;
end;
%
% COMPUTE FIRST INTEGRAL AND ERROR APPROXIMATIONS.
% ------------------------------------------------
%
resabs = 0.0d+00;
for i = 1 : nintmlv;
b1 = pts(i+1);
[f,a1,b1,area1,error1,defabs,resa]=dqk21(f,a1,b1,area1,error1,defabs,resa);
abserr = abserr + error1;
resultf = resultf + area1;
ndin(i) = 0;
if( error1==resa && error1~=0.0d+00 )
ndin(i) = 1;
end;
resabs = resabs + defabs;
level(i) = 0;
elist(i) = error1;
alist(i) = a1;
blist(i) = b1;
rlist(i) = area1;
iord(i) = fix(i);
a1 = b1;
end; i = fix(nintmlv+1);
errsum = 0.0d+00;
for i = 1 : nintmlv;
if( ndin(i)==1 )
elist(i) = abserr;
end;
errsum = errsum + elist(i);
end; i = fix(nintmlv+1);
%
% TEST ON ACCURACY.
%
last = fix(nintmlv);
neval = fix(21.*nintmlv);
dres = abs(resultf);
errbnd = max(epsabs,epsrel.*dres);
if( abserr<=0.1d+03.*epmach.*resabs && abserr>errbnd )
ier = 2;
end;
if( nintmlv~=1 )
for i = 1 : npts;
jlow = fix(i + 1);
ind1 = fix(iord(i));
for j = jlow : nintmlv;
ind2 = fix(iord(j));
if( elist(ind1)<=elist(ind2) )
ind1 = fix(ind2);
k = fix(j);
end;
end; j = fix(nintmlv+1);
if( ind1~=iord(i) )
iord(k) = fix(iord(i));
iord(i) = fix(ind1);
end;
end; i = fix(npts+1);
if( limit<npts2 )
ier = 1;
end;
end;
if( ier~=0 || abserr<=errbnd )
alist_shape=zeros(alist_shape);alist_shape(:)=alist(1:numel(alist_shape));alist=alist_shape;
blist_shape=zeros(blist_shape);blist_shape(:)=blist(1:numel(blist_shape));blist=blist_shape;
elist_shape=zeros(elist_shape);elist_shape(:)=elist(1:numel(elist_shape));elist=elist_shape;
iord_shape=zeros(iord_shape);iord_shape(:)=iord(1:numel(iord_shape));iord=iord_shape;
level_shape=zeros(level_shape);level_shape(:)=level(1:numel(level_shape));level=level_shape;
ndin_shape=zeros(ndin_shape);ndin_shape(:)=ndin(1:numel(ndin_shape));ndin=ndin_shape;
points_shape=zeros(points_shape);points_shape(:)=points(1:numel(points_shape));points=points_shape;
pts_shape=zeros(pts_shape);pts_shape(:)=pts(1:numel(pts_shape));pts=pts_shape;
rlist_shape=zeros(rlist_shape);rlist_shape(:)=rlist(1:numel(rlist_shape));rlist=rlist_shape;
return;
end;
%
% INITIALIZATION
% --------------
%
rlist2(1) = resultf;
maxerr = fix(iord(1));
errmax = elist(maxerr);
area = resultf;
nrmax = 1;
nres = 0;
numrl2 = 1;
ktmin = 0;
extrap = false;
noext = false;
erlarg = errsum;
ertest = errbnd;
levmax = 1;
iroff1 = 0;
iroff2 = 0;
iroff3 = 0;
ierro = 0;
[uflow ]=d1mach(1);
[oflow ]=d1mach(2);
abserr = oflow;
ksgn = -1;
if( dres>=(0.1d+01-0.5d+02.*epmach).*resabs )
ksgn = 1;
end;
%
% MAIN DO-LOOP
% ------------
%
for last = npts2 : limit;
gt(:)=0;
%
% BISECT THE SUBINTERVAL WITH THE NRMAX-TH LARGEST ERROR
% ESTIMATE.
%
levcur = fix(level(maxerr) + 1);
a1 = alist(maxerr);
b1 = 0.5d+00.*(alist(maxerr)+blist(maxerr));
a2 = b1;
b2 = blist(maxerr);
erlast = errmax;
[f,a1,b1,area1,error1,resa,defab1]=dqk21(f,a1,b1,area1,error1,resa,defab1);
[f,a2,b2,area2,error2,resa,defab2]=dqk21(f,a2,b2,area2,error2,resa,defab2);
%
% IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
% AND ERROR AND TEST FOR ACCURACY.
%
neval = fix(neval + 42);
area12 = area1 + area2;
erro12 = error1 + error2;
errsum = errsum + erro12 - errmax;
area = area + area12 - rlist(maxerr);
if( defab1~=error1 && defab2~=error2 )
if( abs(rlist(maxerr)-area12)<=0.1d-04.*abs(area12) &&erro12>=0.99d+00.*errmax )
if( extrap )
iroff2 = fix(iroff2 + 1);
end;
if( ~extrap )
iroff1 = fix(iroff1 + 1);
end;
end;
if( last>10 && erro12>errmax )
iroff3 = fix(iroff3 + 1);
end;
end;
level(maxerr) = fix(levcur);
level(last) = fix(levcur);
rlist(maxerr) = area1;
rlist(last) = area2;
errbnd = max(epsabs,epsrel.*abs(area));
%
% TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG.
%
if( iroff1+iroff2>=10 || iroff3>=20 )
ier = 2;
end;
if( iroff2>=5 )
ierro = 3;
end;
%
% SET ERROR FLAG IN THE CASE THAT THE NUMBER OF
% SUBINTERVALS EQUALS LIMIT.
%
if( last==limit )
ier = 1;
end;
%
% SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
% AT A POINT OF THE INTEGRATION RANGE
%
if( max(abs(a1),abs(b2))<=(0.1d+01+0.1d+03.*epmach).*(abs(a2)+0.1d+04.*uflow) )
ier = 4;
end;
%
% APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
%
if( error2>error1 )
alist(maxerr) = a2;
alist(last) = a1;
blist(last) = b1;
rlist(maxerr) = area2;
rlist(last) = area1;
elist(maxerr) = error2;
elist(last) = error1;
else;
alist(last) = a2;
blist(maxerr) = b1;
blist(last) = b2;
elist(maxerr) = error1;
elist(last) = error2;
end;
%
% CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING
% IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
% WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT).
%
[limit,last,maxerr,errmax,elist,iord,nrmax]=dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax);
% ***JUMP OUT OF DO-LOOP
if( errsum<=errbnd )
gt(2)=1;
break;
end;
% ***JUMP OUT OF DO-LOOP
if( ier~=0 )
break;
end;
if( ~(noext) )
erlarg = erlarg - erlast;
if( levcur+1<=levmax )
erlarg = erlarg + erro12;
end;
if( ~(extrap) )
%
% TEST WHETHER THE INTERVAL TO BE BISECTED NEXT IS THE
% SMALLEST INTERVAL.
%
if( level(maxerr)+1<=levmax )
continue;
end;
extrap = true;
nrmax = 2;
end;
if( ierro~=3 && erlarg>ertest )
%
% THE SMALLEST INTERVAL HAS THE LARGEST ERROR.
% BEFORE BISECTING DECREASE THE SUM OF THE ERRORS OVER
% THE LARGER INTERVALS (ERLARG) AND PERFORM EXTRAPOLATION.
%
id = fix(nrmax);
jupbnd = fix(last);
if( last>(2+fix(limit./2)) )
jupbnd = fix(limit + 3 - last);
end;
for k = id : jupbnd;
maxerr = fix(iord(nrmax));
errmax = elist(maxerr);
% ***JUMP OUT OF DO-LOOP
if( level(maxerr)+1<=levmax )
gt(1)=1;
break;
end;
nrmax = fix(nrmax + 1);
end;
if(gt(1)==1)
continue;
end;
end;
%
% PERFORM EXTRAPOLATION.
%
numrl2 = fix(numrl2 + 1);
rlist2(numrl2) = area;
if( numrl2>2 )
[numrl2,rlist2,reseps,abseps,res3la,nres]=dqelg(numrl2,rlist2,reseps,abseps,res3la,nres);
ktmin = fix(ktmin + 1);
if( ktmin>5 && abserr<0.1d-02.*errsum )
ier = 5;
end;
if( abseps<abserr )
ktmin = 0;
abserr = abseps;
resultf = reseps;
correc = erlarg;
ertest = max(epsabs,epsrel.*abs(reseps));
% ***JUMP OUT OF DO-LOOP
if( abserr<ertest )
break;
end;
end;
%
% PREPARE BISECTION OF THE SMALLEST INTERVAL.
%
if( numrl2==1 )
noext = true;
end;
if( ier>=5 )
break;
end;
end;
maxerr = fix(iord(1));
errmax = elist(maxerr);
nrmax = 1;
extrap = false;
levmax = fix(levmax + 1);
erlarg = errsum;
end;
end;
%
% SET THE FINAL RESULT.
% ---------------------
%
%
while(gt(2)==0);
if( abserr~=oflow )
if((ier+ierro)~=0 )
if( ierro==3 )
abserr = abserr + correc;
end;
if( ier==0 )
ier = 3;
end;
if( resultf==0.0d+00 || area==0.0d+00 )
if( abserr>errsum )
break;
end;
if( area==0.0d+00 )
gt(3)=1;
break;
end;
elseif( abserr./abs(resultf)>errsum./abs(area) ) ;
break;
end;
end;
%
% TEST ON DIVERGENCE.
%
if( ksgn~=(-1) || max(abs(resultf),abs(area))>defabs.*0.1d-01 )
if( 0.1d-01>(resultf./area) ||(resultf./area)>0.1d+03 ||errsum>abs(area) )
ier = 6;
end;
end;
gt(3)=1;
break;
end;
break;
end;
%
% COMPUTE GLOBAL INTEGRAL SUM.
%
if(gt(3)==0)
resultf = 0.0d+00;
for k = 1 : last;
resultf = resultf + rlist(k);
end; k = fix(last+1);
abserr = errsum;
end;
if( ier>2 )
ier = fix(ier - 1);
end;
resultf = resultf.*signmlv;
alist_shape=zeros(alist_shape);alist_shape(:)=alist(1:numel(alist_shape));alist=alist_shape;
blist_shape=zeros(blist_shape);blist_shape(:)=blist(1:numel(blist_shape));blist=blist_shape;
elist_shape=zeros(elist_shape);elist_shape(:)=elist(1:numel(elist_shape));elist=elist_shape;
iord_shape=zeros(iord_shape);iord_shape(:)=iord(1:numel(iord_shape));iord=iord_shape;
level_shape=zeros(level_shape);level_shape(:)=level(1:numel(level_shape));level=level_shape;
ndin_shape=zeros(ndin_shape);ndin_shape(:)=ndin(1:numel(ndin_shape));ndin=ndin_shape;
points_shape=zeros(points_shape);points_shape(:)=points(1:numel(points_shape));points=points_shape;
pts_shape=zeros(pts_shape);pts_shape(:)=pts(1:numel(pts_shape));pts=pts_shape;
rlist_shape=zeros(rlist_shape);rlist_shape(:)=rlist(1:numel(rlist_shape));rlist=rlist_shape;
end %subroutine dqagpe
%DECK DQAGP
|
|
