#To unpack, sh this file cat > READ.ME <<'CUT HERE............' *************************************************************************** * All the software contained in this library is protected by copyright. * * Permission to use, copy, modify, and distribute this software for any * * purpose without fee is hereby granted, provided that this entire notice * * is included in all copies of any software which is or includes a copy * * or modification of this software and in all copies of the supporting * * documentation for such software. * *************************************************************************** * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED * * WARRANTY. IN NO EVENT, NEITHER THE AUTHORS, NOR THE PUBLISHER, NOR ANY * * MEMBER OF THE EDITORIAL BOARD OF THE JOURNAL "NUMERICAL ALGORITHMS", * * NOR ITS EDITOR-IN-CHIEF, BE LIABLE FOR ANY ERROR IN THE SOFTWARE, ANY * * MISUSE OF IT OR ANY DAMAGE ARISING OUT OF ITS USE. THE ENTIRE RISK OF * * USING THE SOFTWARE LIES WITH THE PARTY DOING SO. * *************************************************************************** * ANY USE OF THE SOFTWARE CONSTITUTES ACCEPTANCE OF THE TERMS OF THE * * ABOVE STATEMENT. * *************************************************************************** AUTHORS: T.O. ESPELID AND K.J OVERHOLT UNIVERSITY OF BERGEN - NORWAY REFERENCE: - DQAINF: AN ALGORITHM FOR AUTOMATIC INTEGRATION OF INFINITE OSCILLATING TAILS NUMERICAL ALGORITHMS, 8(1994), PP. 83-101. SOFTWARE REVISION DATE: MARCH 10, 1994 REMARKS: This disk contains a makefile and FORTRAN-77 source programs for integrating infinite oscillating tails: 1) one makefile 2) single precision programs: The driver program for the example: STESTI and seven subroutines SQAINF,SCHINF,SGAINF,SADINF,SGAINT,STRINT,SEXINF + R1MACH (netlib/blas) Warning: the machine constants are set as for IEEE arithmetic machines in r1mach. If you are running a different type of machine then comment these lines out and activate the correct constants. 3)double precision programs: The driver program for the example: DTESTI and seven subroutines DQAINF,DCHINF,DGAINF,DADINF,DGAINT,DTRINT,DEXINF + D1MACH (netlib/blas) Warning: the machine constants are set as for IEEE arithmetic machines in d1mach. If you are running a different type of machine then comment these lines out and activate the correct constants. 4)then run either make stesti or make dtesti Next the single precision or double precision example may be run. Good luck! CUT HERE............ cat > makefile <<'CUT HERE............' # # dtesti.f is a driver designed to illustrate the use of # the double precision code on one example. # Note: d1mach.f (from netlib/blas) is included. # # # stesti.f is a driver designed to illustrate the use of # the single precision code on one example. # Note: r1mach.f (from netlib/blas) is included. # dtesti: dtesti.o dqainf.o dadinf.o dchinf.o dexinf.o\ dtrint.o dgaint.o dgainf.o d1mach.o f77 -o dtesti dtesti.o dqainf.o dadinf.o dchinf.o \ dtrint.o dgaint.o dgainf.o d1mach.o dexinf.o stesti: stesti.o sqainf.o sadinf.o schinf.o sexinf.o\ strint.o sgaint.o sgainf.o r1mach.o f77 -o stesti stesti.o sqainf.o sadinf.o schinf.o \ strint.o sgaint.o sgainf.o r1mach.o sexinf.o dtesti.o: dtesti.f f77 -u -c dtesti.f dqainf.o: dqainf.f f77 -u -c dqainf.f dadinf.o: dadinf.f f77 -u -c dadinf.f dchinf.o: dchinf.f f77 -u -c dchinf.f dtrint.o: dtrint.f f77 -u -c dtrint.f dgaint.o: dgaint.f f77 -u -c dgaint.f dgainf.o: dgainf.f f77 -u -c dgainf.f d1mach.o: d1mach.f f77 -u -c d1mach.f dexinf.o: dexinf.f f77 -u -c dexinf.f stesti.o: stesti.f f77 -u -c stesti.f sqainf.o: sqainf.f f77 -u -c sqainf.f sadinf.o: sadinf.f f77 -u -c sadinf.f schinf.o: schinf.f f77 -u -c schinf.f strint.o: strint.f f77 -u -c strint.f sgaint.o: sgaint.f f77 -u -c sgaint.f sgainf.o: sgainf.f f77 -u -c sgainf.f s1mach.o: s1mach.f f77 -u -c s1mach.f sexinf.o: sexinf.f f77 -u -c sexinf.f CUT HERE............ cat > dtesti.f <<'CUT HERE............' C C The following program will approximate three integrals over the C interval (0,infinity) and print out the values of the estimated C integrals, the estimated errors, the number of function evaluations C and IFAIL. C PROGRAM DTESTI C C This demo-program has calls to the main driver: DQAINF C Through this main driver the rest of the modules are activated: C C level 1: DQAINF C C level 2: DCHINF and DADINF C C level 3: DGAINF, DGAINT, DEXINF, DTRINT C C level 4: D1MACH (from blas in netlib) and FUNSUB (user specified) C EXTERNAL FUNSUB INTEGER NUMFUN,NW,MINPTS,MAXPTS,KEY,RESTAR INTEGER WRKSUB,EMAX,NUMNUL,NUM C C For the user's convenience we give a value to several parameters C through the parameter statement. C The following two parameters the user should not change, C unless the user changes the code appropriately: C PARAMETER (NUM=21,NUMNUL=16) C C The following three parameters the user may choose and these have C to remain unchanged through the entire run: C NUMFUN: the number of functions you want to integrate; C WRKSUB: the maximum number of subregions allowed; C EMAX: the maximum number of extrapolation steps; C PARAMETER (NUMFUN=3,WRKSUB=100,EMAX=25) C C Based on the previous 5 parameters we can now compute the C size of the main working space: NW. C PARAMETER (NW=1+5*NUMFUN+4*WRKSUB*NUMFUN+3*WRKSUB+(EMAX+1)**2+ +(NUMNUL+1+NUM)*NUMFUN) C C Next we choose the method to be used C Three options: KEY = 0(Euler), 1(Euler modified), 2(Overholt's C P-order 2) C PARAMETER (KEY=2) C C We need to choose WRKSUB such that: C WRKSUB >= MAXPTS/NUM (We have NUM=21 fixed in this implementation) C This simply means that we have enough space to achieve MAXPTS C function evaluations. By choosing WRKSUB bigger than C this lower bound we may later change MAXPTS and RESTAR and restart C the computations from the point where the last run stopped. C The reason for stopping may have been that MAXPTS was too small C to achieve the requested error. C With our choice WRKSUB = 100 then any value of MAXPTS C between MINPTS and 2100 will be accepted by the code. Choosing C MAXPTS = 1000 then allows us to restart with a greater value C if necessary. C Note: For a given problem it may happen that the code stops after C reaching say MAXPTS=2100 function evaluations, having used less than C WRKSUB = 100 subregions. The reason for this is that if a given C region needs to be further subdivided then the old region will not C be stored while the function values computed over that region C will be counted. C The example gave the following output on C a SUN Sparc station 10: C C A= 0. B= 3.0000000000000 C GAMMA= 1.0000000000000 PERIOD= 6.2831853071796 C MAXPTS= 150 WRKSUB= 100 NW= 2306 C KEY= 2 EPSREL= 1.0000000000000D-09 C C Result and error: problem no. 1, 2 and 3: C 0.12431084921897E+01 0.71825718819736E+00 0.48923485452174E+00 C 0.59475024208688E-05 0.35781998629104E-06 0.94496089623798E-07 C NEVAL= 147 IFAIL= 1 C C New attempt using the restart feature, with MAXPTS= 300 C Result and error: problem no. 1, 2 and 3: C 0.12431084850988E+01 0.71825716313760E+00 0.48923480269063E+00 C 0.38596743952420E-10 0.44266003080563E-10 0.81963473649869E-10 C NEVAL= 273 IFAIL= 0 C DOUBLE PRECISION A(WRKSUB),B(WRKSUB),EPSABS,EPSREL,GAMMA,PERIOD, +RESULT(NUMFUN),ABSERR(NUMFUN),WORK(NW),PI INTEGER NEVAL,IFAIL,IWORK(3*WRKSUB) LOGICAL DONE(NUMFUN) PI=2*ASIN(1.D0) C C We specify the period and the decay parameter GAMMA>0. If C the value of GAMMA < 0 then the code will attempt to estimate C GAMMA in the case KEY = 2. GAMMA is not used by the other two C methods. C PERIOD=2*PI GAMMA=1.0D0 C C Next we give the left endpoint of integration and the point b C from where we expect the theory to be valid. Advice: Since b C is not uniquely defined please do some experiments with C different values of b to gain experience. C The effect one achieves by increasing b depends on the value of KEY. C A(1)=0 B(1)=3 C C This is a first attempt: set RESTAR, give MINPTS, MAXPTS and the C absolute and relative error requests. C RESTAR=0 MINPTS=42 MAXPTS=150 EPSABS=0 EPSREL=1.D-9 C C Note: these four parameters and RESTAR may be changed in a restart run C WRITE (*,*) 'A=',A(1),' B=',B(1) WRITE (*,*) 'GAMMA=',GAMMA,' PERIOD=',PERIOD WRITE (*,*) 'MAXPTS=',MAXPTS,' WRKSUB=',WRKSUB,' NW=',NW WRITE (*,*) 'KEY=',KEY,' EPSREL=',EPSREL CALL DQAINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,WRKSUB,EMAX,MINPTS, +MAXPTS,EPSABS,EPSREL,NW,KEY,RESTAR,RESULT,ABSERR,NEVAL,IFAIL,DONE, +WORK,IWORK) WRITE (*,*) WRITE (*,*) 'Result and error: problem no. 1, 2 and 3:' WRITE (*,10) RESULT WRITE (*,10) ABSERR WRITE (*,*) 'NEVAL=',NEVAL,' IFAIL=',IFAIL C C Test if this attempt was successful C IF (IFAIL.EQ.1) THEN C C We make a new attempt by increasing MAXPTS and using the C restarting feature by setting RESTAR = 1. We keep MINPTS C EPSABS and EPSREL unchanged in this example.. C RESTAR=1 MAXPTS=300 WRITE (*,*) WRITE (*,*) 'New attempt using the restart feature, with MAXPTS +=',MAXPTS CALL DQAINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,WRKSUB,EMAX,MINPTS, + MAXPTS,EPSABS,EPSREL,NW,KEY,RESTAR,RESULT,ABSERR,NEVAL,IFAIL, + DONE,WORK,IWORK) WRITE (*,*) 'Result and error: problem no. 1, 2 and 3:' WRITE (*,10) RESULT WRITE (*,10) ABSERR WRITE (*,*) 'NEVAL=',NEVAL,' IFAIL=',IFAIL END IF 10 FORMAT (' ',3E22.14) END SUBROUTINE FUNSUB (X,NUMFUN,NP,FUNVLS) INTEGER NUMFUN,I,NP,J DOUBLE PRECISION X(NP),FUNVLS(NP,1),Y C C This subroutine must be provided by the user for computing C all components of the integrand at the given C evaluation points. C It must have parameters (X,NUMFUN,NP,FUNVLS) C Input parameters: C X Real array of dimension NP (= 21) defining the C x-coordinates of the evaluation points. C NUMFUN Integer that defines the number of C components of the integral I. C NP Integer that gives the number of evaluation points C in the quadrature rule used (Gauss, 21 point rule). C Output parameter: C FUNVLS Real array of dimension (NP,NUMFUN) that defines C the function values in the 21 evaluation points C for the NUMFUN components of the integrand. C C Note that we may save computer time when integrating C several functions simultaneously over the same interval C if we take advantage of the functions' similarities. C DO 10 I=1,NP Y=(2*SIN(X(I))+X(I)*COS(X(I))) DO 10 J=1,NUMFUN FUNVLS(I,J)=Y/(J+X(I)**2) 10 CONTINUE RETURN END CUT HERE............ cat > dqainf.f <<'CUT HERE............' SUBROUTINE DQAINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,WRKSUB,EMAX, +MINPTS,MAXPTS,EPSABS,EPSREL,NW,KEY,RESTAR,RESULT,ABSERR,NEVAL, +IFAIL,DONE,WORK,IWORK) C***BEGIN PROLOGUE DQAINF C***DATE WRITTEN 930515 (YYMMDD) C***REVISION DATE 940310 (YYMMDD) C***CATEGORY NO. H2B2A1/H2B2A2 C***AUTHORS C Terje O. Espelid, Department of Informatics, C University of Bergen, Hoyteknologisenteret. C N-5020 Bergen, Norway C Email.. terje@ii.uib.no C C Kjell J. Overholt, Department of Informatics, C University of Bergen, Hoyteknologisenteret. C N-5020 Bergen, Norway C Email.. kjell@ii.uib.no C C***PURPOSE To compute several one-dimensional integrals over C an infinite region. All of the integrands are assumed to be C oscillatory, with the same asymptotic periodic behavior. C This routine is the driving routine for the integration C routine DADINF. C C***KEYWORDS Quadrature, one dimensional, oscillatory integrands, C infinite region, extrapolation, globally adaptive, C periodic functions. C C***DESCRIPTION C C The routine calculates an approximation to a given vector of C infinite integrals C C I (F ,F ,...,F ) DX, C 1 2 NUMFUN C C where the interval of integration is [A,infinity), where A is C supplied by the user, and with C B as a possible subdivision point also supplied by the user: C A <= B. We assume that all functions have the same oscillatory C behavior for values of X => B. C C Here F = F (X), J = 1,2,...,NUMFUN. C J J C C A globally adaptive strategy, combined with extrapolation, C is applied in order to compute approximations RESULT(J), C for each component J of I, the following claim for accuracy: C C ABS(I(J)-RESULT(J)).LE.MAX(EPSABS,EPSREL*ABS(I(J))) C C DQAINF is a driver for the integration routine DADINF. C C DADINF either (1) repeatedly subdivides the interval with greatest C estimated errors; then estimates the integrals and the errors over C the new sub-intervals and finally updates the extrapolation tableau, C until the error request is met or MAXPTS function C evaluations have been used, or (2) decides instead to compute C a new term: the integral over the next half-period and then C perform a new extrapolation step. C C Only one basic integration rule (21 points) is offered: Gauss rule. C C Three different extrapolation methods are offered through the C parameter KEY:: Euler's method, a modification of Eulers method and C Overholt's P-order 2-method. C C If the values of the input parameters EPSABS or EPSREL are selected C great enough, the chosen integration rule is applied over the two C first intervals [A,B] and [B,B+PERIOD/2] for each C function in the vector and then one extrapolation step is performed C to give approximations RESULT(J), J = 1, 2, ..., NUMFUN. C No further subdivisions and no C more new terms will then be computed. C C When DQAINF computes estimates to a vector of integrals, C then all components of the vector are given the same treatment. C That is, I(F ) and I(F ) for J not equal to K, are estimated with C J K C the same subdivision of the original interval of integration C and the same extrapolation scheme is applied. C For integrals with enough similarity, we may save time by applying C DQAINF to all integrands in one call. For integrals that varies C continuously as functions of some parameter, the estimates C produced by DQAINF will also vary continuously when the same C subdivision is applied to all components. This will generally not C be the case when the different components are given separate C treatment. C C On the other hand this feature should be used with caution when the C different components of the integrals require clearly different C subdivisions. C C ON ENTRY C C NUMFUN Integer. C Number of components of the integral. C FUNSUB Externally declared subroutine for computing C all components of the integrand at the given C evaluation point. C It must have parameters (X,NUMFUN,NP,FUNVLS) C Input parameters: C X Real array of dimension 21 defining the C x-coordinates of the evaluation points. C NUMFUN Integer that defines the number of C components of the integral I. C NP Integer that gives the number of evaluation points C Output parameter: C FUNVLS Real array of dimension (21,NUMFUN) that defines C the function values in the 21 evaluation points C for the NUMFUN components of the integrand. C PERIOD Real. C All functions are assumed to have the same asymptotic C period > 0. C GAMMA Real. C All functions are assumed to decay as c/x**GAMMA, C when x >> 1 and we go in steps of length PERIOD, C (GAMMA > 0). C A Real array of dimension WRKSUB. C B Real array of dimension WRKSUB. C A(1) and B(1) are the left endpoint and right endpoint C of the first interval, A(1) <= B(1). B(1) is chosen by the C user and it is assumed that the integrand oscillates for C X >= B(1). Asymptotic period is PERIOD. Thus oscillations C are expected to be observed for points of distance PERIOD/2 C A(2), B(2), etc. are defined by the code, and as a warning: C the code may change the value of A(1) and B(1). C WRKSUB Integer. C Defines the length of the arrays A and B. (And thus the max C number of subregions allowed.) C EMAX Integer. C The maximum number of extrapolation steps. C MINPTS Integer. C Minimum number of function evaluations. C MAXPTS Integer. C Maximum number of function evaluations. C The number of function evaluations over each sub-interval C is 21. C EPSABS Real. C Requested absolute error. C EPSREL Real. C Requested relative error. C NW Integer. C Defines the length of the working array WORK. C We let C KEY Integer C Choice of extrapolation method: C KEY = 0 : Euler's method. C KEY = 1 : Modified Euler. C KEY = 2 : Overholt's polynomial order 2-method. C This last method is the only one that needs C the value of GAMMA. C RESTAR Integer. C If RESTAR = 0, this is the first attempt to compute C the integral. C If RESTAR = 1, C then we restart a previous attempt. C In this case the only parameters for DQAINF that may C be changed (with respect to the previous call of DQAINF) C are MINPTS, MAXPTS, EPSABS, EPSREL and RESTAR. C C ON RETURN C C RESULT Real array of dimension NUMFUN. C Approximations to all components of the integral. C ABSERR Real array of dimension NUMFUN. C Estimates of absolute errors. C NEVAL Integer. C Number of function evaluations used by DQAINF. C IFAIL Integer. C IFAIL = 0 for normal exit, when C C ABSERR(J) <= EPSABS or C ABSERR(J) <= ABS(RESULT(J))*EPSREL with MAXPTS or less C function evaluations for all values of J, C from 1 to NUMFUN. C C IFAIL = +1 if MAXPTS was too small C to obtain the required accuracy. In this case DQAINF C returns values of RESULT with estimated absolute C errors ABSERR. C IFAIL = -1 if EMAX was too small C to obtain the required accuracy. In this case DQAINF C returns values of RESULT with estimated absolute C errors ABSERR. C IFAIL = 2 if NUMFUN is less than 1. C IFAIL = 3 if B(1) < A(1). C IFAIL = 4 if unlegal PERIOD. C IFAIL = 5 if MAXPTS is less than MINPTS or MINPTS < 21. C IFAIL = 6 if EPSABS <= 0 and EPSREL <= 0. C IFAIL = 7 if WRKSUB is too small. C IFAIL = 8 if unlegal value of EMAX. C IFAIL = 9 if unlegal RESTAR. C IFAIL = 10 if unlegal value of KEY. C IFAIL = 11 if NW is less than LIMIT (See the code DCHINF). C IFAIL = 12 if either PERIOD is wrong or B(1) is too small. C IFAIL = 13 if unable to estimate GAMMA (In case KEY=2 only) C A Real array of dimension WRKSUB. C B Real array of dimension WRKSUB. C On exit A and B contains the endpoints of the intervals C used to produce the approximations to the integrals. C DONE Logical array of dimension NUMFUN. C On exit : DONE(I)=.TRUE. if function number I has been C integrated to the specified precision, else DONE(I)=.FALSE. C (Note that IFAIL = 1 just tells you C that something is wrong, however most of the integrals in C the vector may have been computed to the specified precisio C WORK Real array of dimension NW. C Used as working storage. C WORK(NW) = NSUB, the number of sub-intervals in the data C structure. C Let NW >= 1+5*NUMFUN+4*WRKSUB*NUMFUN C +3*WRKSUB+(EMAX+1)**2 + (NUMNUL+3+NUM)*NUMFUN) C Then C WORK(1),...,WORK(NUMFUN*WRKSUB) contain the estimated C components of the integrals over the sub-intervals. C WORK(NUMFUN*WRKSUB+1),...,WORK(2*NUMFUN*WRKSUB) contain C the estimated errors over the sub-intervals. C WORK(2*NUMFUN*WRKSUB+1),...,WORK(2*NUMFUN* C WRKSUB+WRKSUB) contain the greatest errors C in each sub-interval. C WORK(2*WRKSUB*NUMFUN+WRKSUB+1),...,WORK(2*WRKSUB*NUMFUN C +2*WRKSUB) contain the left sub-division point in each C sub-interval. C WORK(2*WRKSUB*NUMFUN+2*WRKSUB+1),...,WORK(2*WRKSUB*NUMFUN C +3*WRKSUB) contain the right sub-division point in each C sub-interval. C WORK((2*NUMFUN+3)*WRKSUB+1),...,WORK(NW) is used as C temporary storage in DADINF. C IWORK Integer array of dimension 3*WRKSUB. C Used as working storage. C C***LONG DESCRIPTION C C The information for each interval is contained in the data C structures A,B, WORK and IWORK. A and B contain the endpoints of C the intervals. When passed on to DADINF, WORK is split into four C arrays VALUES, ERRORS, GREATE and WORK2. VALUES contains the C estimated values of the integrals. ERRORS contains the estimated C errors. GREATE contains the greatest estimated error for each C interval. The data structures for the intervals are in DADINF C organized as a heap, and the size of GREATE(I) defines the C position of interval I in the heap. The heap is maintained by the C program DTRINF and we use a partially ordered list of pointers, C saved in IWORK. When passed on to DADINF, IWORK is split into three C arrays WORST, LIST and UPOINT. LIST is a partially ordered list C such that GREATE(LIST(1)) is the maximum worst case error estimate C for all sub-intervals in our data-structure. In WORST the index to C the integral with the greatest error-estimate is kept. UPOINT is C an array of pointers to where in the U-sequence a region belongs. C This information is used when updating the corresponding U-term C after a subdivision. C C The subroutine for estimating the integral and the error over each C sub-interval, DGAINT, uses WORK2 as a work array. C C***REFERENCES C C T.O.Espelid and K.J.Overholt, DQAINF: An Algorithm for Automatic C Integration of Infinite Oscillating Tails, submitted for publ. 1993. C C***ROUTINES CALLED DCHINF,DADINF C***END PROLOGUE DQAINF C C Parameters C C C Global variables. C EXTERNAL FUNSUB INTEGER WRKSUB INTEGER NUMFUN,MINPTS,MAXPTS,NW,RESTAR,NEVAL,IFAIL,IWORK(3*WRKSUB) +,KEY,EMAX DOUBLE PRECISION A(WRKSUB),B(WRKSUB),EPSABS,EPSREL,RESULT(NUMFUN), +ABSERR(NUMFUN),WORK(NW),PERIOD,GAMMA LOGICAL DONE(NUMFUN) C C Local variables. C INTEGER NSUB,LENW,NUM,NUMNUL INTEGER I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12,I13,I14 PARAMETER (NUM=21,NUMNUL=16) C C***FIRST EXECUTABLE STATEMENT DQAINF C C Check the input parameters. C CALL DCHINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,MINPTS,MAXPTS,EPSABS, +EPSREL,WRKSUB,NW,EMAX,KEY,RESTAR,NEVAL,IFAIL) IF (IFAIL.NE.0) THEN RETURN END IF C C Split up the work space. C I1=1 I2=I1+WRKSUB*NUMFUN I3=I2+WRKSUB*NUMFUN I4=I3+WRKSUB I5=I4+WRKSUB I6=I5+WRKSUB I7=I6+WRKSUB*NUMFUN I8=I7+WRKSUB*NUMFUN I9=I8+NUMFUN I10=I9+NUMFUN I11=I10+NUMFUN I12=I11+(EMAX+1)*(EMAX+1) I13=I12+NUMFUN I14=I13+NUMFUN C C On restart runs the number of sub-intervals from the C previous call is assigned to NSUB. C IF (RESTAR.EQ.1) THEN NSUB=WORK(NW) END IF C C Compute the size of the temporary work space needed in DADINF. C LENW=(NUMNUL+1+NUM)*NUMFUN CALL DADINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,WRKSUB,EMAX,MINPTS, +MAXPTS,EPSABS,EPSREL,RESTAR,LENW,KEY,RESULT,ABSERR,NEVAL,NSUB, +IFAIL,DONE,WORK(I1),WORK(I2),WORK(I3),WORK(I4),WORK(I5),WORK(I6), +WORK(I7),WORK(I8),WORK(I9),WORK(I10),WORK(I11),WORK(I12),WORK(I13) +,WORK(I14),IWORK(1),IWORK(1+WRKSUB),IWORK(1+2*WRKSUB)) WORK(NW)=NSUB RETURN C C***END DQAINF C END CUT HERE............ cat > dadinf.f <<'CUT HERE............' SUBROUTINE DADINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,WRKSUB,EMAX, +MINPTS,MAXPTS,EPSABS,EPSREL,RESTAR,LENW,KEY,RESULT,ABSERR,NEVAL, +NSUB,IFAIL,DONE,VALUES,ERRORS,GREATE,C,D,U,E,T,CT,EXTERR,BETA, +UNEW,UERR,WORK2,WORST,LIST,UPOINT) C***BEGIN PROLOGUE DADINF C***KEYWORDS Quasi-linear extrapolation, infinite integration, C oscillating functions, adaptive strategy. C***PURPOSE To compute better estimates to a vector of approximations C to integrals of oscillating functions over an infinite C interval. C C***LAST MODIFICATION 94-03-10 C C***REFER TO DQAINF C C***DESCRIPTION Computation of integrals over an infinite interval C by combining extrapolation and adaptive strategies. C C DADINF uses extrapolation on a sum of integrals over subintervals C Each subinterval has length P = PERIOD/2 and the function involved i C assumed to be asymptotically periodic with period PERIOD. C This routine will EITHER take another extrapolation step C in order to reduce the pure extrapolation error, OR C subdivide the interval in the collection with greatest estimated C error. In both cases estimate(s) of the integral(s) and the error(s) C to the new sub-interval(s) are computed; this process continues C until the error request is met or MAXPTS function evaluations have C been used. C C The 21 point rule is: Gauss. C C ON ENTRY C C NUMFUN Integer. C The number of components of the integral. C FUNSUB Externally declared subroutine for computing C all components of the integrand at all NP C evaluation points. C It must have parameters (X,NUMFUN,FUNVLS) C Input parameters: C X The x-coordinates of the NP evaluation points. C NUMFUN Integer that defines the number of C components of I. C NP Integer that gives the number of evaluation points C Output parameter: C FUNVLS Real array of dimension NUMFUN C that defines NUMFUN components of the integrand. C C PERIOD Real. C All functions are assumed to have the same asymptotic C period. C GAMMA Real. C All functions are assumed to decay as c/x**GAMMA, C when x >> 1 and we go in steps of length PERIOD, C (GAMMA > 0). C A Real array of dimension WRKSUB. C B Real array of dimension WRKSUB. C A(1) is the left endpoint of integration. C B(1) >= A(1) is a user specified subdivision point. C We assume that the assumptions on which this code is C based hold for x >= B(1). We will not extrapolate based C on the function values between A(1) and B(1). C WRKSUB Integer. C Defines the length of the arrays A and B . C EMAX Integer. C The maximum number of extrapolation steps. C MINPTS Integer. C Minimum number of function evaluations. C MAXPTS C Maximum number of function evaluations. C EPSABS Real. C Requested absolute error. C EPSREL Real. C Requested relative error. C C RESTAR Integer. C If RESTAR = 0, this is the first attempt to compute C the integral. C If RESTAR = 1, C then we restart a previous attempt. C In this case the only parameters for DQAINF that may C be changed (with respect to the previous call of DQAINF) C are MINPTS, MAXPTS, EPSABS, EPSREL and RESTAR. C LENW Integer. C Length of the workspace WORK2. C KEY Integer C Choice of extrapolation method: C KEY = 0 : Euler's method. C KEY = 1 : Modified Euler. C KEY = 2 : Overholt's polynomial order 2-method. C This last method is the only one that needs C GAMMA. C C ON RETURN C C RESULT Real array of dimension NUMFUN. C Approximations to all components of the integral. C ABSERR Real array of dimension NUMFUN. C Estimates of absolute errors. C NEVAL Integer. C The number of function evaluations. C NSUB Integer. C The number of intervals in the data structure. C A Real array of dimension WRKSUB. C B Real array of dimension WRKSUB. C On exit A and B contains the endpoints of the intervals C used to produce the approximations to the integrals. C Warning: A(1) and B(1) may have changed due to C adaptive subdivision of the first interval. C IFAIL Integer. C IFAIL = 0 for normal exit. C IFAIL = +1 if MAXPTS was too small for DADINF C to obtain the required accuracy. In this case C DADINF returns values of RESULT with estimated C absolute errors ABSERR. C IFAIL = -1 if EMAX was too small for DADINF C to obtain the required accuracy. In this case C DADINF returns values of RESULT with estimated C absolute errors ABSERR. C C DONE Logical array of dimension NUMFUN. C On exit : DONE(I)=.TRUE. if function number I has been C integrated to the specified precision, else DONE(I)=.FALSE. C C VALUES Real array of dimension (NUMFUN,WRKSUB). C The estimated components of the integrals over the C sub-intervals. C ERRORS Real array of dimension (NUMFUN,WRKSUB). C The estimated errors over the sub-intervals. C GREATE Real array of dimension WRKSUB. C The greatest error in each sub-interval. C C Real array of dimension WRKSUB. C The left sub-division points of all intervals in the heap. C D Real array of dimension WRKSUB. C The right sub-division points of all intervals in the heap. C U Real array of dimension (NUMFUN,WRKSUB) C contains the terms in series. C E Real array of dimension (NUMFUN,WRKSUB) C contains the estimated errors in each U-term. C T Real array of dimension NUMFUN. C Dummy parameter. C CT Real array of dimension NUMFUN. C Dummy parameter. C EXTERR Real array of dimension NUMFUN. C These errors are associated with the region (LEFT,infinity) C and they are the pure extrapolation errors. C BETA Real array of dimension (0:EMAX,0:EMAX). C Dummy parameter. C UNEW Real array of dimension NUMFUN. C This gives the next terms in the series (new extrapolation C step) else it is the correction to the U-values (updating). C UERR Real array of dimension NUMFUN. C The estimated errors of all U-terms in the series. C WORK2 Real array of dimension LENW. C Work array used in DGAINT. C WORST Integer array of dimension WRKSUB. C Index to the greatest error in each sub-interval. C LIST Integer array used in DTRINT of dimension WRKSUB. C Is a partially sorted list, where LIST(1) is the top C element in a heap of sub-intervals. C UPOINT Integer array of dimension WRKSUB. C Is an array of pointers to where in the U-sequence C a region belongs. This information is used when updating C the corresponding U-term after a subdivision. C C***ROUTINES CALLED DTRINT, DGAINT, DEXINF, D1MACH. C***END PROLOGUE DADINF C C Global variables. C EXTERNAL FUNSUB,D1MACH INTEGER NUMFUN,LENW,RESTAR,NEVAL,WRKSUB,NSUB,IFAIL, +LIST(WRKSUB),KEY,WORST(WRKSUB),UPOINT(WRKSUB),EMAX, +MAXPTS,MINPTS DOUBLE PRECISION A(WRKSUB),B(WRKSUB),EPSABS,EPSREL,RESULT(NUMFUN), +ABSERR(NUMFUN),VALUES(NUMFUN,WRKSUB),PERIOD,ERRORS(NUMFUN,WRKSUB), +GREATE(WRKSUB),WORK2(LENW),C(WRKSUB),D(WRKSUB),U(NUMFUN,WRKSUB), +E(NUMFUN,WRKSUB),GAMMA,EXTERR(NUMFUN),T(NUMFUN),CT(NUMFUN), +BETA(0:EMAX,0:EMAX),UNEW(NUMFUN),UERR(NUMFUN),D1MACH LOGICAL DONE(NUMFUN) C C Local variables. C C SBRGNS is the number of stored sub-intervals. C NDIV The number of sub-intervals to be divided in each main step. C POINTR Pointer to the position in the data structure where C the new sub-intervals are to be stored. C G is the homogeneous coordinates of the integration rule. C W is the weights of the integration rule and the null rules. C TOP is a pointer to the top element in the heap of sub-intervals. C AOLD Left-endpoint of the interval to be processed. C BOLD Right-endpoint of the interval to be processed. C FLAG Signals small interval. C MORE Signals more work to do. C P Half a period. C NUMU The number of U-terms. C NUM The number of points in the basic rule. C NUMNUL The number nullrules used. C UPDATE Index to U-term to be updated. If UPDATE < 0, then this means C a new extrapolation step. C VACANT Index to vacant position in the data-structure. C EPMACH Machine epsilon C CC Ratio: the user specified subdivision point B(1)/PERIOD. INTEGER I,J,K,JJ,NUMINT,ONE,NUMU,VACANT,UPDATE INTEGER SBRGNS,I1,I2,I3,L1,L2,L3 INTEGER NDIV,POINTR,INDEX,TOP,FLAG,NUM,NUMNUL DOUBLE PRECISION AOLD,BOLD,GREAT,P,LEFT,MAXEER,EPMACH,CC LOGICAL MORE SAVE MAXEER,NUMU,LEFT,P,CC PARAMETER (NUMINT=2,NUM=21,NUMNUL=16) C C***FIRST EXECUTABLE STATEMENT DADINF C C Define machine epsilon. C EPMACH=D1MACH(4) C C Set some pointers for WORK2. C L1=1 L2=L1+NUMNUL*NUMFUN L3=L2+NUM*NUMFUN C IF (RESTAR.EQ.1) THEN SBRGNS=NSUB GO TO 90 ELSE FLAG=0 CC=B(1)/PERIOD END IF C C Initiate SBRGNS, DONE, MORE, P, A and B. C MORE=.TRUE. SBRGNS=0 DO 10 J=1,NUMFUN DONE(J)=.FALSE. 10 CONTINUE P=PERIOD/2 A(2)=B(1) B(2)=A(2)+P LEFT=B(2) C C Initialize E and U C DO 20 J=1,NUMFUN DO 20 I=1,WRKSUB E(J,I)=0 20 U(J,I)=0 C C Apply the rule procedure over the two first intervals. C IF (A(1).EQ.B(1)) THEN DO 30 I=1,NUMFUN VALUES(I,1)=0 30 ERRORS(I,1)=0 GREATE(1)=0 WORST(1)=1 ONE=2 SBRGNS=1 UPOINT(1)=1 ELSE ONE=1 END IF DO 40 I=ONE,2 CALL DGAINT (A(I),B(I),NUMFUN,FUNSUB,DONE,MORE,EPMACH, + WORK2(L1),WORK2(L2),WORK2(L3),FLAG,VALUES(1,I), + ERRORS(1,I),GREATE(I),WORST(I),C(I),D(I)) NEVAL=NEVAL+NUM SBRGNS=SBRGNS+1 UPOINT(I)=I 40 CONTINUE C C Initialize RESULT, U, E and NUMU. C DO 60 I=1,2 DO 50 J=1,NUMFUN U(J,I)=VALUES(J,I) E(J,I)=ERRORS(J,I) 50 CONTINUE 60 CONTINUE NUMU=2 C C Store results in heap. C DO 80 I=1,NUMINT GREAT=0.D0 DO 70 J=1,NUMFUN IF (ERRORS(J,I).GT.GREAT) THEN GREAT=ERRORS(J,I) WORST(I)=J END IF 70 CONTINUE GREATE(I)=GREAT C(I)=A(I)+(B(I)-A(I))/3 D(I)=B(I)-(B(I)-A(I))/3 INDEX=I CALL DTRINT (2,INDEX,GREATE,LIST,I) 80 CONTINUE C C Extrapolate for the first time based on these 2 terms. C This will initialize the global error as well, and C MAXEER: the greatest extrapolation error. C Finally RESULT will be initialized. C UPDATE=-1 CALL DEXINF (NUMFUN,GAMMA,CC,KEY,NUMU-1,T,CT,UPDATE,EMAX,U,E, +RESULT,ABSERR,EXTERR,BETA,MAXEER) C C Check for termination. C IF (MORE.OR.(NEVAL.LT.MINPTS)) THEN GO TO 90 END IF IFAIL=0 GO TO 170 C C***End initialization. C C***Begin loop while the error is too great, C C First we check the ranking of the TOP element. C Then we determine if we will create a new term by integrating C from LEFT to LEFT + P and then extrapolate, or subdivide the TOP C element in the heap of subregions. C C 90 TOP=LIST(1) C C New term? C IF (GREATE(TOP).LT.(MAXEER)) THEN C C We want to extrapolate. Check if we C are allowed to take a new extrapolation step. C IF (NUMU-1.GT.EMAX) THEN IFAIL=-1 GO TO 170 ELSE VACANT=0 NDIV=1 END IF ELSE VACANT=TOP NDIV=3 END IF C C Check if NEVAL+NDIV*NUM is less than or equal to MAXPTS: C MAXPTS is the maximum number of function evaluations that are C allowed to be computed. C IF (NEVAL+NDIV*NUM.LE.MAXPTS) THEN C C When we pick an interval to divide in three, then one of the new C sub-intervals is stored in the original interval's position in the C data structure. C C Let POINTR point to the first free position in the data structure. C C Determine if this is a new extrapolation step or an update. C UPDATE will point to the element in the C U-series to be corrected or created. C IF (VACANT.EQ.0) THEN POINTR=SBRGNS+1 NUMU=NUMU+1 UPDATE=NUMU INDEX=POINTR A(INDEX)=LEFT LEFT=LEFT+P B(INDEX)=LEFT TOP=INDEX DO 100 J=1,NUMFUN UERR(J)=0 100 UNEW(J)=0 ELSE UPDATE=UPOINT(VACANT) C C Save the endpoints. C AOLD=A(TOP) BOLD=B(TOP) C C Remove the top element from the heap.(The value of K does not matter C POINTR=SBRGNS+2 CALL DTRINT (1,SBRGNS,GREATE,LIST,K) C C Compute the three new intervals. C I1=TOP I2=POINTR-1 I3=POINTR A(I1)=AOLD B(I1)=C(TOP) A(I2)=B(I1) B(I2)=D(TOP) A(I3)=B(I2) B(I3)=BOLD INDEX=VACANT DO 110 J=1,NUMFUN UERR(J)=-ERRORS(J,INDEX) 110 UNEW(J)=-VALUES(J,INDEX) END IF C C Apply basic rule. C DO 130 I=1,NDIV CALL DGAINT (A(INDEX),B(INDEX),NUMFUN,FUNSUB,DONE,MORE, + EPMACH,WORK2(L1),WORK2(L2),WORK2(L3),FLAG, + VALUES(1,INDEX),ERRORS(1,INDEX),GREATE(INDEX), + WORST(INDEX),C(INDEX),D(INDEX)) UPOINT(INDEX)=UPDATE C C Compute the U-term C DO 120 J=1,NUMFUN UNEW(J)=UNEW(J)+VALUES(J,INDEX) UERR(J)=UERR(J)+ERRORS(J,INDEX) 120 CONTINUE INDEX=SBRGNS+I+1 130 CONTINUE NEVAL=NEVAL+NDIV*NUM C C Compute the E and U terms (These may be new terms or terms that C have to be updated), C DO 140 J=1,NUMFUN U(J,UPDATE)=U(J,UPDATE)+UNEW(J) E(J,UPDATE)=E(J,UPDATE)+UERR(J) 140 CONTINUE C C Do the extrapolation and compute the global results and errors. C UPDATE is used to signal an extrapolation step. C IF (VACANT.EQ.0) THEN UPDATE=-1 END IF CALL DEXINF (NUMFUN,GAMMA,CC,KEY,NUMU-1,T,CT,UPDATE,EMAX,U,E, + RESULT,ABSERR,EXTERR,BETA,MAXEER) C C Check each integration problem and set DONE(J)=.TRUE. if they C are finished. MORE will signal if there is more to do. C MORE=.FALSE. DO 150 J=1,NUMFUN IF (ABSERR(J).LE.EPSREL*ABS(RESULT(J)).OR.ABSERR(J).LE. + EPSABS) THEN DONE(J)=.TRUE. ELSE MORE=.TRUE. DONE(J)=.FALSE. END IF 150 CONTINUE C C Store results in heap. C INDEX=TOP DO 160 I=1,NDIV JJ=SBRGNS+I CALL DTRINT (2,JJ,GREATE,LIST,INDEX) INDEX=SBRGNS+I+1 160 CONTINUE SBRGNS=SBRGNS+NDIV C C Check for termination. C IF (MORE.OR.(NEVAL.LT.MINPTS)) THEN GO TO 90 END IF C C Else we did not succeed with the C given value of MAXPTS. Note: We do not use IFAIL to signal C FLAG in the current implementation. C ELSE IFAIL=+1 END IF C 170 NSUB=SBRGNS RETURN C C***END DADINF C END CUT HERE............ cat > dchinf.f <<'CUT HERE............' SUBROUTINE DCHINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,MINPTS,MAXPTS, +EPSABS,EPSREL,WRKSUB,NW,EMAX,KEY,RESTAR,NEVAL,IFAIL) C***BEGIN PROLOGUE DCHINF C***REFER TO DQAINF C***PURPOSE DCHINF checks the validity of the input parameters to C DQAINF. C C***LAST MODIFICATION 93-05-05 C C***DESCRIPTION C DCHINF computes IFAIL as a function of the input C parameters to DQAINF, and checks the validity of the input C parameters to DQAINF. C C ON ENTRY C C NUMFUN Integer. C Number of components of the integral. C FUNSUB Externally declared subroutine for computing C all components of the integrand at the given C evaluation point. C It must have parameters (X,NUMFUN,FUNVLS) C Input parameters: C X Real array of dimension 21 defining the C x-coordinates of the evaluation points. C NUMFUN Integer that defines the number of C components of I. C NP Integer that gives the number of evaluation points C Output parameter: C FUNVLS Real array of dimension (21,NUMFUN) that defines C the function values in the 21 evaluation points C for the NUMFUN components of the integrand. C PERIOD Real. C All functions are assumed to have the same asymptotic C period. C GAMMA Real. C All functions are assumed to decay as c/x**GAMMA, C when x >> 1 and we go in steps of length PERIOD, C (GAMMA > 0). C A Real array of dimension WRKSUB. C B Real array of dimension WRKSUB. C A(1) and B(1) are the left endpoint and right endpoint C of the first interval, A(1) <= B(1). B(1) is chosen by the C user and it is assumed that the integrand oscillates for C X >= B(1). Asymptotic period is PERIOD. Thus oscillations C are expected to be observed for points of distance PERIOD/2 C MINPTS Integer. C Minimum number of function evaluations. C MAXPTS Integer. C Maximum number of function evaluations. C The number of function evaluations over each sub-interval C is NUM = 21. C EPSABS Real. C Requested absolute error. C EPSREL Real. C Requested relative error. C WRKSUB Integer. C Defines the length of the arrays A and B. (And thus the C maximum number of subregions allowed.) C NW Integer. C Defines the length of the working array WORK. C EMAX Integer. C The maximum number of extrapolation steps. (At least one C step!) C KEY Integer. C Choice of extrapolation method: C KEY = 0 : Euler's method. C KEY = 1 : Modified Euler. C KEY = 2 : Overholt's polynomial order 2-method. C This last method is the only one that needs C GAMMA. C RESTAR Integer. C If RESTAR = 0, this is the first attempt to compute C the integral. C If RESTAR = 1, C then we restart a previous attempt. C RESTAR is allowed to take these two values only. C ON RETURN C C NEVAL Integer. C The number of function evaluations. C IFAIL Integer. C IFAIL = 0 for normal exit. C IFAIL = 2 if NUMFUN is less than 1. C IFAIL = 3 if B(1) < A(1). C IFAIL = 4 if unlegal PERIOD. C IFAIL = 5 if MAXPTS is less than MINPTS or MINPTS < 21. C IFAIL = 6 if EPSABS <= 0 and EPSREL <= 0. C IFAIL = 7 if WRKSUB is too small. C IFAIL = 8 if unlegal value of EMAX. C IFAIL = 9 if unlegal RESTAR. C IFAIL = 10 if unlegal value of KEY. C IFAIL = 11 if NW is less than LIMIT (See the code DCHINF). C IFAIL = 12 if either PERIOD is wrong or B(1) is too small. C IFAIL = 13 if unable to estimate GAMMA (In case KEY=2 only) C C***ROUTINES CALLED DGAINF C***END PROLOGUE DCHINF C C Global variables. C INTEGER NUMFUN,MINPTS,MAXPTS,NW,RESTAR,EMAX,NEVAL,WRKSUB,KEY, +IFAIL DOUBLE PRECISION A(WRKSUB),B(WRKSUB),EPSABS,EPSREL,PERIOD,GAMMA EXTERNAL DGAINF,FUNSUB C C Local variables. C INTEGER LIMIT,NUMNUL,NUM DOUBLE PRECISION WIDTH PARAMETER (NUMNUL=16,NUM=21) C C***FIRST EXECUTABLE STATEMENT DCHINF C IFAIL=0 C C Check on legal NUMFUN. C IF (NUMFUN.LT.1) THEN IFAIL=2 RETURN END IF C C Check on legal length of the first interval of integration. C WIDTH=B(1)-A(1) IF (WIDTH.LT.0) THEN IFAIL=3 RETURN END IF C C Check on legal sign of PERIOD. C IF (PERIOD.LE.0) THEN IFAIL=4 RETURN END IF C C Check on MAXPTS >= MINPTS. C IF ((MAXPTS.LT.MINPTS).OR.(MINPTS.LT.NUM)) THEN IFAIL=5 RETURN END IF C C Check on legal accuracy requests. C IF (EPSABS.LE.0.AND.EPSREL.LE.0) THEN IFAIL=6 RETURN END IF C C Check on legal WRKSUB. C IF (NUM*WRKSUB.LT.MAXPTS) THEN IFAIL=7 RETURN END IF C C Check on legal value of EMAX. C IF (EMAX.LT.1) THEN IFAIL=8 RETURN END IF C C Check on legal RESTAR. C IF (RESTAR.NE.0.AND.RESTAR.NE.1) THEN IFAIL=9 RETURN ELSE IF (RESTAR.EQ.0) THEN NEVAL=0 END IF C C Check on legal KEY. C IF (KEY.NE.0.AND.KEY.NE.1.AND.KEY.NE.2) THEN IFAIL=10 RETURN END IF C C Check on big enough NW. C LIMIT=1+5*NUMFUN+4*WRKSUB*NUMFUN+3*WRKSUB+(EMAX+1)**2+(NUMNUL+1+ +NUM)*NUMFUN IF (NW.LT.LIMIT) THEN IFAIL=11 RETURN END IF C C If the user gives GAMMA <= 0 and KEY = 2 then DQAINF will try to C estimate GAMMA. C If the subroutine estimates GAMMA successfully within its error C bound, IFAIL will remain 0 while GAMMA is given the new value. C The routine DGAINF may detect errors: C C IFAIL = 12 indicates that either C B is too small or the PERIOD is wrong, or maybe the function does no C have the desired properties. (NOTE: The code checks ONLY the first C function in the vector, assuming that if this is correct then the C other functions share its properties.) C C IFAIL = 13 : No errors where detected, but the code was unable to C compute GAMMA sufficiently accurate. Remedy: use KEY = 1 instead. C IF ((GAMMA.LE.0).AND.(KEY.EQ.2)) THEN CALL DGAINF (NUMFUN,FUNSUB,B(1),PERIOD,GAMMA,NEVAL,IFAIL) RETURN END IF C RETURN C C***END DCHINF C END CUT HERE............ cat > dtrint.f <<'CUT HERE............' SUBROUTINE DTRINT (DVFLAG,SBRGNS,GREATE,LIST,NEW) C***BEGIN PROLOGUE DTRINT C***REFER TO DQAINF C***PURPOSE DTRINT maintains a heap of subregions. C***DESCRIPTION DTRINT maintains a heap of subregions. The subregions C are stored in a partially sorted binary tree, ordered to the C size of the greatest error estimates of each subregion(GREATE). C The subregion with greatest error estimate is in the first C position of the heap. C C PARAMETERS C C DVFLAG Integer. C If DVFLAG = 1, we remove the subregion with C greatest error from the heap. C If DVFLAG = 2, we insert a new subregion in the heap. C If DVFLAG = 3, we move the top element to a new position. C SBRGNS Integer. C Number of subregions in the heap. C GREATE Real array of dimension SBRGNS. C Used to store the greatest estimated errors in C all subregions. C LIST Integer array of dimension SBRGNS. C Used as a partially ordered list of pointers to the C different subregions. This list is a heap where the C element on top of the list is the subregion with the C greatest error estimate. C NEW Integer. C Index to the new region to be inserted in the heap. C***ROUTINES CALLED-NONE C***END PROLOGUE DTRINT C C Global variables. C INTEGER DVFLAG,NEW,SBRGNS,LIST(*) DOUBLE PRECISION GREATE(*) C C Local variables. C C GREAT is used as intermediate storage for the greatest error of a C subregion. C PNTR Pointer. C SUBRGN Position of child/parent subregion in the heap. C SUBTMP Position of parent/child subregion in the heap. INTEGER SUBRGN,SUBTMP,PNTR DOUBLE PRECISION GREAT C C***FIRST EXECUTABLE STATEMENT DTRINT C C If DVFLAG = 1, we will reduce the partial ordered list by the C element with greatest estimated error. Thus the element in C in the heap with index LIST(1) is vacant and can be used later. C Reducing the heap by one element implies that the last element C should be re-positioned. C If DVFLAG = 3, we will keep the size of the partially ordered C list but re-position the first element. C IF (DVFLAG.NE.2) THEN IF (DVFLAG.EQ.1) THEN PNTR=LIST(SBRGNS) GREAT=GREATE(PNTR) SBRGNS=SBRGNS-1 ELSE IF (DVFLAG.EQ.3) THEN PNTR=LIST(1) GREAT=GREATE(PNTR) END IF SUBRGN=1 10 SUBTMP=2*SUBRGN IF (SUBTMP.LE.SBRGNS) THEN IF (SUBTMP.NE.SBRGNS) THEN C C Find max. of left and right child. C IF (GREATE(LIST(SUBTMP)).LT.GREATE(LIST(SUBTMP+1))) THEN SUBTMP=SUBTMP+1 END IF END IF C C Compare max.child with parent. C If parent is max., then done. C IF (GREAT.LT.GREATE(LIST(SUBTMP))) THEN C C Move the pointer at position subtmp up the heap. C LIST(SUBRGN)=LIST(SUBTMP) SUBRGN=SUBTMP GO TO 10 END IF END IF C C Update the pointer. C IF (SBRGNS.GT.0) THEN LIST(SUBRGN)=PNTR END IF ELSE IF (DVFLAG.EQ.2) THEN C C If DVFLAG = 2, find the position for the NEW region in the heap. C GREAT=GREATE(NEW) SUBRGN=SBRGNS 20 SUBTMP=SUBRGN/2 IF (SUBTMP.GE.1) THEN C C Compare max.child with parent. C If parent is max, then done. C IF (GREAT.GT.GREATE(LIST(SUBTMP))) THEN C C Move the pointer at position subtmp down the heap. C LIST(SUBRGN)=LIST(SUBTMP) SUBRGN=SUBTMP GO TO 20 END IF END IF C C Set the pointer to the new region in the heap. C LIST(SUBRGN)=NEW END IF C C***END DTRINT C RETURN END CUT HERE............ cat > dgaint.f <<'CUT HERE............' SUBROUTINE DGAINT (A,B,NUMFUN,FUNSUB,DONE,MORE,EPMACH,NULL, +FUNVAL,BASABS,FLAG,BASVAL,RGNERR,GREATE,WORST,C,D) C***BEGIN PROLOGUE DGAINT C***PURPOSE To compute basic integration rule values and C corresponding error estimates. C C***LAST MODIFICATION 94-03-10 C C***REFER TO DQAINF C C***DESCRIPTION DGAINT computes basic integration rule values C for a vector of integrands over an interval. C DGAINT also computes estimates for the errors by C using several null rule approximations. C ON ENTRY C C A Real. C B Real. C The left and right endpoints of the interval. C NUMFUN Integer. C Number of components of the vector integrand. C FUNSUB Externally declared subroutine for computing C all components of the integrand at the given C evaluation point. C It must have parameters (X,NUMFUN,NP,FUNVLS) C Input parameters: C X The x-coordinates of all 21 evaluation points. C NUMFUN Integer that defines the number of C components of I. C NP Integer that gives the number of evaluation points C Output parameter: C FUNVLS Real array of dimension (21,NUMFUN) C that defines NUMFUN components of the integrand C evaluated in all 21 evaluation points. C C DONE Logical array of dimension NUMFUN. C DONE(I)=.TRUE. if function number I has been C integrated to the specified precision, else DONE(I)=.FALSE. C MORE Logical. C Information about the fact that there is still work to do. C EPMACH Real. C Machine epsilon. C NULL Real array of dimension (16,NUMFUN). C A work array. C FUNVAL Real array of dimension (21,NUMFUN). C A work array. C BASABS Real array of dimension NUMFUN. C A work array. C ON RETURN C C FLAG Integer. C Signals that at least one interval has become too small C to distinguish between a rule evaluation point and an C endpoint. C C BASVAL Real array of dimension NUMFUN. C The values for the basic rule for each component C of the integrand. C RGNERR Real array of dimension NUMFUN. C The error estimates for each component of the integrand. C GREATE Real. C The greatest error component of the integrand. C WORST Index to the greatest error in each sub-interval. C C C Real. C D Real C C and D are subdivision points based on C uniform trisection: C = A + (B-A)/3 and D = B - (B-a)/3. C C***REFERENCES Espelid, T. O., Integration Rules, Null Rules and C Error Estimation, Report no 33, Informatics, Univ. of Bergen, C 1988. C C Berntsen,J. and Espelid,T.O., Error estimation in Automatic C Quadrature Routines, ACM Trans. Math. Softw., C 17, 2, 233-252, 1991. C C Espelid, T. O., DQAINT: An Algorithm for Adaptive Quadrature C over a Collection of Finite Intervals, in Numerical C Integration, Recent Developments, Software and Applications, C Espelid and Genz (eds), NATO ASI Series C: C Mathematical and Physical Sciences - Vol. 357, C Kluwer, Dordrecht, The Netherlands, pages 341-342, 1992. C C***ROUTINES CALLED FUNSUB C***END PROLOGUE DGAINT C C Global variables. C EXTERNAL FUNSUB LOGICAL DONE(NUMFUN),MORE INTEGER NUMFUN,FLAG,WORST DOUBLE PRECISION A,B,BASVAL(NUMFUN),RGNERR(NUMFUN),NULL(16, +NUMFUN),GREATE,C,D,BASABS(NUMFUN),FUNVAL(21,NUMFUN),EPMACH C C Local variables. C C WTLENG The number of weights of the integration rule. C G Real array of dimension WTLENG. C The coordinates for the generators of C the evaluation points. C The integration rule is using symmetric evaluation C points and has the structure (1,6,3). That is, C 1 point of multiplicity 1, C 6 sets of points of multiplicity 3 and C 3 sets of points of multiplicity 6. C This gives totally 37 evaluation points. C In order to reduce the number of loops in DGAINT, C the 3 loops for the sets of multiplicity 6 are split C into 6 loops and added to the loops for the sets of C multiplicity 3. C The number of weights we have to give with C this splitting is 13(WTLENG). C C W Real array of dimension (9,WTLENG). C The weights of the basic rule and the null rules. C W(1),...,W(WTLENG) are weights for the basic rule. C C X Real array of dimension 21. C Evaluation points in (A,B). C HLEN Real. C Half of intervals length INTEGER I,I1,I2,J,K,WTLENG,DGE1,DEGREE,NUMNUL,DIMNUL PARAMETER (NUMNUL=8,DIMNUL=16) DOUBLE PRECISION X(21),RR(NUMNUL/2-1),R,NOISE,E(NUMNUL/2),EMAX, +ALPHA,CONST,W(11),ABSCIS,HLEN,SAFETY,RCRIT,ABSSUM,SUM,DIFF,CENTR, +G(11),FACTOR,NULLW(16,11) PARAMETER (WTLENG=11,FACTOR=0.02D+0,SAFETY=10.D0,RCRIT=0.5D0,DGE1= +18,DEGREE=41,ALPHA=(DEGREE-DGE1)/2.D0) C C The abscissas to the 21 point Gauss integration rule. C The abscissas are given in -1,1: due to the symmetry we only specify C values in 0,1. C DATA (G(I),I=1,11)/0.9937521706203895D0,0.9672268385663062D0,0. +9200993341504008D0,0.8533633645833172D0,0.7684399634756779D0,0. +6671388041974123D0,0.5516188358872198D0,0.4243421202074387D0,0. +2880213168024010D0,0.1455618541608950D0,0.000000000000000D0/ C C Weights of the 21 point Gauss quadrature rule. Same remark C about symmetry. C DATA (W(I),I=1,11)/0.01601722825777433D0,0.03695378977085249D0,0. +05713442542685720D0,0.07610011362837930D0,0.09344442345603386D0,0. +1087972991671483D0,0.1218314160537285D0,0.1322689386333374D0,0. +1398873947910731D0,0.1445244039899700D0,0.1460811336496904D0/ C C Weights of the 5 symmetric nullrules of degrees 19,17,15,13,11 C Weights of the 5 asymmetric nullrules of degrees 18,16,14,12,10 C Same remark about symmetry as above. The nullrules are all C orthogonal and has the same norm as the basic rule and are C given with decreasingly degrees. C DATA NULLW(1,1)/0.5859224694730026D-02/ DATA NULLW(1,2)/-0.2024707000741622D-01/ DATA NULLW(1,3)/0.3883580247121445D-01/ DATA NULLW(1,4)/-0.5965412595242497D-01/ DATA NULLW(1,5)/0.8114279498343020D-01/ DATA NULLW(1,6)/-0.1019231472030305D+00/ DATA NULLW(1,7)/0.1207652963052454D+00/ DATA NULLW(1,8)/-0.1366043796711610D+00/ DATA NULLW(1,9)/0.1485698262567817D+00/ DATA NULLW(1,10)/-0.1560150459859118D+00/ DATA NULLW(1,11)/0.1585416482170856D+00/ DATA NULLW(2,1)/0.1301976799706014D-01/ DATA NULLW(2,2)/-0.4379005851020851D-01/ DATA NULLW(2,3)/0.7990096087086482D-01/ DATA NULLW(2,4)/-0.1138307201442027D+00/ DATA NULLW(2,5)/0.1394263659262734D+00/ DATA NULLW(2,6)/-0.1520456605731098D+00/ DATA NULLW(2,7)/0.1489588260146731D+00/ DATA NULLW(2,8)/-0.1296181347851887D+00/ DATA NULLW(2,9)/0.9568420420614478D-01/ DATA NULLW(2,10)/-0.5078074459100106D-01/ DATA NULLW(2,11)/0.0000000000000000D+00/ DATA NULLW(3,1)/0.2158184987561927D-01/ DATA NULLW(3,2)/-0.6965227115767195D-01/ DATA NULLW(3,3)/0.1174438965053943D+00/ DATA NULLW(3,4)/-0.1473794001233916D+00/ DATA NULLW(3,5)/0.1481989091733945D+00/ DATA NULLW(3,6)/-0.1168273220215079D+00/ DATA NULLW(3,7)/0.5890214560028095D-01/ DATA NULLW(3,8)/0.1273585156460484D-01/ DATA NULLW(3,9)/-0.8133037350927629D-01/ DATA NULLW(3,10)/0.1304777802379009D+00/ DATA NULLW(3,11)/-0.1483021322906934D+00/ DATA NULLW(4,1)/0.3119657617222001D-01/ DATA NULLW(4,2)/-0.9516116459787523D-01/ DATA NULLW(4,3)/0.1431705707841666D+00/ DATA NULLW(4,4)/-0.1462171986707822D+00/ DATA NULLW(4,5)/0.9677919508585969D-01/ DATA NULLW(4,6)/-0.1075583592960879D-01/ DATA NULLW(4,7)/-0.7936141880173113D-01/ DATA NULLW(4,8)/0.1380749552472930D+00/ DATA NULLW(4,9)/-0.1417577117227090D+00/ DATA NULLW(4,10)/0.8867798725104829D-01/ DATA NULLW(4,11)/0.0000000000000000D+00/ DATA NULLW(5,1)/0.4157639307601386D-01/ DATA NULLW(5,2)/-0.1179114894020921D+00/ DATA NULLW(5,3)/0.1511566572815612D+00/ DATA NULLW(5,4)/-0.1073644825716617D+00/ DATA NULLW(5,5)/0.4174411212397235D-02/ DATA NULLW(5,6)/0.1012057375471417D+00/ DATA NULLW(5,7)/-0.1472858866418607D+00/ DATA NULLW(5,8)/0.1063772962797608D+00/ DATA NULLW(5,9)/-0.2323645744823691D-02/ DATA NULLW(5,10)/-0.1030910117645103D+00/ DATA NULLW(5,11)/0.1469720414561505D+00/ DATA NULLW(6,1)/0.5246625962340516D-01/ DATA NULLW(6,2)/-0.1358182960749584D+00/ DATA NULLW(6,3)/0.1386355746642898D+00/ DATA NULLW(6,4)/-0.3967547044517777D-01/ DATA NULLW(6,5)/-0.8983329656061084D-01/ DATA NULLW(6,6)/0.1471801958801420D+00/ DATA NULLW(6,7)/-0.8524048604745531D-01/ DATA NULLW(6,8)/-0.4617298114857220D-01/ DATA NULLW(6,9)/0.1397282757969823D+00/ DATA NULLW(6,10)/-0.1185867937861861D+00/ DATA NULLW(6,11)/0.0000000000000000D+00/ DATA NULLW(7,1)/0.6363031447247006D-01/ DATA NULLW(7,2)/-0.1472028208379138D+00/ DATA NULLW(7,3)/0.1063757528761779D+00/ DATA NULLW(7,4)/0.3881687506910375D-01/ DATA NULLW(7,5)/-0.1432999618142209D+00/ DATA NULLW(7,6)/0.9698969424297650D-01/ DATA NULLW(7,7)/0.5209450556829039D-01/ DATA NULLW(7,8)/-0.1455658951943161D+00/ DATA NULLW(7,9)/0.8343286549711612D-01/ DATA NULLW(7,10)/0.6800562635441401D-01/ DATA NULLW(7,11)/-0.1465539124681842D+00/ DATA NULLW(8,1)/0.7484599067063009D-01/ DATA NULLW(8,2)/-0.1508776892345244D+00/ DATA NULLW(8,3)/0.5853283458897407D-01/ DATA NULLW(8,4)/0.1062457136342151D+00/ DATA NULLW(8,5)/-0.1318732622123368D+00/ DATA NULLW(8,6)/-0.1673118490576078D-01/ DATA NULLW(8,7)/0.1428979325476485D+00/ DATA NULLW(8,8)/-0.7818432195969258D-01/ DATA NULLW(8,9)/-0.9112601052788798D-01/ DATA NULLW(8,10)/0.1382849496064090D+00/ DATA NULLW(8,11)/0.0000000000000000D+00/ DATA NULLW(9,1)/0.8590167508061712D-01/ DATA NULLW(9,2)/-0.1462121283895959D+00/ DATA NULLW(9,3)/0.1973066649848703D-02/ DATA NULLW(9,4)/0.1434120884358229D+00/ DATA NULLW(9,5)/-0.6050045998747565D-01/ DATA NULLW(9,6)/-0.1192968264851738D+00/ DATA NULLW(9,7)/0.1063582979588903D+00/ DATA NULLW(9,8)/0.7871971420591506D-01/ DATA NULLW(9,9)/-0.1360664606736734D+00/ DATA NULLW(9,10)/-0.2747426951466367D-01/ DATA NULLW(9,11)/0.1463706054390223D+00/ DATA NULLW(10,1)/0.9659618304508728D-01/ DATA NULLW(10,2)/-0.1331667766592828D+00/ DATA NULLW(10,3)/-0.5483608118503819D-01/ DATA NULLW(10,4)/0.1395396581193282D+00/ DATA NULLW(10,5)/0.3842219337177220D-01/ DATA NULLW(10,6)/-0.1430606163131484D+00/ DATA NULLW(10,7)/-0.2498840938693774D-01/ DATA NULLW(10,8)/0.1451753725543706D+00/ DATA NULLW(10,9)/0.1236834040097046D-01/ DATA NULLW(10,10)/-0.1461902970879641D+00/ DATA NULLW(10,11)/0.0000000000000000D+00/ DATA NULLW(11,1)/0.1067392105869384D+00/ DATA NULLW(11,2)/-0.1122928178247447D+00/ DATA NULLW(11,3)/-0.1031959020477783D+00/ DATA NULLW(11,4)/0.9558143541939021D-01/ DATA NULLW(11,5)/0.1196951864405313D+00/ DATA NULLW(11,6)/-0.7225984000378730D-01/ DATA NULLW(11,7)/-0.1339424732473705D+00/ DATA NULLW(11,8)/0.4492456833975673D-01/ DATA NULLW(11,9)/0.1431238351576778D+00/ DATA NULLW(11,10)/-0.1523606417516131D-01/ DATA NULLW(11,11)/-0.1462742772906911D+00/ DATA NULLW(12,1)/0.1161523191050730D+00/ DATA NULLW(12,2)/-0.8469391287377601D-01/ DATA NULLW(12,3)/-0.1355896186086413D+00/ DATA NULLW(12,4)/0.2408868272651161D-01/ DATA NULLW(12,5)/0.1460359069105494D+00/ DATA NULLW(12,6)/0.4632194803727984D-01/ DATA NULLW(12,7)/-0.1231814607655799D+00/ DATA NULLW(12,8)/-0.1068762140630544D+00/ DATA NULLW(12,9)/0.7029919038187181D-01/ DATA NULLW(12,10)/0.1416700606593806D+00/ DATA NULLW(12,11)/0.0000000000000000D+00/ DATA NULLW(13,1)/0.1246701955330830D+00/ DATA NULLW(13,2)/-0.5195253287985397D-01/ DATA NULLW(13,3)/-0.1469123277046623D+00/ DATA NULLW(13,4)/-0.5433985749387003D-01/ DATA NULLW(13,5)/0.1052913655334742D+00/ DATA NULLW(13,6)/0.1341759283651172D+00/ DATA NULLW(13,7)/-0.2310968036052471D-02/ DATA NULLW(13,8)/-0.1358135230198954D+00/ DATA NULLW(13,9)/-0.1024826424015055D+00/ DATA NULLW(13,10)/0.5656562071840163D-01/ DATA NULLW(13,11)/0.1462174827723311D+00/ DATA NULLW(14,1)/0.1321420280297885D+00/ DATA NULLW(14,2)/-0.1602500237425218D-01/ DATA NULLW(14,3)/-0.1353193782985104D+00/ DATA NULLW(14,4)/-0.1170027382391999D+00/ DATA NULLW(14,5)/0.1614011431557236D-01/ DATA NULLW(14,6)/0.1330641979525424D+00/ DATA NULLW(14,7)/0.1205891076794731D+00/ DATA NULLW(14,8)/-0.8640919108146020D-02/ DATA NULLW(14,9)/-0.1294253464460428D+00/ DATA NULLW(14,10)/-0.1251272318395094D+00/ DATA NULLW(14,11)/0.0000000000000000D+00/ DATA NULLW(15,1)/0.1384328909795413D+00/ DATA NULLW(15,2)/0.2088816813445404D-01/ DATA NULLW(15,3)/-0.1025551817987029D+00/ DATA NULLW(15,4)/-0.1456993480020755D+00/ DATA NULLW(15,5)/-0.8041833842953963D-01/ DATA NULLW(15,6)/0.4369891359834745D-01/ DATA NULLW(15,7)/0.1355713757017371D+00/ DATA NULLW(15,8)/0.1284341564046552D+00/ DATA NULLW(15,9)/0.2777799655524739D-01/ DATA NULLW(15,10)/-0.9304002613430708D-01/ DATA NULLW(15,11)/-0.1461812140165950D+00/ DATA NULLW(16,1)/0.1434250647895144D+00/ DATA NULLW(16,2)/0.5648832390790171D-01/ DATA NULLW(16,3)/-0.5370731005946019D-01/ DATA NULLW(16,4)/-0.1320553898020212D+00/ DATA NULLW(16,5)/-0.1399117916675364D+00/ DATA NULLW(16,6)/-0.7464504070837483D-01/ DATA NULLW(16,7)/0.2922259459390760D-01/ DATA NULLW(16,8)/0.1177993871727999D+00/ DATA NULLW(16,9)/0.1454239155303997D+00/ DATA NULLW(16,10)/0.9797588771824906D-01/ DATA NULLW(16,11)/0.0000000000000000D+00/ C C***FIRST EXECUTABLE STATEMENT DGAINT C C Define constants C CONST=SAFETY*RCRIT**(1-ALPHA) C C C Compute half-length and center of interval. C HLEN=(B-A)/2 CENTR=A+HLEN X(11)=CENTR C C Compute all abscissas C DO 10 I=1,10 ABSCIS=HLEN*G(I) X(I)=CENTR-ABSCIS X(22-I)=CENTR+ABSCIS 10 CONTINUE IF ((X(21).EQ.X(20)).OR.(X(1).EQ.X(2))) THEN FLAG=1 END IF C C Evaluate all NUMFUN functions in the 21 evaluation points. C CALL FUNSUB (X,NUMFUN,21,FUNVAL) C C Apply the basic rule: first center point C DO 20 J=1,NUMFUN BASVAL(J)=W(11)*FUNVAL(11,J) BASABS(J)=ABS(BASVAL(J)) 20 CONTINUE C C Apply all nullrules: first center point C DO 40 J=1,NUMFUN DO 30 I=1,NUMNUL NULL(I,J)=NULLW(I,11)*FUNVAL(11,J) 30 CONTINUE 40 CONTINUE C C Compute the basic rule contributions from the other points. C DO 60 J=1,NUMFUN DO 50 I=1,10 SUM=FUNVAL(I,J)+FUNVAL(22-I,J) ABSSUM=ABS(FUNVAL(I,J))+ABS(FUNVAL(22-I,J)) BASVAL(J)=W(I)*SUM+BASVAL(J) BASABS(J)=BASABS(J)+W(I)*ABSSUM 50 CONTINUE 60 CONTINUE C C Compute the null rule contributions from the other points. C Recall that one half of the nullrules is symmetric and the other C half is asymmetric. C DO 90 J=1,NUMFUN DO 80 K=1,10 SUM=FUNVAL(K,J)+FUNVAL(22-K,J) DIFF=FUNVAL(K,J)-FUNVAL(22-K,J) DO 70 I=1,NUMNUL/2 I2=2*I I1=I2-1 NULL(I1,J)=NULL(I1,J)+NULLW(I1,K)*SUM NULL(I2,J)=NULL(I2,J)+NULLW(I2,K)*DIFF 70 CONTINUE 80 CONTINUE 90 CONTINUE C C Scale the results with the length of the interval C DO 100 J=1,NUMFUN BASVAL(J)=HLEN*BASVAL(J) BASABS(J)=HLEN*BASABS(J) 100 CONTINUE C C Compute errors. C GREATE=0 DO 130 J=1,NUMFUN EMAX=0 DO 110 I=1,NUMNUL/2 I2=2*I I1=I2-1 E(I)=HLEN*DSQRT(NULL(I1,J)**2+NULL(I2,J)**2) EMAX=MAX(EMAX,E(I)) 110 CONTINUE R=0 DO 120 I=1,NUMNUL/2-1 IF (E(I+1).NE.0) THEN RR(I)=E(I)/E(I+1) ELSE RR(I)=1 END IF R=MAX(R,RR(I)) 120 CONTINUE IF (R.GE.1) THEN RGNERR(J)=SAFETY*EMAX ELSE IF (R.GE.RCRIT) THEN RGNERR(J)=SAFETY*R*E(1) ELSE RGNERR(J)=CONST*(R**ALPHA)*E(1) END IF C C Check the noise level. C NOISE=50*EPMACH*BASABS(J) IF ((E(1).LT.NOISE).AND.(E(2).LT.NOISE)) THEN RGNERR(J)=NOISE ELSE RGNERR(J)=MAX(NOISE,RGNERR(J)) END IF C C Check if this is the greatest error so far among the remaining C problems. One exception: If the user wants to force the code C to do extra work (which seems unnecessary). C IF (.NOT.(MORE.AND.DONE(J)).AND.(RGNERR(J).GT.GREATE)) THEN GREATE=RGNERR(J) WORST=J END IF 130 CONTINUE C C Compute the new subdivision points. C C=A+(B-A)/3 D=B-(B-A)/3 RETURN C C***END DGAINT C END CUT HERE............ cat > dgainf.f <<'CUT HERE............' SUBROUTINE DGAINF (NUMFUN,F,B,PERIOD,GAMMA,NEVAL,IFAIL) C C***BEGIN PROLOGUE DGAINF C***KEYWORDS Linear extrapolation, C oscillating functions, asymptotic behavior. C***PURPOSE To compute an estimate of the exponent GAMMA which governs C the asymptotic behavior of an oscillating function: C F(X+PERIOD*K) appr. C(X)/(X+PERIOD*K)**GAMMA, for all X C and all integers K sufficiently great. C The function is assumed to C have asymptotical oscillating behavior and C F(X+PERIOD*(K+1/2)) approx. - F(X+PERIOD*K). C C***LAST MODIFICATION 93-01-31 C C***REFER TO DQAINF C C***DESCRIPTION Computation of the exponent GAMMA > 0 which governs C the asymptotic behavior of the functions F. All functions in the C vector F are assumed to be asymptotically periodic with the same C PERIOD and the same asymptotical behavior. Furthermore we assume C that for all functions in the vector we have C F(X+PERIOD*K) approx. C(X)/(X+PERIOD*K)**GAMMA, for all values of X C and K, such that (X+PERIOD*K) >> B. C may be zero for some values of C X. Assuming that C is non-zero then we have that C F(X+2*PERIOD*K)/F(X+PERIOD*K) approx. (1 + D)**GAMMA, with C D = PERIOD*K/(X+PERIOD*K). C C ON ENTRY C C NUMFUN Integer. C The number of components of the integral. C F Externally declared subroutine for computing C all components of the integrand at all NP C evaluation points. C It must have parameters (X,NUMFUN,NP,FUNVLS) C Input parameters: C X The NP x-coordinates of the evaluation points. C NUMFUN Integer that defines the number of C components of I. C NP Integer that gives the number of evaluation points C Output parameter: C FUNVLS Real array of dimension (NUMFUN,NP) C that defines NUMFUN components of the integrand fo C each evaluation point. C C B Real. C Asymptotic behavior assumed satisfied for all X >= B and C for all the integrands. C PERIOD Real. C All functions are assumed to be asymptotically periodic, C and with the same period. C C ON RETURN C C GAMMA Real C Estimated value of the rate of decay of the functions C when X tends to +infinity. C NEVAL Integer C The number of function evaluations. C IFAIL Integer. C IFAIL = 0 for normal exit. C IFAIL = 12 if either PERIOD is wrong or B(1) is too small. C IFAIL = 13 if unable to estimate GAMMA (In case KEY=2 only) C C***ROUTINES CALLED F,D1MACH C***END PROLOGUE DGAINF C C Global variables. C EXTERNAL F,D1MACH DOUBLE PRECISION PERIOD,GAMMA,B,D1MACH INTEGER NUMFUN,IFAIL,NEVAL C C Local variables. C DOUBLE PRECISION P,X,H,XB,FX,FXP,XR,FR,EPMACH,FOLD,FNEW,SIGN,XOLD, +XNEW INTEGER I,IMX,M PARAMETER (IMX=15) DOUBLE PRECISION A(0:IMX),DA,FACTOR(0:IMX),FI(0:IMX),DF,LAST C C***FIRST EXECUTABLE STATEMENT DGAINF C EPMACH=D1MACH(4) C C First: compute a starting point X >= B where the function value F(X) C is reasonably large in absolute value. Remark: We simply want to C avoid that F(X) becomes approximately 0. C XOLD=B CALL F (XOLD,1,1,FOLD) P=PERIOD/4 XNEW=XOLD DO 10 I=1,3,1 XNEW=XNEW+P CALL F (XNEW,1,1,FNEW) IF (ABS(FOLD).LT.ABS(FNEW)) THEN XOLD=XNEW IF ((FOLD*FNEW).GT.0) THEN P=-P END IF FOLD=FNEW ELSE IF ((FOLD*FNEW).GT.0) THEN XNEW=XOLD END IF END IF P=P/2 10 CONTINUE X=XOLD FX=FOLD C C Check that the function values do change sign C and that the absolute value of the function is decreasing C when evaluated at points which differ by PERIOD*integer/2. C CALL F (X+PERIOD/2,1,1,FXP) NEVAL=5 LAST=FXP SIGN=FX FI(0)=X IF (((SIGN*FXP).LT.0).AND.(ABS(FX).GT.ABS(FXP))) THEN FR=-FXP/FX XR=X/(X+PERIOD/2) A(0)=LOG(FR)/LOG(XR) FACTOR(0)=A(0)*(1/ABS(FR*LOG(FR))+1/ABS(XR*LOG(XR)))*EPMACH ELSE IFAIL=12 RETURN END IF C C Compute a sequence of estimates to GAMMA and extrapolate. C C FACTOR is meant to keep track of the increase in the error due to C the extrapolation. We assume that PERIOD is known correctly to full C double precision. If X is an exact local extremum of F(X) then it is C not important to know PERIOD with such a high precision. C However, since X is just a point where ABS(F(X)) is not too small C the quality of PERIOD becomes essential when we estimate GAMMA. C H=PERIOD/2 XB=X DO 30 I=1,IMX H=2*H X=XB+H FI(I)=X CALL F (X+PERIOD/2,1,1,FXP) CALL F (X,1,1,FX) NEVAL=NEVAL+2 C C Check the computed function values. C IF (((SIGN*FX).GT.0).AND.((SIGN*FXP).LT.0).AND.(ABS(FX).GT. + ABS(FXP)).AND.(ABS(LAST).GT.ABS(FX))) THEN LAST=FXP FR=-FXP/FX XR=X/(X+PERIOD/2) A(I)=LOG(FR)/LOG(XR) FACTOR(I)=A(0)*(2**I/ABS(FR*LOG(FR))+1/ABS(XR*LOG(XR)))* + EPMACH ELSE C C Either PERIOD is wrong or B is too small. C IFAIL=12 RETURN END IF DO 20 M=I-1,0,-1 DA=(A(M+1)-A(M))/((FI(I)-FI(M))/FI(M)) DF=(FACTOR(M+1)+FACTOR(M))/((FI(I)-FI(M))/FI(M)) A(M)=A(M+1)+DA FACTOR(M)=FACTOR(M+1)+DF 20 CONTINUE IF (ABS(DA).LE.FACTOR(0)) THEN GAMMA=A(0) C C Normal return C RETURN END IF 30 CONTINUE C C We did not succeed in estimating GAMMA to sufficient precision. C GAMMA=A(0) IFAIL=13 RETURN C C***END DGAINF C END CUT HERE............ cat > dexinf.f <<'CUT HERE............' SUBROUTINE DEXINF (NUMFUN,GAMMA,C,KEY,N,T,CT,UPDATE,EMAX,U,E, +RESULT,ABSERR,EXTERR,BETA,MAXEER) C***BEGIN PROLOGUE DEXINF C***KEYWORDS Quasi-linear extrapolation, infinite integration, C oscillating functions, error estimation. C***PURPOSE To compute better estimates to a vector of approximations C to integrals of oscillating functions over an infinite C interval and to provide new and updated error estimates. C C***LAST MODIFICATION 94-03-11 C C***DESCRIPTION The routine uses linear extrapolation to compute better C approximations to each component in a vector of C integrals. All components are assumed to be integrals of C oscillating functions which asymptotically have the same C behavior at infinity. C A series approach is used, assuming C that the terms are given with estimates of the error in C each term. New error estimates are computed too. The C routine have two options: EITHER a new extrapolation term C is provided and we take a new extrapolation step, C OR we update one previously computed term in the series and C therefore have to update the extrapolation tableau. C C ON ENTRY C C NUMFUN Integer. C The number of components of the integral. C GAMMA Real. C Asymptotic decay of the functions; 1/x**(GAMMA), GAMMA>0. C C Real. C Reference point: User specified subdivision point/period. C KEY Integer C Choice of extrapolation method: C KEY = 0 : Euler's method. C KEY = 1 : Modified Euler. C KEY = 2 : Overholt's polynomial order 2-method. C This last method is the only one that needs C the value of GAMMA. C N Integer C The number of U-terms in the series: 0,1,2,...,N. C Makes it possible to perform N extrapolation steps. C T Real array of dimension (NUMFUN,0:EMAX) C Contains the last row in the extrapolation tableau for C each function in the vector. C CT Real array of dimension NUMFUN. C Contains the element to the left of the diagonal element C in the extrapolation tableau. C UPDATE Integer C If < 0 then this is a new extrapolation step. C If >= 1 then this is a step where we have to correct the C existing tableau. The value of UPDATE gives the index to C the U-element that has been modified. C EMAX Integer C The maximum allowed number of extrapolation steps. C U Real array of dimension (NUMFUN,0:N) C C E Real array of dimension (NUMFUN,0:N) C The estimated errors of all U-terms in the series. C ON RETURN C C T Real array of dimension (NUMFUN:0:EMAX) C Contains the diagonal element in the extrapolation tableau C for each function in the vector after the extrapolation. C CT Real array of dimension NUMFUN. C Contains the element to the left of the diagonal element C in the extrapolation tableau. C RESULT Real array of dimension NUMFUN C Contains the new approximations for the NUMFUN components C of the integral. C ABSERR Real array of dimension NUMFUN. C Returns the global errors for all components. C This includes both the pure extrapolation C error and the effect of not knowing the U-terms exactly. C EXTERR Real array of dimension NUMFUN. C These errors are the pure extrapolation errors. C BETA Real Array of dimension (EMAX +1)*(EMAX +1) C A table of coefficients to be used in the extrapolation C and the error estimation. C MAXEER Real. C The maximum extrapolation error. C***REFERENCES C C T.O.Espelid and K.J.Overholt, DQAINF: An Algorithm for Automatic C Integration of Infinite Oscillating Tails, C submitted for publication 1993. C C K.J.Overholt: The P-algorithms for extrapolation, C Department of Informatics, Report No. 36, 1989. C C***ROUTINES CALLED NONE C***END PROLOGUE DEXINF C C Global variables. C INTEGER N,EMAX,UPDATE,NUMFUN,KEY DOUBLE PRECISION GAMMA,T(NUMFUN),CT(NUMFUN),C,U(NUMFUN,0:N), +E(NUMFUN,0:N),BETA(0:EMAX,0:EMAX),MAXEER,RESULT(NUMFUN), +ABSERR(NUMFUN),EXTERR(NUMFUN) C C Local variables. C INTEGER I,J DOUBLE PRECISION SAVE1,SAVE2,CONST,BASE1,BASE2,P,Q PARAMETER (CONST=1.D0) C C CONST heuristic constant used in the error estimation. C C***FIRST EXECUTABLE STATEMENT DEXINF C C Initialize the beta-factors. C BETA(0,0)=1 IF (KEY.EQ.2) THEN BASE1=GAMMA BASE2=MAX(GAMMA-1,4*C) P=2.D0 ELSE BASE1=0.D0 BASE2=MAX(0.D0,4*C) P=1.D0 END IF C C Compute the new extrapolation coefficients. C IF (UPDATE.LT.0) THEN BETA(0,N)=1 BETA(N,0)=1 DO 20 I=1,N-1 SAVE1=1 DO 10 J=N-I,N-1 IF (KEY.EQ.0) THEN Q=0.5D0 ELSE Q=(1+(-BASE1-P*J)/(2*N+BASE2))/2.D0 END IF SAVE2=SAVE1-(SAVE1-BETA(I,J))*Q BETA(I,J)=SAVE1 10 SAVE1=SAVE2 BETA(I,N)=SAVE1 20 CONTINUE DO 30 J=1,N IF (KEY.EQ.0) THEN Q=0.5D0 ELSE Q=(1+(-BASE1-P*(J-1))/(2*N+BASE2))/2.D0 END IF 30 BETA(N,J)=(1-Q)*BETA(N,J-1) END IF C C Extrapolate; use the extrapolation coefficients. C DO 40 J=1,NUMFUN T(J)=0 CT(J)=0 DO 40 I=N,0,-1 CT(J)=CT(J)+BETA(I,N-1)*U(J,I) 40 T(J)=T(J)+BETA(I,N)*U(J,I) C C Then compute the error estimates. C First the extrapolation error and then the U-effect. C The maximum extrapolation error is computed too. C MAXEER=0 DO 60 J=1,NUMFUN EXTERR(J)=CONST*ABS(T(J)-CT(J))/Q MAXEER=MAX(EXTERR(J),MAXEER) C C Note: The last U-errors are effected by the extrapolation-process C Accumulate the total error in exterr. C DO 50 I=0,N 50 EXTERR(J)=EXTERR(J)+BETA(I,N)*E(J,I) 60 CONTINUE C C Define the results: update the results only if the error is improved C DO 70 J=1,NUMFUN IF ((EXTERR(J).LE.ABSERR(J)).OR.(ABSERR(J).EQ.0)) THEN RESULT(J)=T(J) ABSERR(J)=EXTERR(J) END IF 70 CONTINUE RETURN C C***END DEXINF C END CUT HERE............ cat > d1mach.f <<'CUT HERE............' DOUBLE PRECISION FUNCTION D1MACH(I) C C DOUBLE-PRECISION MACHINE CONSTANTS C C D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. C C D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. C C D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING. C C D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING. C C D1MACH( 5) = LOG10(B) C C TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT, C THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY C REMOVING THE C FROM COLUMN 1. C ON RARE MACHINES A STATIC STATEMENT MAY NEED TO BE ADDED. C (BUT PROBABLY MORE SYSTEMS PROHIBIT IT THAN REQUIRE IT.) C C FOR IEEE-ARITHMETIC MACHINES (BINARY STANDARD), ONE OF THE FIRST C TWO SETS OF CONSTANTS BELOW SHOULD BE APPROPRIATE. C C WHERE POSSIBLE, DECIMAL, OCTAL OR HEXADECIMAL CONSTANTS ARE USED C TO SPECIFY THE CONSTANTS EXACTLY. SOMETIMES THIS REQUIRES USING C EQUIVALENT INTEGER ARRAYS. IF YOUR COMPILER USES HALF-WORD C INTEGERS BY DEFAULT (SOMETIMES CALLED INTEGER*2), YOU MAY NEED TO C CHANGE INTEGER TO INTEGER*4 OR OTHERWISE INSTRUCT YOUR COMPILER C TO USE FULL-WORD INTEGERS IN THE NEXT 5 DECLARATIONS. C INTEGER SMALL(4) INTEGER LARGE(4) INTEGER RIGHT(4) INTEGER DIVER(4) INTEGER LOG10(4) INTEGER SC,I C DOUBLE PRECISION DMACH(5) C EQUIVALENCE (DMACH(1),SMALL(1)) EQUIVALENCE (DMACH(2),LARGE(1)) EQUIVALENCE (DMACH(3),RIGHT(1)) EQUIVALENCE (DMACH(4),DIVER(1)) EQUIVALENCE (DMACH(5),LOG10(1)) C C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T C 3B SERIES AND MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T C PC 7300), IN WHICH THE MOST SIGNIFICANT BYTE IS STORED FIRST. C DATA SMALL(1),SMALL(2) / 1048576, 0 / DATA LARGE(1),LARGE(2) / 2146435071, -1 / DATA RIGHT(1),RIGHT(2) / 1017118720, 0 / DATA DIVER(1),DIVER(2) / 1018167296, 0 / DATA LOG10(1),LOG10(2) / 1070810131, 1352628735 /, SC/987/ C C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES AND 8087-BASED C MICROS, SUCH AS THE IBM PC AND AT&T 6300, IN WHICH THE LEAST C SIGNIFICANT BYTE IS STORED FIRST. C C DATA SMALL(1),SMALL(2) / 0, 1048576 / C DATA LARGE(1),LARGE(2) / -1, 2146435071 / C DATA RIGHT(1),RIGHT(2) / 0, 1017118720 / C DATA DIVER(1),DIVER(2) / 0, 1018167296 / C DATA LOG10(1),LOG10(2) / 1352628735, 1070810131 /, SC/987/ C C MACHINE CONSTANTS FOR AMDAHL MACHINES. C C DATA SMALL(1),SMALL(2) / 1048576, 0 / C DATA LARGE(1),LARGE(2) / 2147483647, -1 / C DATA RIGHT(1),RIGHT(2) / 856686592, 0 / C DATA DIVER(1),DIVER(2) / 873463808, 0 / C DATA LOG10(1),LOG10(2) / 1091781651, 1352628735 /, SC/987/ C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA SMALL(1) / ZC00800000 / C DATA SMALL(2) / Z000000000 / C C DATA LARGE(1) / ZDFFFFFFFF / C DATA LARGE(2) / ZFFFFFFFFF / C C DATA RIGHT(1) / ZCC5800000 / C DATA RIGHT(2) / Z000000000 / C C DATA DIVER(1) / ZCC6800000 / C DATA DIVER(2) / Z000000000 / C C DATA LOG10(1) / ZD00E730E7 / C DATA LOG10(2) / ZC77800DC0 /, SC/987/ C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O0000000000000000 / C C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O0007777777777777 / C C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 /, SC/987/ C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O7770000000000000 / C C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O7777777777777777 / C C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 /, SC/987/ C C MACHINE CONSTANTS FOR FTN4 ON THE CDC 6000/7000 SERIES. C C DATA SMALL(1) / 00564000000000000000B / C DATA SMALL(2) / 00000000000000000000B / C C DATA LARGE(1) / 37757777777777777777B / C DATA LARGE(2) / 37157777777777777774B / C C DATA RIGHT(1) / 15624000000000000000B / C DATA RIGHT(2) / 00000000000000000000B / C C DATA DIVER(1) / 15634000000000000000B / C DATA DIVER(2) / 00000000000000000000B / C C DATA LOG10(1) / 17164642023241175717B / C DATA LOG10(2) / 16367571421742254654B /, SC/987/ C C MACHINE CONSTANTS FOR FTN5 ON THE CDC 6000/7000 SERIES. C C DATA SMALL(1) / O"00564000000000000000" / C DATA SMALL(2) / O"00000000000000000000" / C C DATA LARGE(1) / O"37757777777777777777" / C DATA LARGE(2) / O"37157777777777777774" / C C DATA RIGHT(1) / O"15624000000000000000" / C DATA RIGHT(2) / O"00000000000000000000" / C C DATA DIVER(1) / O"15634000000000000000" / C DATA DIVER(2) / O"00000000000000000000" / C C DATA LOG10(1) / O"17164642023241175717" / C DATA LOG10(2) / O"16367571421742254654" /, SC/987/ C C MACHINE CONSTANTS FOR CONVEX C-1 C C DATA SMALL(1),SMALL(2) / '00100000'X, '00000000'X / C DATA LARGE(1),LARGE(2) / '7FFFFFFF'X, 'FFFFFFFF'X / C DATA RIGHT(1),RIGHT(2) / '3CC00000'X, '00000000'X / C DATA DIVER(1),DIVER(2) / '3CD00000'X, '00000000'X / C DATA LOG10(1),LOG10(2) / '3FF34413'X, '509F79FF'X /, SC/987/ C C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. C C DATA SMALL(1) / 201354000000000000000B / C DATA SMALL(2) / 000000000000000000000B / C C DATA LARGE(1) / 577767777777777777777B / C DATA LARGE(2) / 000007777777777777776B / C C DATA RIGHT(1) / 376434000000000000000B / C DATA RIGHT(2) / 000000000000000000000B / C C DATA DIVER(1) / 376444000000000000000B / C DATA DIVER(2) / 000000000000000000000B / C C DATA LOG10(1) / 377774642023241175717B / C DATA LOG10(2) / 000007571421742254654B /, SC/987/ C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE - C STATIC DMACH(5) C C DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/ C DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/ C DATA LOG10/40423K,42023K,50237K,74776K/, SC/987/ C C MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7 C C DATA SMALL(1),SMALL(2) / '20000000, '00000201 / C DATA LARGE(1),LARGE(2) / '37777777, '37777577 / C DATA RIGHT(1),RIGHT(2) / '20000000, '00000333 / C DATA DIVER(1),DIVER(2) / '20000000, '00000334 / C DATA LOG10(1),LOG10(2) / '23210115, '10237777 /, SC/987/ C C MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. C C DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 / C DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 / C DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 / C DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 / C DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 /, SC/987/ C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86. C C DATA SMALL(1),SMALL(2) / Z00100000, Z00000000 / C DATA LARGE(1),LARGE(2) / Z7FFFFFFF, ZFFFFFFFF / C DATA RIGHT(1),RIGHT(2) / Z33100000, Z00000000 / C DATA DIVER(1),DIVER(2) / Z34100000, Z00000000 / C DATA LOG10(1),LOG10(2) / Z41134413, Z509F79FF /, SC/987/ C C MACHINE CONSTANTS FOR THE INTERDATA 8/32 C WITH THE UNIX SYSTEM FORTRAN 77 COMPILER. C C FOR THE INTERDATA FORTRAN VII COMPILER REPLACE C THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S. C C DATA SMALL(1),SMALL(2) / Z'00100000', Z'00000000' / C DATA LARGE(1),LARGE(2) / Z'7EFFFFFF', Z'FFFFFFFF' / C DATA RIGHT(1),RIGHT(2) / Z'33100000', Z'00000000' / C DATA DIVER(1),DIVER(2) / Z'34100000', Z'00000000' / C DATA LOG10(1),LOG10(2) / Z'41134413', Z'509F79FF' /, SC/987/ C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). C C DATA SMALL(1),SMALL(2) / "033400000000, "000000000000 / C DATA LARGE(1),LARGE(2) / "377777777777, "344777777777 / C DATA RIGHT(1),RIGHT(2) / "113400000000, "000000000000 / C DATA DIVER(1),DIVER(2) / "114400000000, "000000000000 / C DATA LOG10(1),LOG10(2) / "177464202324, "144117571776 /, SC/987/ C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). C C DATA SMALL(1),SMALL(2) / "000400000000, "000000000000 / C DATA LARGE(1),LARGE(2) / "377777777777, "377777777777 / C DATA RIGHT(1),RIGHT(2) / "103400000000, "000000000000 / C DATA DIVER(1),DIVER(2) / "104400000000, "000000000000 / C DATA LOG10(1),LOG10(2) / "177464202324, "047674776746 /, SC/987/ C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1),SMALL(2) / 8388608, 0 / C DATA LARGE(1),LARGE(2) / 2147483647, -1 / C DATA RIGHT(1),RIGHT(2) / 612368384, 0 / C DATA DIVER(1),DIVER(2) / 620756992, 0 / C DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 /, SC/987/ C C DATA SMALL(1),SMALL(2) / O00040000000, O00000000000 / C DATA LARGE(1),LARGE(2) / O17777777777, O37777777777 / C DATA RIGHT(1),RIGHT(2) / O04440000000, O00000000000 / C DATA DIVER(1),DIVER(2) / O04500000000, O00000000000 / C DATA LOG10(1),LOG10(2) / O07746420232, O20476747770 /, SC/987/ C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1),SMALL(2) / 128, 0 / C DATA SMALL(3),SMALL(4) / 0, 0 / C C DATA LARGE(1),LARGE(2) / 32767, -1 / C DATA LARGE(3),LARGE(4) / -1, -1 / C C DATA RIGHT(1),RIGHT(2) / 9344, 0 / C DATA RIGHT(3),RIGHT(4) / 0, 0 / C C DATA DIVER(1),DIVER(2) / 9472, 0 / C DATA DIVER(3),DIVER(4) / 0, 0 / C C DATA LOG10(1),LOG10(2) / 16282, 8346 / C DATA LOG10(3),LOG10(4) / -31493, -12296 /, SC/987/ C C DATA SMALL(1),SMALL(2) / O000200, O000000 / C DATA SMALL(3),SMALL(4) / O000000, O000000 / C C DATA LARGE(1),LARGE(2) / O077777, O177777 / C DATA LARGE(3),LARGE(4) / O177777, O177777 / C C DATA RIGHT(1),RIGHT(2) / O022200, O000000 / C DATA RIGHT(3),RIGHT(4) / O000000, O000000 / C C DATA DIVER(1),DIVER(2) / O022400, O000000 / C DATA DIVER(3),DIVER(4) / O000000, O000000 / C C DATA LOG10(1),LOG10(2) / O037632, O020232 / C DATA LOG10(3),LOG10(4) / O102373, O147770 /, SC/987/ C C MACHINE CONSTANTS FOR THE PRIME 50 SERIES SYSTEMS C WITH 32-BIT INTEGERS AND 64V MODE INSTRUCTIONS, C SUPPLIED BY IGOR BRAY. C C DATA SMALL(1),SMALL(2) / :10000000000, :00000100001 / C DATA LARGE(1),LARGE(2) / :17777777777, :37777677775 / C DATA RIGHT(1),RIGHT(2) / :10000000000, :00000000122 / C DATA DIVER(1),DIVER(2) / :10000000000, :00000000123 / C DATA LOG10(1),LOG10(2) / :11504046501, :07674600177 /, SC/987/ C C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000 C C DATA SMALL(1),SMALL(2) / $00000000, $00100000 / C DATA LARGE(1),LARGE(2) / $FFFFFFFF, $7FEFFFFF / C DATA RIGHT(1),RIGHT(2) / $00000000, $3CA00000 / C DATA DIVER(1),DIVER(2) / $00000000, $3CB00000 / C DATA LOG10(1),LOG10(2) / $509F79FF, $3FD34413 /, SC/987/ C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. C C DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 / C DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 / C DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 / C DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 / C DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 /, SC/987/ C C MACHINE CONSTANTS FOR THE VAX UNIX F77 COMPILER C C DATA SMALL(1),SMALL(2) / 128, 0 / C DATA LARGE(1),LARGE(2) / -32769, -1 / C DATA RIGHT(1),RIGHT(2) / 9344, 0 / C DATA DIVER(1),DIVER(2) / 9472, 0 / C DATA LOG10(1),LOG10(2) / 546979738, -805796613 /, SC/987/ C C MACHINE CONSTANTS FOR THE VAX-11 WITH C FORTRAN IV-PLUS COMPILER C C DATA SMALL(1),SMALL(2) / Z00000080, Z00000000 / C DATA LARGE(1),LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1),RIGHT(2) / Z00002480, Z00000000 / C DATA DIVER(1),DIVER(2) / Z00002500, Z00000000 / C DATA LOG10(1),LOG10(2) / Z209A3F9A, ZCFF884FB /, SC/987/ C C MACHINE CONSTANTS FOR VAX/VMS VERSION 2.2 C C DATA SMALL(1),SMALL(2) / '80'X, '0'X / C DATA LARGE(1),LARGE(2) / 'FFFF7FFF'X, 'FFFFFFFF'X / C DATA RIGHT(1),RIGHT(2) / '2480'X, '0'X / C DATA DIVER(1),DIVER(2) / '2500'X, '0'X / C DATA LOG10(1),LOG10(2) / '209A3F9A'X, 'CFF884FB'X /, SC/987/ C C *** ISSUE STOP 779 IF ALL DATA STATEMENTS ARE COMMENTED... IF (SC .NE. 987) STOP 779 IF (I .LT. 1 .OR. I .GT. 5) GOTO 999 D1MACH = DMACH(I) RETURN 999 WRITE(*,1999) I 1999 FORMAT(' D1MACH - I OUT OF BOUNDS',I10) STOP END CUT HERE............ cat > stesti.f <<'CUT HERE............' C C The following program will approximate three integrals over the C interval (0,infinity) and print out the values of the estimated C integrals, the estimated errors, the number of function evaluations C and IFAIL. C PROGRAM STESTI C C This demo-program has calls to the main driver: SQAINF C Through this main driver the rest of the modules are activated: C C level 1: SQAINF C C level 2: SCHINF and SADINF C C level 3: SGAINF, SGAINT, SEXINF, STRINT C C level 4: R1MACH (from blas in netlib) and FUNSUB (user specified) C EXTERNAL FUNSUB INTEGER NUMFUN,NW,MINPTS,MAXPTS,KEY,RESTAR INTEGER WRKSUB,EMAX,NUMNUL,NUM C C For the user's convenience we give a value to several parameters C through the parameter statement. C The following two parameters the user should not change, C unless the user changes the code appropriately: C PARAMETER (NUM=21,NUMNUL=16) C C The following three parameters the user may choose and these have C to remain unchanged through the entire run: C NUMFUN: the number of functions you want to integrate; C WRKSUB: the maximum number of subregions allowed; C EMAX: the maximum number of extrapolation steps; C PARAMETER (NUMFUN=3,WRKSUB=100,EMAX=25) C C Based on the previous 5 parameters we can now compute the C size of the main working space: NW. C PARAMETER (NW=1+5*NUMFUN+4*WRKSUB*NUMFUN+3*WRKSUB+(EMAX+1)**2+ +(NUMNUL+1+NUM)*NUMFUN) C C Next we choose the method to be used: C Three options: KEY = 0(Euler), 1(Euler modified), 2(Overholt's C P-order 2) C PARAMETER (KEY=2) C C We need to choose WRKSUB such that: C WRKSUB >= MAXPTS/NUM (We have NUM=21 fixed in this implementation) C This simply means that we have enough space to achieve MAXPTS C function evaluations. By choosing WRKSUB bigger than C this lower bound we may later change MAXPTS and RESTAR and restart C the computations from the point where the last run stopped. C The reason for stopping may have been that MAXPTS was too small C to achieve the requested error. C With our choice WRKSUB = 100 then any value of MAXPTS C between MINPTS and 2100 will be accepted by the code. Choosing C MAXPTS = 1000 then allows us to restart with a greater value C if necessary. C Note: For a given problem it may happen that the code stops after C reaching say MAXPTS=2100 function evaluations, having used less than C WRKSUB = 100 subregions. The reason for this is that if a given C region needs to be further subdivided then the old region will not C be stored while the function values computed over that region C will be counted. C The example gave the following output on C a SUN Sparc station 10: C C A= 0. B= 3.00000 C GAMMA= 1.00000 PERIOD= 6.28319 C MAXPTS= 300 WRKSUB= 100 NW= 2306 C KEY= 2 EPSREL= 1.00000E-04 C C Result and error: problem no. 1, 2 and 3: C 0.12431083E+01 0.71825701E+00 0.48923481E+00 C 0.25305973E-04 0.27054797E-04 0.26731776E-04 C NEVAL= 126 IFAIL= 0 C REAL A(WRKSUB),B(WRKSUB),EPSABS,EPSREL,GAMMA,PERIOD, +RESULT(NUMFUN),ABSERR(NUMFUN),WORK(NW),PI INTEGER NEVAL,IFAIL,IWORK(3*WRKSUB) LOGICAL DONE(NUMFUN) PI=2*ASIN(1.E0) C C We specify the period and the decay parameter GAMMA>0. If C the value of GAMMA < 0 then the code will attempt to estimate C GAMMA in the case KEY = 2. GAMMA is not used by the other two C methods. C PERIOD=2*PI GAMMA=1.0E0 C C Next we give the left endpoint of integration and the point b C from where we expect the theory to be valid. Advice: Since b C is not uniquely defined please do some experiments with C different values of b to gain experience. C The effect one achieves by increasing b depends on the value of KEY. C A(1)=0 B(1)=3 C C This is a first attempt: set RESTAR, give MINPTS, MAXPTS and the C absolute and relative error requests. C RESTAR=0 MINPTS=42 MAXPTS=300 EPSABS=0 EPSREL=1.E-4 C C Note: these four parameters and RESTAR may be changed in a restart run C WRITE (*,*) 'A=',A(1),' B=',B(1) WRITE (*,*) 'GAMMA=',GAMMA,' PERIOD=',PERIOD WRITE (*,*) 'MAXPTS=',MAXPTS,' WRKSUB=',WRKSUB,' NW=',NW WRITE (*,*) 'KEY=',KEY,' EPSREL=',EPSREL CALL SQAINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,WRKSUB,EMAX,MINPTS, +MAXPTS,EPSABS,EPSREL,NW,KEY,RESTAR,RESULT,ABSERR,NEVAL,IFAIL,DONE, +WORK,IWORK) WRITE (*,*) WRITE (*,*) 'Result and error: problem no. 1, 2 and 3:' WRITE (*,10) RESULT WRITE (*,10) ABSERR WRITE (*,*) 'NEVAL=',NEVAL,' IFAIL=',IFAIL 10 FORMAT (' ',3E16.8) END SUBROUTINE FUNSUB (X,NUMFUN,NP,FUNVLS) INTEGER NUMFUN,I,NP,J REAL X(NP),FUNVLS(NP,1),Y C C This subroutine must be provided by the user for computing C all components of the integrand at the given C evaluation points. C It must have parameters (X,NUMFUN,NP,FUNVLS) C Input parameters: C X Real array of dimension NP (= 21) defining the C x-coordinates of the evaluation points. C NUMFUN Integer that defines the number of C components of the integral I. C NP Integer that gives the number of evaluation points C in the quadrature rule used (Gauss, 21 point rule). C Output parameter: C FUNVLS Real array of dimension (NP,NUMFUN) that defines C the function values in the 21 evaluation points C for the NUMFUN components of the integrand. C C Note that we may save computer time when integrating C several functions simultaneously over the same interval C if we take advantage of the functions' similarities. C DO 10 I=1,NP Y=(2*SIN(X(I))+X(I)*COS(X(I))) DO 10 J=1,NUMFUN FUNVLS(I,J)=Y/(J+X(I)**2) 10 CONTINUE RETURN END CUT HERE............ cat > sqainf.f <<'CUT HERE............' SUBROUTINE SQAINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,WRKSUB,EMAX, +MINPTS,MAXPTS,EPSABS,EPSREL,NW,KEY,RESTAR,RESULT,ABSERR,NEVAL, +IFAIL,DONE,WORK,IWORK) C***BEGIN PROLOGUE SQAINF C***DATE WRITTEN 930515 (YYMMDD) C***REVISION DATE 940310 (YYMMDD) C***CATEGORY NO. H2B2A1/H2B2A2 C***AUTHORS C Terje O. Espelid, Department of Informatics, C University of Bergen, Hoyteknologisenteret. C N-5020 Bergen, Norway C Email.. terje@ii.uib.no C C Kjell J. Overholt, Department of Informatics, C University of Bergen, Hoyteknologisenteret. C N-5020 Bergen, Norway C Email.. kjell@ii.uib.no C C***PURPOSE To compute several one-dimensional integrals over C an infinite region. All of the integrands are assumed to be C oscillatory, with the same asymptotic periodic behavior. C This routine is the driving routine for the integration C routine SADINF. C C***KEYWORDS Quadrature, one dimensional, oscillatory integrands, C infinite region, extrapolation, globally adaptive, C periodic functions. C C***DESCRIPTION C C The routine calculates an approximation to a given vector of C infinite integrals C C I (F ,F ,...,F ) DX, C 1 2 NUMFUN C C where the interval of integration is [A,infinity), where A is C supplied by the user, and with C B as a possible subdivision point also supplied by the user: C A <= B. We assume that all functions have the same oscillatory C behavior for values of X => B. C C Here F = F (X), J = 1,2,...,NUMFUN. C J J C C A globally adaptive strategy, combined with extrapolation, C is applied in order to compute approximations RESULT(J), C for each component J of I, the following claim for accuracy: C C ABS(I(J)-RESULT(J)).LE.MAX(EPSABS,EPSREL*ABS(I(J))) C C SQAINF is a driver for the integration routine SADINF. C C SADINF either (1) repeatedly subdivides the interval with greatest C estimated errors; then estimates the integrals and the errors over C the new sub-intervals and finally updates the extrapolation tableau, C until the error request is met or MAXPTS function C evaluations have been used, or (2) decides instead to compute C a new term: the integral over the next half-period and then C perform a new extrapolation step. C C Only one basic integration rule (21 points) is offered: Gauss rule. C C Three different extrapolation methods are offered through the C parameter KEY:: Euler's method, a modification of Eulers method and C Overholt's P-order 2-method. C C If the values of the input parameters EPSABS or EPSREL are selected C great enough, the chosen integration rule is applied over the two C first intervals [A,B] and [B,B+PERIOD/2] for each C function in the vector and then one extrapolation step is performed C to give approximations RESULT(J), J = 1, 2, ..., NUMFUN. C No further subdivisions and no C more new terms will then be computed. C C When SQAINF computes estimates to a vector of integrals, C then all components of the vector are given the same treatment. C That is, I(F ) and I(F ) for J not equal to K, are estimated with C J K C the same subdivision of the original interval of integration C and the same extrapolation scheme is applied. C For integrals with enough similarity, we may save time by applying C SQAINF to all integrands in one call. For integrals that varies C continuously as functions of some parameter, the estimates C produced by SQAINF will also vary continuously when the same C subdivision is applied to all components. This will generally not C be the case when the different components are given separate C treatment. C C On the other hand this feature should be used with caution when the C different components of the integrals require clearly different C subdivisions. C C ON ENTRY C C NUMFUN Integer. C Number of components of the integral. C FUNSUB Externally declared subroutine for computing C all components of the integrand at the given C evaluation point. C It must have parameters (X,NUMFUN,NP,FUNVLS) C Input parameters: C X Real array of dimension 21 defining the C x-coordinates of the evaluation points. C NUMFUN Integer that defines the number of C components of the integral I. C NP Integer that gives the number of evaluation points C Output parameter: C FUNVLS Real array of dimension (21,NUMFUN) that defines C the function values in the 21 evaluation points C for the NUMFUN components of the integrand. C PERIOD Real. C All functions are assumed to have the same asymptotic C period > 0. C GAMMA Real. C All functions are assumed to decay as c/x**GAMMA, C when x >> 1 and we go in steps of length PERIOD, C (GAMMA > 0). C A Real array of dimension WRKSUB. C B Real array of dimension WRKSUB. C A(1) and B(1) are the left endpoint and right endpoint C of the first interval, A(1) <= B(1). B(1) is chosen by the C user and it is assumed that the integrand oscillates for C X >= B(1). Asymptotic period is PERIOD. Thus oscillations C are expected to be observed for points of distance PERIOD/2 C A(2), B(2), etc. are defined by the code, and as a warning: C the code may change the value of A(1) and B(1). C WRKSUB Integer. C Defines the length of the arrays A and B. (And thus the max C number of subregions allowed.) C EMAX Integer. C The maximum number of extrapolation steps. C MINPTS Integer. C Minimum number of function evaluations. C MAXPTS Integer. C Maximum number of function evaluations. C The number of function evaluations over each sub-interval C is 21. C EPSABS Real. C Requested absolute error. C EPSREL Real. C Requested relative error. C NW Integer. C Defines the length of the working array WORK. C We let C KEY Integer C Choice of extrapolation method: C KEY = 0 : Euler's method. C KEY = 1 : Modified Euler. C KEY = 2 : Overholt's polynomial order 2-method. C This last method is the only one that needs C the value of GAMMA. C RESTAR Integer. C If RESTAR = 0, this is the first attempt to compute C the integral. C If RESTAR = 1, C then we restart a previous attempt. C In this case the only parameters for SQAINF that may C be changed (with respect to the previous call of SQAINF) C are MINPTS, MAXPTS, EPSABS, EPSREL and RESTAR. C C ON RETURN C C RESULT Real array of dimension NUMFUN. C Approximations to all components of the integral. C ABSERR Real array of dimension NUMFUN. C Estimates of absolute errors. C NEVAL Integer. C Number of function evaluations used by SQAINF. C IFAIL Integer. C IFAIL = 0 for normal exit, when C C ABSERR(J) <= EPSABS or C ABSERR(J) <= ABS(RESULT(J))*EPSREL with MAXPTS or less C function evaluations for all values of J, C from 1 to NUMFUN. C C IFAIL = +1 if MAXPTS was too small C to obtain the required accuracy. In this case SQAINF C returns values of RESULT with estimated absolute C errors ABSERR. C IFAIL = -1 if EMAX was too small C to obtain the required accuracy. In this case SQAINF C returns values of RESULT with estimated absolute C errors ABSERR. C IFAIL = 2 if NUMFUN is less than 1. C IFAIL = 3 if B(1) < A(1). C IFAIL = 4 if unlegal PERIOD. C IFAIL = 5 if MAXPTS is less than MINPTS or MINPTS < 21. C IFAIL = 6 if EPSABS <= 0 and EPSREL <= 0. C IFAIL = 7 if WRKSUB is too small. C IFAIL = 8 if unlegal value of EMAX. C IFAIL = 9 if unlegal RESTAR. C IFAIL = 10 if unlegal value of KEY. C IFAIL = 11 if NW is less than LIMIT (See the code DCHINF). C IFAIL = 12 if either PERIOD is wrong or B(1) is too small. C IFAIL = 13 if unable to estimate GAMMA (In case KEY=2 only) C A Real array of dimension WRKSUB. C B Real array of dimension WRKSUB. C On exit A and B contains the endpoints of the intervals C used to produce the approximations to the integrals. C DONE Logical array of dimension NUMFUN. C On exit : DONE(I)=.TRUE. if function number I has been C integrated to the specified precision, else DONE(I)=.FALSE. C (Note that IFAIL = 1 just tells you C that something is wrong, however most of the integrals in C the vector may have been computed to the specified precisio C WORK Real array of dimension NW. C Used as working storage. C WORK(NW) = NSUB, the number of sub-intervals in the data C structure. C Let NW >= 1+5*NUMFUN+4*WRKSUB*NUMFUN C +3*WRKSUB+(EMAX+1)**2 + (NUMNUL+3+NUM)*NUMFUN) C Then C WORK(1),...,WORK(NUMFUN*WRKSUB) contain the estimated C components of the integrals over the sub-intervals. C WORK(NUMFUN*WRKSUB+1),...,WORK(2*NUMFUN*WRKSUB) contain C the estimated errors over the sub-intervals. C WORK(2*NUMFUN*WRKSUB+1),...,WORK(2*NUMFUN* C WRKSUB+WRKSUB) contain the greatest errors C in each sub-interval. C WORK(2*WRKSUB*NUMFUN+WRKSUB+1),...,WORK(2*WRKSUB*NUMFUN C +2*WRKSUB) contain the left sub-division point in each C sub-interval. C WORK(2*WRKSUB*NUMFUN+2*WRKSUB+1),...,WORK(2*WRKSUB*NUMFUN C +3*WRKSUB) contain the right sub-division point in each C sub-interval. C WORK((2*NUMFUN+3)*WRKSUB+1),...,WORK(NW) is used as C temporary storage in SADINF. C IWORK Integer array of dimension 3*WRKSUB. C Used as working storage. C C***LONG DESCRIPTION C C The information for each interval is contained in the data C structures A,B, WORK and IWORK. A and B contain the endpoints of C the intervals. When passed on to SADINF, WORK is split into four C arrays VALUES, ERRORS, GREATE and WORK2. VALUES contains the C estimated values of the integrals. ERRORS contains the estimated C errors. GREATE contains the greatest estimated error for each C interval. The data structures for the intervals are in SADINF C organized as a heap, and the size of GREATE(I) defines the C position of interval I in the heap. The heap is maintained by the C program STRINF and we use a partially ordered list of pointers, C saved in IWORK. When passed on to SADINF, IWORK is split into three C arrays WORST, LIST and UPOINT. LIST is a partially ordered list C such that GREATE(LIST(1)) is the maximum worst case error estimate C for all sub-intervals in our data-structure. In WORST the index to C the integral with the greatest error-estimate is kept. UPOINT is C an array of pointers to where in the U-sequence a region belongs. C This information is used when updating the corresponding U-term C after a subdivision. C C The subroutine for estimating the integral and the error over each C sub-interval, SGAINT, uses WORK2 as a work array. C C***REFERENCES C C T.O.Espelid and K.J.Overholt, DQAINF: An Algorithm for Automatic C Integration of Infinite Oscillating Tails, submitted for publ. 1993. C C***ROUTINES CALLED SCHINF,SADINF C***END PROLOGUE SQAINF C C Parameters C C C Global variables. C EXTERNAL FUNSUB INTEGER WRKSUB INTEGER NUMFUN,MINPTS,MAXPTS,NW,RESTAR,NEVAL,IFAIL,IWORK(3*WRKSUB) +,KEY,EMAX REAL A(WRKSUB),B(WRKSUB),EPSABS,EPSREL,RESULT(NUMFUN), +ABSERR(NUMFUN),WORK(NW),PERIOD,GAMMA LOGICAL DONE(NUMFUN) C C Local variables. C INTEGER NSUB,LENW,NUM,NUMNUL INTEGER I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12,I13,I14 PARAMETER (NUM=21,NUMNUL=16) C C***FIRST EXECUTABLE STATEMENT DQAINF C C Check the input parameters. C CALL SCHINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,MINPTS,MAXPTS,EPSABS, +EPSREL,WRKSUB,NW,EMAX,KEY,RESTAR,NEVAL,IFAIL) IF (IFAIL.NE.0) THEN RETURN END IF C C Split up the work space. C I1=1 I2=I1+WRKSUB*NUMFUN I3=I2+WRKSUB*NUMFUN I4=I3+WRKSUB I5=I4+WRKSUB I6=I5+WRKSUB I7=I6+WRKSUB*NUMFUN I8=I7+WRKSUB*NUMFUN I9=I8+NUMFUN I10=I9+NUMFUN I11=I10+NUMFUN I12=I11+(EMAX+1)*(EMAX+1) I13=I12+NUMFUN I14=I13+NUMFUN C C On restart runs the number of sub-intervals from the C previous call is assigned to NSUB. C IF (RESTAR.EQ.1) THEN NSUB=WORK(NW) END IF C C Compute the size of the temporary work space needed in DADINF. C LENW=(NUMNUL+1+NUM)*NUMFUN CALL SADINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,WRKSUB,EMAX,MINPTS, +MAXPTS,EPSABS,EPSREL,RESTAR,LENW,KEY,RESULT,ABSERR,NEVAL,NSUB, +IFAIL,DONE,WORK(I1),WORK(I2),WORK(I3),WORK(I4),WORK(I5),WORK(I6), +WORK(I7),WORK(I8),WORK(I9),WORK(I10),WORK(I11),WORK(I12),WORK(I13) +,WORK(I14),IWORK(1),IWORK(1+WRKSUB),IWORK(1+2*WRKSUB)) WORK(NW)=NSUB RETURN C C***END SQAINF C END CUT HERE............ cat > sadinf.f <<'CUT HERE............' SUBROUTINE SADINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,WRKSUB,EMAX, +MINPTS,MAXPTS,EPSABS,EPSREL,RESTAR,LENW,KEY,RESULT,ABSERR,NEVAL, +NSUB,IFAIL,DONE,VALUES,ERRORS,GREATE,C,D,U,E,T,CT,EXTERR,BETA, +UNEW,UERR,WORK2,WORST,LIST,UPOINT) C***BEGIN PROLOGUE SADINF C***KEYWORDS Quasi-linear extrapolation, infinite integration, C oscillating functions, adaptive strategy. C***PURPOSE To compute better estimates to a vector of approximations C to integrals of oscillating functions over an infinite C interval. C C***LAST MODIFICATION 94-03-10 C C***REFER TO SAQINF C C***DESCRIPTION Computation of integrals over an infinite interval C by combining extrapolation and adaptive strategies. C C SADINF uses extrapolation on a sum of integrals over subintervals C Each subinterval has length P = PERIOD/2 and the function involved i C assumed to be asymptotically periodic with period PERIOD. C This routine will EITHER take another extrapolation step C in order to reduce the pure extrapolation error, OR C subdivide the interval in the collection with greatest estimated C error. In both cases estimate(s) of the integral(s) and the error(s) C to the new sub-interval(s) are computed; this process continues C until the error request is met or MAXPTS function evaluations have C been used. C C The 21 point rule is: Gauss. C C ON ENTRY C C NUMFUN Integer. C The number of components of the integral. C FUNSUB Externally declared subroutine for computing C all components of the integrand at all NP C evaluation points. C It must have parameters (X,NUMFUN,FUNVLS) C Input parameters: C X The x-coordinates of the NP evaluation points. C NUMFUN Integer that defines the number of C components of I. C NP Integer that gives the number of evaluation points C Output parameter: C FUNVLS Real array of dimension NUMFUN C that defines NUMFUN components of the integrand. C C PERIOD Real. C All functions are assumed to have the same asymptotic C period. C GAMMA Real. C All functions are assumed to decay as c/x**GAMMA, C when x >> 1 and we go in steps of length PERIOD, C (GAMMA > 0). C A Real array of dimension WRKSUB. C B Real array of dimension WRKSUB. C A(1) is the left endpoint of integration. C B(1) >= A(1) is a user specified subdivision point. C We assume that the assumptions on which this code is C based hold for x >= B(1). We will not extrapolate based C on the function values between A(1) and B(1). C WRKSUB Integer. C Defines the length of the arrays A and B . C EMAX Integer. C The maximum number of extrapolation steps. C MINPTS Integer. C Minimum number of function evaluations. C MAXPTS C Maximum number of function evaluations. C EPSABS Real. C Requested absolute error. C EPSREL Real. C Requested relative error. C C RESTAR Integer. C If RESTAR = 0, this is the first attempt to compute C the integral. C If RESTAR = 1, C then we restart a previous attempt. C In this case the only parameters for SAQINF that may C be changed (with respect to the previous call of SAQINF) C are MINPTS, MAXPTS, EPSABS, EPSREL and RESTAR. C LENW Integer. C Length of the workspace WORK2. C KEY Integer C Choice of extrapolation method: C KEY = 0 : Euler's method. C KEY = 1 : Modified Euler. C KEY = 2 : Overholt's polynomial order 2-method. C This last method is the only one that needs C GAMMA. C C ON RETURN C C RESULT Real array of dimension NUMFUN. C Approximations to all components of the integral. C ABSERR Real array of dimension NUMFUN. C Estimates of absolute errors. C NEVAL Integer. C The number of function evaluations. C NSUB Integer. C The number of intervals in the data structure. C A Real array of dimension WRKSUB. C B Real array of dimension WRKSUB. C On exit A and B contains the endpoints of the intervals C used to produce the approximations to the integrals. C Warning: A(1) and B(1) may have changed due to C adaptive subdivision of the first interval. C IFAIL Integer. C IFAIL = 0 for normal exit. C IFAIL = +1 if MAXPTS was too small for SADINF C to obtain the required accuracy. In this case C SADINF returns values of RESULT with estimated C absolute errors ABSERR. C IFAIL = -1 if EMAX was too small for SADINF C to obtain the required accuracy. In this case C SADINF returns values of RESULT with estimated C absolute errors ABSERR. C C DONE Logical array of dimension NUMFUN. C On exit : DONE(I)=.TRUE. if function number I has been C integrated to the specified precision, else DONE(I)=.FALSE. C C VALUES Real array of dimension (NUMFUN,WRKSUB). C The estimated components of the integrals over the C sub-intervals. C ERRORS Real array of dimension (NUMFUN,WRKSUB). C The estimated errors over the sub-intervals. C GREATE Real array of dimension WRKSUB. C The greatest error in each sub-interval. C C Real array of dimension WRKSUB. C The left sub-division points of all intervals in the heap. C D Real array of dimension WRKSUB. C The right sub-division points of all intervals in the heap. C U Real array of dimension (NUMFUN,WRKSUB) C contains the terms in series. C E Real array of dimension (NUMFUN,WRKSUB) C contains the estimated errors in each U-term. C T Real array of dimension NUMFUN. C Dummy parameter. C CT Real array of dimension NUMFUN. C Dummy parameter. C EXTERR Real array of dimension NUMFUN. C These errors are associated with the region (LEFT,infinity) C and they are the pure extrapolation errors. C BETA Real array of dimension (0:EMAX,0:EMAX). C Dummy parameter. C UNEW Real array of dimension NUMFUN. C This gives the next terms in the series (new extrapolation C step) else it is the correction to the U-values (updating). C UERR Real array of dimension NUMFUN. C The estimated errors of all U-terms in the series. C WORK2 Real array of dimension LENW. C Work array used in SGAINT. C WORST Integer array of dimension WRKSUB. C Index to the greatest error in each sub-interval. C LIST Integer array used in STRINT of dimension WRKSUB. C Is a partially sorted list, where LIST(1) is the top C element in a heap of sub-intervals. C UPOINT Integer array of dimension WRKSUB. C Is an array of pointers to where in the U-sequence C a region belongs. This information is used when updating C the corresponding U-term after a subdivision. C C***ROUTINES CALLED STRINT, SGAINT, SEXINF, R1MACH. C***END PROLOGUE SADINF C C Global variables. C EXTERNAL FUNSUB,R1MACH INTEGER NUMFUN,LENW,RESTAR,NEVAL,WRKSUB,NSUB,IFAIL, +LIST(WRKSUB),KEY,WORST(WRKSUB),UPOINT(WRKSUB),EMAX, +MAXPTS,MINPTS REAL A(WRKSUB),B(WRKSUB),EPSABS,EPSREL,RESULT(NUMFUN), +ABSERR(NUMFUN),VALUES(NUMFUN,WRKSUB),PERIOD,ERRORS(NUMFUN,WRKSUB), +GREATE(WRKSUB),WORK2(LENW),C(WRKSUB),D(WRKSUB),U(NUMFUN,WRKSUB), +E(NUMFUN,WRKSUB),GAMMA,EXTERR(NUMFUN),T(NUMFUN),CT(NUMFUN), +BETA(0:EMAX,0:EMAX),UNEW(NUMFUN),UERR(NUMFUN),R1MACH LOGICAL DONE(NUMFUN) C C Local variables. C C SBRGNS is the number of stored sub-intervals. C NDIV The number of sub-intervals to be divided in each main step. C POINTR Pointer to the position in the data structure where C the new sub-intervals are to be stored. C G is the homogeneous coordinates of the integration rule. C W is the weights of the integration rule and the null rules. C TOP is a pointer to the top element in the heap of sub-intervals. C AOLD Left-endpoint of the interval to be processed. C BOLD Right-endpoint of the interval to be processed. C FLAG Signals small interval. C MORE Signals more work to do. C P Half a period. C NUMU The number of U-terms. C NUM The number of points in the basic rule. C NUMNUL The number nullrules used. C UPDATE Index to U-term to be updated. If UPDATE < 0, then this means C a new extrapolation step. C VACANT Index to vacant position in the data-structure. C EPMACH Machine epsilon C CC Ratio: the user specified subdivision point B(1)/PERIOD. INTEGER I,J,K,JJ,NUMINT,ONE,NUMU,VACANT,UPDATE INTEGER SBRGNS,I1,I2,I3,L1,L2,L3 INTEGER NDIV,POINTR,INDEX,TOP,FLAG,NUM,NUMNUL REAL AOLD,BOLD,GREAT,P,LEFT,MAXEER,EPMACH,CC LOGICAL MORE SAVE MAXEER,NUMU,LEFT,P,CC PARAMETER (NUMINT=2,NUM=21,NUMNUL=16) C C***FIRST EXECUTABLE STATEMENT SADINF C C Define machine epsilon. C EPMACH=R1MACH(4) C C Set some pointers for WORK2. C L1=1 L2=L1+NUMNUL*NUMFUN L3=L2+NUM*NUMFUN C IF (RESTAR.EQ.1) THEN SBRGNS=NSUB GO TO 90 ELSE FLAG=0 CC=B(1)/PERIOD END IF C C Initiate SBRGNS, DONE, MORE, P, A and B. C MORE=.TRUE. SBRGNS=0 DO 10 J=1,NUMFUN DONE(J)=.FALSE. 10 CONTINUE P=PERIOD/2 A(2)=B(1) B(2)=A(2)+P LEFT=B(2) C C Initialize E and U C DO 20 J=1,NUMFUN DO 20 I=1,WRKSUB E(J,I)=0 20 U(J,I)=0 C C Apply the rule procedure over the two first intervals. C IF (A(1).EQ.B(1)) THEN DO 30 I=1,NUMFUN VALUES(I,1)=0 30 ERRORS(I,1)=0 GREATE(1)=0 WORST(1)=1 ONE=2 SBRGNS=1 UPOINT(1)=1 ELSE ONE=1 END IF DO 40 I=ONE,2 CALL SGAINT (A(I),B(I),NUMFUN,FUNSUB,DONE,MORE,EPMACH, + WORK2(L1),WORK2(L2),WORK2(L3),FLAG,VALUES(1,I), + ERRORS(1,I),GREATE(I),WORST(I),C(I),D(I)) NEVAL=NEVAL+NUM SBRGNS=SBRGNS+1 UPOINT(I)=I 40 CONTINUE C C Initialize RESULT, U, E and NUMU. C DO 60 I=1,2 DO 50 J=1,NUMFUN U(J,I)=VALUES(J,I) E(J,I)=ERRORS(J,I) 50 CONTINUE 60 CONTINUE NUMU=2 C C Store results in heap. C DO 80 I=1,NUMINT GREAT=0.D0 DO 70 J=1,NUMFUN IF (ERRORS(J,I).GT.GREAT) THEN GREAT=ERRORS(J,I) WORST(I)=J END IF 70 CONTINUE GREATE(I)=GREAT C(I)=A(I)+(B(I)-A(I))/3 D(I)=B(I)-(B(I)-A(I))/3 INDEX=I CALL STRINT (2,INDEX,GREATE,LIST,I) 80 CONTINUE C C Extrapolate for the first time based on these 2 terms. C This will initialize the global error as well, and C MAXEER: the greatest extrapolation error. C Finally RESULT will be initialized. C UPDATE=-1 CALL SEXINF (NUMFUN,GAMMA,CC,KEY,NUMU-1,T,CT,UPDATE,EMAX,U,E, +RESULT,ABSERR,EXTERR,BETA,MAXEER) C C Check for termination. C IF (MORE.OR.(NEVAL.LT.MINPTS)) THEN GO TO 90 END IF IFAIL=0 GO TO 170 C C***End initialization. C C***Begin loop while the error is too great, C C First we check the ranking of the TOP element. C Then we determine if we will create a new term by integrating C from LEFT to LEFT + P and then extrapolate, or subdivide the TOP C element in the heap of subregions. C C 90 TOP=LIST(1) C C New term? C IF (GREATE(TOP).LT.(MAXEER)) THEN C C We want to extrapolate. Check if we C are allowed to take a new extrapolation step. C IF (NUMU-1.GT.EMAX) THEN IFAIL=-1 GO TO 170 ELSE VACANT=0 NDIV=1 END IF ELSE VACANT=TOP NDIV=3 END IF C C Check if NEVAL+NDIV*NUM is less than or equal to MAXPTS: C MAXPTS is the maximum number of function evaluations that are C allowed to be computed. C IF (NEVAL+NDIV*NUM.LE.MAXPTS) THEN C C When we pick an interval to divide in three, then one of the new C sub-intervals is stored in the original interval's position in the C data structure. C C Let POINTR point to the first free position in the data structure. C C Determine if this is a new extrapolation step or an update. C UPDATE will point to the element in the C U-series to be corrected or created. C IF (VACANT.EQ.0) THEN POINTR=SBRGNS+1 NUMU=NUMU+1 UPDATE=NUMU INDEX=POINTR A(INDEX)=LEFT LEFT=LEFT+P B(INDEX)=LEFT TOP=INDEX DO 100 J=1,NUMFUN UERR(J)=0 100 UNEW(J)=0 ELSE UPDATE=UPOINT(VACANT) C C Save the endpoints. C AOLD=A(TOP) BOLD=B(TOP) C C Remove the top element from the heap.(The value of K does not matter C POINTR=SBRGNS+2 CALL STRINT (1,SBRGNS,GREATE,LIST,K) C C Compute the three new intervals. C I1=TOP I2=POINTR-1 I3=POINTR A(I1)=AOLD B(I1)=C(TOP) A(I2)=B(I1) B(I2)=D(TOP) A(I3)=B(I2) B(I3)=BOLD INDEX=VACANT DO 110 J=1,NUMFUN UERR(J)=-ERRORS(J,INDEX) 110 UNEW(J)=-VALUES(J,INDEX) END IF C C Apply basic rule. C DO 130 I=1,NDIV CALL SGAINT (A(INDEX),B(INDEX),NUMFUN,FUNSUB,DONE,MORE, + EPMACH,WORK2(L1),WORK2(L2),WORK2(L3),FLAG, + VALUES(1,INDEX),ERRORS(1,INDEX),GREATE(INDEX), + WORST(INDEX),C(INDEX),D(INDEX)) UPOINT(INDEX)=UPDATE C C Compute the U-term C DO 120 J=1,NUMFUN UNEW(J)=UNEW(J)+VALUES(J,INDEX) UERR(J)=UERR(J)+ERRORS(J,INDEX) 120 CONTINUE INDEX=SBRGNS+I+1 130 CONTINUE NEVAL=NEVAL+NDIV*NUM C C Compute the E and U terms (These may be new terms or terms that C have to be updated), C DO 140 J=1,NUMFUN U(J,UPDATE)=U(J,UPDATE)+UNEW(J) E(J,UPDATE)=E(J,UPDATE)+UERR(J) 140 CONTINUE C C Do the extrapolation and compute the global results and errors. C UPDATE is used to signal an extrapolation step. C IF (VACANT.EQ.0) THEN UPDATE=-1 END IF CALL SEXINF (NUMFUN,GAMMA,CC,KEY,NUMU-1,T,CT,UPDATE,EMAX,U,E, + RESULT,ABSERR,EXTERR,BETA,MAXEER) C C Check each integration problem and set DONE(J)=.TRUE. if they C are finished. MORE will signal if there is more to do. C MORE=.FALSE. DO 150 J=1,NUMFUN IF (ABSERR(J).LE.EPSREL*ABS(RESULT(J)).OR.ABSERR(J).LE. + EPSABS) THEN DONE(J)=.TRUE. ELSE MORE=.TRUE. DONE(J)=.FALSE. END IF 150 CONTINUE C C Store results in heap. C INDEX=TOP DO 160 I=1,NDIV JJ=SBRGNS+I CALL STRINT (2,JJ,GREATE,LIST,INDEX) INDEX=SBRGNS+I+1 160 CONTINUE SBRGNS=SBRGNS+NDIV C C Check for termination. C IF (MORE.OR.(NEVAL.LT.MINPTS)) THEN GO TO 90 END IF C C Else we did not succeed with the C given value of MAXPTS. Note: We do not use IFAIL to signal C FLAG in the current implementation. C ELSE IFAIL=+1 END IF C 170 NSUB=SBRGNS RETURN C C***END SADINF C END CUT HERE............ cat > schinf.f <<'CUT HERE............' SUBROUTINE SCHINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,MINPTS,MAXPTS, +EPSABS,EPSREL,WRKSUB,NW,EMAX,KEY,RESTAR,NEVAL,IFAIL) C***BEGIN PROLOGUE SCHINF C***REFER TO SQAINF C***PURPOSE SCHINF checks the validity of the input parameters to C SQAINF. C C***LAST MODIFICATION 93-05-05 C C***DESCRIPTION C SCHINF computes IFAIL as a function of the input C parameters to SQAINF, and checks the validity of the input C parameters to SQAINF. C C ON ENTRY C C NUMFUN Integer. C Number of components of the integral. C FUNSUB Externally declared subroutine for computing C all components of the integrand at the given C evaluation point. C It must have parameters (X,NUMFUN,FUNVLS) C Input parameters: C X Real array of dimension 21 defining the C x-coordinates of the evaluation points. C NUMFUN Integer that defines the number of C components of I. C NP Integer that gives the number of evaluation points C Output parameter: C FUNVLS Real array of dimension (21,NUMFUN) that defines C the function values in the 21 evaluation points C for the NUMFUN components of the integrand. C PERIOD Real. C All functions are assumed to have the same asymptotic C period. C GAMMA Real. C All functions are assumed to decay as c/x**GAMMA, C when x >> 1 and we go in steps of length PERIOD, C (GAMMA > 0). C A Real array of dimension WRKSUB. C B Real array of dimension WRKSUB. C A(1) and B(1) are the left endpoint and right endpoint C of the first interval, A(1) <= B(1). B(1) is chosen by the C user and it is assumed that the integrand oscillates for C X >= B(1). Asymptotic period is PERIOD. Thus oscillations C are expected to be observed for points of distance PERIOD/2 C MINPTS Integer. C Minimum number of function evaluations. C MAXPTS Integer. C Maximum number of function evaluations. C The number of function evaluations over each sub-interval C is NUM = 21. C EPSABS Real. C Requested absolute error. C EPSREL Real. C Requested relative error. C WRKSUB Integer. C Defines the length of the arrays A and B. (And thus the C maximum number of subregions allowed.) C NW Integer. C Defines the length of the working array WORK. C EMAX Integer. C The maximum number of extrapolation steps. (At least one C step!) C KEY Integer. C Choice of extrapolation method: C KEY = 0 : Euler's method. C KEY = 1 : Modified Euler. C KEY = 2 : Overholt's polynomial order 2-method. C This last method is the only one that needs C GAMMA. C RESTAR Integer. C If RESTAR = 0, this is the first attempt to compute C the integral. C If RESTAR = 1, C then we restart a previous attempt. C RESTAR is allowed to take these two values only. C ON RETURN C C NEVAL Integer. C The number of function evaluations. C IFAIL Integer. C IFAIL = 0 for normal exit. C IFAIL = 2 if NUMFUN is less than 1. C IFAIL = 3 if B(1) < A(1). C IFAIL = 4 if unlegal PERIOD. C IFAIL = 5 if MAXPTS is less than MINPTS or MINPTS < 21. C IFAIL = 6 if EPSABS <= 0 and EPSREL <= 0. C IFAIL = 7 if WRKSUB is too small. C IFAIL = 8 if unlegal value of EMAX. C IFAIL = 9 if unlegal RESTAR. C IFAIL = 10 if unlegal value of KEY. C IFAIL = 11 if NW is less than LIMIT (See the code SCHINF). C IFAIL = 12 if either PERIOD is wrong or B(1) is too small. C IFAIL = 13 if unable to estimate GAMMA (In case KEY=2 only) C C***ROUTINES CALLED SGAINF C***END PROLOGUE SCHINF C C Global variables. C INTEGER NUMFUN,MINPTS,MAXPTS,NW,RESTAR,EMAX,NEVAL,WRKSUB,KEY, +IFAIL REAL A(WRKSUB),B(WRKSUB),EPSABS,EPSREL,PERIOD,GAMMA EXTERNAL SGAINF,FUNSUB C C Local variables. C INTEGER LIMIT,NUMNUL,NUM REAL WIDTH PARAMETER (NUMNUL=16,NUM=21) C C***FIRST EXECUTABLE STATEMENT SCHINF C IFAIL=0 C C Check on legal NUMFUN. C IF (NUMFUN.LT.1) THEN IFAIL=2 RETURN END IF C C Check on legal length of the first interval of integration. C WIDTH=B(1)-A(1) IF (WIDTH.LT.0) THEN IFAIL=3 RETURN END IF C C Check on legal sign of PERIOD. C IF (PERIOD.LE.0) THEN IFAIL=4 RETURN END IF C C Check on MAXPTS >= MINPTS. C IF ((MAXPTS.LT.MINPTS).OR.(MINPTS.LT.NUM)) THEN IFAIL=5 RETURN END IF C C Check on legal accuracy requests. C IF (EPSABS.LE.0.AND.EPSREL.LE.0) THEN IFAIL=6 RETURN END IF C C Check on legal WRKSUB. C IF (NUM*WRKSUB.LT.MAXPTS) THEN IFAIL=7 RETURN END IF C C Check on legal value of EMAX. C IF (EMAX.LT.1) THEN IFAIL=8 RETURN END IF C C Check on legal RESTAR. C IF (RESTAR.NE.0.AND.RESTAR.NE.1) THEN IFAIL=9 RETURN ELSE IF (RESTAR.EQ.0) THEN NEVAL=0 END IF C C Check on legal KEY. C IF (KEY.NE.0.AND.KEY.NE.1.AND.KEY.NE.2) THEN IFAIL=10 RETURN END IF C C Check on big enough NW. C LIMIT=1+5*NUMFUN+4*WRKSUB*NUMFUN+3*WRKSUB+(EMAX+1)**2+(NUMNUL+1+ +NUM)*NUMFUN IF (NW.LT.LIMIT) THEN IFAIL=11 RETURN END IF C C If the user gives GAMMA <= 0 and KEY = 2 then SQAINF will try to C estimate GAMMA. C If the subroutine estimates GAMMA successfully within its error C bound, IFAIL will remain 0 while GAMMA is given the new value. C The routine SGAINF may detect errors: C C IFAIL = 12 indicates that either C B is too small or the PERIOD is wrong, or maybe the function does no C have the desired properties. (NOTE: The code checks ONLY the first C function in the vector, assuming that if this is correct then the C other functions share its properties.) C C IFAIL = 13 : No errors where detected, but the code was unable to C compute GAMMA sufficiently accurate. Remedy: use KEY = 1 instead. C IF ((GAMMA.LE.0).AND.(KEY.EQ.2)) THEN CALL SGAINF (NUMFUN,FUNSUB,B(1),PERIOD,GAMMA,NEVAL,IFAIL) RETURN END IF C RETURN C C***END SCHINF C END CUT HERE............ cat > strint.f <<'CUT HERE............' SUBROUTINE STRINT (DVFLAG,SBRGNS,GREATE,LIST,NEW) C***BEGIN PROLOGUE STRINT C***REFER TO SQAINF C***PURPOSE STRINT maintains a heap of subregions. C***DESCRIPTION STRINT maintains a heap of subregions. The subregions C are stored in a partially sorted binary tree, ordered to the C size of the greatest error estimates of each subregion(GREATE). C The subregion with greatest error estimate is in the first C position of the heap. C C PARAMETERS C C DVFLAG Integer. C If DVFLAG = 1, we remove the subregion with C greatest error from the heap. C If DVFLAG = 2, we insert a new subregion in the heap. C If DVFLAG = 3, we move the top element to a new position. C SBRGNS Integer. C Number of subregions in the heap. C GREATE Real array of dimension SBRGNS. C Used to store the greatest estimated errors in C all subregions. C LIST Integer array of dimension SBRGNS. C Used as a partially ordered list of pointers to the C different subregions. This list is a heap where the C element on top of the list is the subregion with the C greatest error estimate. C NEW Integer. C Index to the new region to be inserted in the heap. C***ROUTINES CALLED-NONE C***END PROLOGUE STRINT C C Global variables. C INTEGER DVFLAG,NEW,SBRGNS,LIST(*) REAL GREATE(*) C C Local variables. C C GREAT is used as intermediate storage for the greatest error of a C subregion. C PNTR Pointer. C SUBRGN Position of child/parent subregion in the heap. C SUBTMP Position of parent/child subregion in the heap. INTEGER SUBRGN,SUBTMP,PNTR REAL GREAT C C***FIRST EXECUTABLE STATEMENT STRINT C C If DVFLAG = 1, we will reduce the partial ordered list by the C element with greatest estimated error. Thus the element in C in the heap with index LIST(1) is vacant and can be used later. C Reducing the heap by one element implies that the last element C should be re-positioned. C If DVFLAG = 3, we will keep the size of the partially ordered C list but re-position the first element. C IF (DVFLAG.NE.2) THEN IF (DVFLAG.EQ.1) THEN PNTR=LIST(SBRGNS) GREAT=GREATE(PNTR) SBRGNS=SBRGNS-1 ELSE IF (DVFLAG.EQ.3) THEN PNTR=LIST(1) GREAT=GREATE(PNTR) END IF SUBRGN=1 10 SUBTMP=2*SUBRGN IF (SUBTMP.LE.SBRGNS) THEN IF (SUBTMP.NE.SBRGNS) THEN C C Find max. of left and right child. C IF (GREATE(LIST(SUBTMP)).LT.GREATE(LIST(SUBTMP+1))) THEN SUBTMP=SUBTMP+1 END IF END IF C C Compare max.child with parent. C If parent is max., then done. C IF (GREAT.LT.GREATE(LIST(SUBTMP))) THEN C C Move the pointer at position subtmp up the heap. C LIST(SUBRGN)=LIST(SUBTMP) SUBRGN=SUBTMP GO TO 10 END IF END IF C C Update the pointer. C IF (SBRGNS.GT.0) THEN LIST(SUBRGN)=PNTR END IF ELSE IF (DVFLAG.EQ.2) THEN C C If DVFLAG = 2, find the position for the NEW region in the heap. C GREAT=GREATE(NEW) SUBRGN=SBRGNS 20 SUBTMP=SUBRGN/2 IF (SUBTMP.GE.1) THEN C C Compare max.child with parent. C If parent is max, then done. C IF (GREAT.GT.GREATE(LIST(SUBTMP))) THEN C C Move the pointer at position subtmp down the heap. C LIST(SUBRGN)=LIST(SUBTMP) SUBRGN=SUBTMP GO TO 20 END IF END IF C C Set the pointer to the new region in the heap. C LIST(SUBRGN)=NEW END IF C C***END STRINT C RETURN END CUT HERE............ cat > sgaint.f <<'CUT HERE............' SUBROUTINE SGAINT (A,B,NUMFUN,FUNSUB,DONE,MORE,EPMACH,NULL, +FUNVAL,BASABS,FLAG,BASVAL,RGNERR,GREATE,WORST,C,D) C***BEGIN PROLOGUE SGAINT C***PURPOSE To compute basic integration rule values and C corresponding error estimates. C C***LAST MODIFICATION 94-03-10 C C***REFER TO SQAINF C C***DESCRIPTION SGAINT computes basic integration rule values C for a vector of integrands over an interval. C SGAINT also computes estimates for the errors by C using several null rule approximations. C ON ENTRY C C A Real. C B Real. C The left and right endpoints of the interval. C NUMFUN Integer. C Number of components of the vector integrand. C FUNSUB Externally declared subroutine for computing C all components of the integrand at the given C evaluation point. C It must have parameters (X,NUMFUN,NP,FUNVLS) C Input parameters: C X The x-coordinates of all 21 evaluation points. C NUMFUN Integer that defines the number of C components of I. C NP Integer that gives the number of evaluation points C Output parameter: C FUNVLS Real array of dimension (21,NUMFUN) C that defines NUMFUN components of the integrand C evaluated in all 21 evaluation points. C C DONE Logical array of dimension NUMFUN. C DONE(I)=.TRUE. if function number I has been C integrated to the specified precision, else DONE(I)=.FALSE. C MORE Logical. C Information about the fact that there is still work to do. C EPMACH Real. C Machine epsilon. C NULL Real array of dimension (16,NUMFUN). C A work array. C FUNVAL Real array of dimension (21,NUMFUN). C A work array. C BASABS Real array of dimension NUMFUN. C A work array. C ON RETURN C C FLAG Integer. C Signals that at least one interval has become too small C to distinguish between a rule evaluation point and an C endpoint. C C BASVAL Real array of dimension NUMFUN. C The values for the basic rule for each component C of the integrand. C RGNERR Real array of dimension NUMFUN. C The error estimates for each component of the integrand. C GREATE Real. C The greatest error component of the integrand. C WORST Index to the greatest error in each sub-interval. C C C Real. C D Real C C and D are subdivision points based on C uniform trisection: C = A + (B-A)/3 and D = B - (B-a)/3. C C***REFERENCES Espelid, T. O., Integration Rules, Null Rules and C Error Estimation, Report no 33, Informatics, Univ. of Bergen, C 1988. C C Berntsen,J. and Espelid,T.O., Error estimation in Automatic C Quadrature Routines, ACM Trans. Math. Softw., C 17, 2, 233-252, 1991. C C Espelid, T. O., DQAINT: An Algorithm for Adaptive Quadrature C over a Collection of Finite Intervals, in Numerical C Integration, Recent Developments, Software and Applications, C Espelid and Genz (eds), NATO ASI Series C: C Mathematical and Physical Sciences - Vol. 357, C Kluwer, Dordrecht, The Netherlands, pages 341-342, 1992. C C***ROUTINES CALLED FUNSUB C***END PROLOGUE SGAINT C C Global variables. C EXTERNAL FUNSUB LOGICAL DONE(NUMFUN),MORE INTEGER NUMFUN,FLAG,WORST REAL A,B,BASVAL(NUMFUN),RGNERR(NUMFUN),NULL(16, +NUMFUN),GREATE,C,D,BASABS(NUMFUN),FUNVAL(21,NUMFUN),EPMACH C C Local variables. C C WTLENG The number of weights of the integration rule. C G Real array of dimension WTLENG. C The coordinates for the generators of C the evaluation points. C The integration rule is using symmetric evaluation C points and has the structure (1,6,3). That is, C 1 point of multiplicity 1, C 6 sets of points of multiplicity 3 and C 3 sets of points of multiplicity 6. C This gives totally 37 evaluation points. C In order to reduce the number of loops in SGAINT, C the 3 loops for the sets of multiplicity 6 are split C into 6 loops and added to the loops for the sets of C multiplicity 3. C The number of weights we have to give with C this splitting is 13(WTLENG). C C W Real array of dimension (9,WTLENG). C The weights of the basic rule and the null rules. C W(1),...,W(WTLENG) are weights for the basic rule. C C X Real array of dimension 21. C Evaluation points in (A,B). C HLEN Real. C Half of intervals length INTEGER I,I1,I2,J,K,WTLENG,DGE1,DEGREE,NUMNUL,DIMNUL PARAMETER (NUMNUL=8,DIMNUL=16) REAL X(21),RR(NUMNUL/2-1),R,NOISE,E(NUMNUL/2),EMAX, +ALPHA,CONST,W(11),ABSCIS,HLEN,SAFETY,RCRIT,ABSSUM,SUM,DIFF,CENTR, +G(11),FACTOR,NULLW(16,11) PARAMETER (WTLENG=11,FACTOR=0.02D+0,SAFETY=10.E0,RCRIT=0.5E0,DGE1= +18,DEGREE=41,ALPHA=(DEGREE-DGE1)/2.E0) C C The abscissas to the 21 point Gauss integration rule. C The abscissas are given in -1,1: due to the symmetry we only specify C values in 0,1. C DATA (G(I),I=1,11)/0.9937521706203895D0,0.9672268385663062D0,0. +9200993341504008D0,0.8533633645833172D0,0.7684399634756779D0,0. +6671388041974123D0,0.5516188358872198D0,0.4243421202074387D0,0. +2880213168024010D0,0.1455618541608950D0,0.000000000000000D0/ C C Weights of the 21 point Gauss quadrature rule. Same remark C about symmetry. C DATA (W(I),I=1,11)/0.01601722825777433D0,0.03695378977085249D0,0. +05713442542685720D0,0.07610011362837930D0,0.09344442345603386D0,0. +1087972991671483D0,0.1218314160537285D0,0.1322689386333374D0,0. +1398873947910731D0,0.1445244039899700D0,0.1460811336496904D0/ C C Weights of the 5 symmetric nullrules of degrees 19,17,15,13,11 C Weights of the 5 asymmetric nullrules of degrees 18,16,14,12,10 C Same remark about symmetry as above. The nullrules are all C orthogonal and has the same norm as the basic rule and are C given with decreasingly degrees. C DATA NULLW(1,1)/0.5859224694730026D-02/ DATA NULLW(1,2)/-0.2024707000741622D-01/ DATA NULLW(1,3)/0.3883580247121445D-01/ DATA NULLW(1,4)/-0.5965412595242497D-01/ DATA NULLW(1,5)/0.8114279498343020D-01/ DATA NULLW(1,6)/-0.1019231472030305D+00/ DATA NULLW(1,7)/0.1207652963052454D+00/ DATA NULLW(1,8)/-0.1366043796711610D+00/ DATA NULLW(1,9)/0.1485698262567817D+00/ DATA NULLW(1,10)/-0.1560150459859118D+00/ DATA NULLW(1,11)/0.1585416482170856D+00/ DATA NULLW(2,1)/0.1301976799706014D-01/ DATA NULLW(2,2)/-0.4379005851020851D-01/ DATA NULLW(2,3)/0.7990096087086482D-01/ DATA NULLW(2,4)/-0.1138307201442027D+00/ DATA NULLW(2,5)/0.1394263659262734D+00/ DATA NULLW(2,6)/-0.1520456605731098D+00/ DATA NULLW(2,7)/0.1489588260146731D+00/ DATA NULLW(2,8)/-0.1296181347851887D+00/ DATA NULLW(2,9)/0.9568420420614478D-01/ DATA NULLW(2,10)/-0.5078074459100106D-01/ DATA NULLW(2,11)/0.0000000000000000D+00/ DATA NULLW(3,1)/0.2158184987561927D-01/ DATA NULLW(3,2)/-0.6965227115767195D-01/ DATA NULLW(3,3)/0.1174438965053943D+00/ DATA NULLW(3,4)/-0.1473794001233916D+00/ DATA NULLW(3,5)/0.1481989091733945D+00/ DATA NULLW(3,6)/-0.1168273220215079D+00/ DATA NULLW(3,7)/0.5890214560028095D-01/ DATA NULLW(3,8)/0.1273585156460484D-01/ DATA NULLW(3,9)/-0.8133037350927629D-01/ DATA NULLW(3,10)/0.1304777802379009D+00/ DATA NULLW(3,11)/-0.1483021322906934D+00/ DATA NULLW(4,1)/0.3119657617222001D-01/ DATA NULLW(4,2)/-0.9516116459787523D-01/ DATA NULLW(4,3)/0.1431705707841666D+00/ DATA NULLW(4,4)/-0.1462171986707822D+00/ DATA NULLW(4,5)/0.9677919508585969D-01/ DATA NULLW(4,6)/-0.1075583592960879D-01/ DATA NULLW(4,7)/-0.7936141880173113D-01/ DATA NULLW(4,8)/0.1380749552472930D+00/ DATA NULLW(4,9)/-0.1417577117227090D+00/ DATA NULLW(4,10)/0.8867798725104829D-01/ DATA NULLW(4,11)/0.0000000000000000D+00/ DATA NULLW(5,1)/0.4157639307601386D-01/ DATA NULLW(5,2)/-0.1179114894020921D+00/ DATA NULLW(5,3)/0.1511566572815612D+00/ DATA NULLW(5,4)/-0.1073644825716617D+00/ DATA NULLW(5,5)/0.4174411212397235D-02/ DATA NULLW(5,6)/0.1012057375471417D+00/ DATA NULLW(5,7)/-0.1472858866418607D+00/ DATA NULLW(5,8)/0.1063772962797608D+00/ DATA NULLW(5,9)/-0.2323645744823691D-02/ DATA NULLW(5,10)/-0.1030910117645103D+00/ DATA NULLW(5,11)/0.1469720414561505D+00/ DATA NULLW(6,1)/0.5246625962340516D-01/ DATA NULLW(6,2)/-0.1358182960749584D+00/ DATA NULLW(6,3)/0.1386355746642898D+00/ DATA NULLW(6,4)/-0.3967547044517777D-01/ DATA NULLW(6,5)/-0.8983329656061084D-01/ DATA NULLW(6,6)/0.1471801958801420D+00/ DATA NULLW(6,7)/-0.8524048604745531D-01/ DATA NULLW(6,8)/-0.4617298114857220D-01/ DATA NULLW(6,9)/0.1397282757969823D+00/ DATA NULLW(6,10)/-0.1185867937861861D+00/ DATA NULLW(6,11)/0.0000000000000000D+00/ DATA NULLW(7,1)/0.6363031447247006D-01/ DATA NULLW(7,2)/-0.1472028208379138D+00/ DATA NULLW(7,3)/0.1063757528761779D+00/ DATA NULLW(7,4)/0.3881687506910375D-01/ DATA NULLW(7,5)/-0.1432999618142209D+00/ DATA NULLW(7,6)/0.9698969424297650D-01/ DATA NULLW(7,7)/0.5209450556829039D-01/ DATA NULLW(7,8)/-0.1455658951943161D+00/ DATA NULLW(7,9)/0.8343286549711612D-01/ DATA NULLW(7,10)/0.6800562635441401D-01/ DATA NULLW(7,11)/-0.1465539124681842D+00/ DATA NULLW(8,1)/0.7484599067063009D-01/ DATA NULLW(8,2)/-0.1508776892345244D+00/ DATA NULLW(8,3)/0.5853283458897407D-01/ DATA NULLW(8,4)/0.1062457136342151D+00/ DATA NULLW(8,5)/-0.1318732622123368D+00/ DATA NULLW(8,6)/-0.1673118490576078D-01/ DATA NULLW(8,7)/0.1428979325476485D+00/ DATA NULLW(8,8)/-0.7818432195969258D-01/ DATA NULLW(8,9)/-0.9112601052788798D-01/ DATA NULLW(8,10)/0.1382849496064090D+00/ DATA NULLW(8,11)/0.0000000000000000D+00/ DATA NULLW(9,1)/0.8590167508061712D-01/ DATA NULLW(9,2)/-0.1462121283895959D+00/ DATA NULLW(9,3)/0.1973066649848703D-02/ DATA NULLW(9,4)/0.1434120884358229D+00/ DATA NULLW(9,5)/-0.6050045998747565D-01/ DATA NULLW(9,6)/-0.1192968264851738D+00/ DATA NULLW(9,7)/0.1063582979588903D+00/ DATA NULLW(9,8)/0.7871971420591506D-01/ DATA NULLW(9,9)/-0.1360664606736734D+00/ DATA NULLW(9,10)/-0.2747426951466367D-01/ DATA NULLW(9,11)/0.1463706054390223D+00/ DATA NULLW(10,1)/0.9659618304508728D-01/ DATA NULLW(10,2)/-0.1331667766592828D+00/ DATA NULLW(10,3)/-0.5483608118503819D-01/ DATA NULLW(10,4)/0.1395396581193282D+00/ DATA NULLW(10,5)/0.3842219337177220D-01/ DATA NULLW(10,6)/-0.1430606163131484D+00/ DATA NULLW(10,7)/-0.2498840938693774D-01/ DATA NULLW(10,8)/0.1451753725543706D+00/ DATA NULLW(10,9)/0.1236834040097046D-01/ DATA NULLW(10,10)/-0.1461902970879641D+00/ DATA NULLW(10,11)/0.0000000000000000D+00/ DATA NULLW(11,1)/0.1067392105869384D+00/ DATA NULLW(11,2)/-0.1122928178247447D+00/ DATA NULLW(11,3)/-0.1031959020477783D+00/ DATA NULLW(11,4)/0.9558143541939021D-01/ DATA NULLW(11,5)/0.1196951864405313D+00/ DATA NULLW(11,6)/-0.7225984000378730D-01/ DATA NULLW(11,7)/-0.1339424732473705D+00/ DATA NULLW(11,8)/0.4492456833975673D-01/ DATA NULLW(11,9)/0.1431238351576778D+00/ DATA NULLW(11,10)/-0.1523606417516131D-01/ DATA NULLW(11,11)/-0.1462742772906911D+00/ DATA NULLW(12,1)/0.1161523191050730D+00/ DATA NULLW(12,2)/-0.8469391287377601D-01/ DATA NULLW(12,3)/-0.1355896186086413D+00/ DATA NULLW(12,4)/0.2408868272651161D-01/ DATA NULLW(12,5)/0.1460359069105494D+00/ DATA NULLW(12,6)/0.4632194803727984D-01/ DATA NULLW(12,7)/-0.1231814607655799D+00/ DATA NULLW(12,8)/-0.1068762140630544D+00/ DATA NULLW(12,9)/0.7029919038187181D-01/ DATA NULLW(12,10)/0.1416700606593806D+00/ DATA NULLW(12,11)/0.0000000000000000D+00/ DATA NULLW(13,1)/0.1246701955330830D+00/ DATA NULLW(13,2)/-0.5195253287985397D-01/ DATA NULLW(13,3)/-0.1469123277046623D+00/ DATA NULLW(13,4)/-0.5433985749387003D-01/ DATA NULLW(13,5)/0.1052913655334742D+00/ DATA NULLW(13,6)/0.1341759283651172D+00/ DATA NULLW(13,7)/-0.2310968036052471D-02/ DATA NULLW(13,8)/-0.1358135230198954D+00/ DATA NULLW(13,9)/-0.1024826424015055D+00/ DATA NULLW(13,10)/0.5656562071840163D-01/ DATA NULLW(13,11)/0.1462174827723311D+00/ DATA NULLW(14,1)/0.1321420280297885D+00/ DATA NULLW(14,2)/-0.1602500237425218D-01/ DATA NULLW(14,3)/-0.1353193782985104D+00/ DATA NULLW(14,4)/-0.1170027382391999D+00/ DATA NULLW(14,5)/0.1614011431557236D-01/ DATA NULLW(14,6)/0.1330641979525424D+00/ DATA NULLW(14,7)/0.1205891076794731D+00/ DATA NULLW(14,8)/-0.8640919108146020D-02/ DATA NULLW(14,9)/-0.1294253464460428D+00/ DATA NULLW(14,10)/-0.1251272318395094D+00/ DATA NULLW(14,11)/0.0000000000000000D+00/ DATA NULLW(15,1)/0.1384328909795413D+00/ DATA NULLW(15,2)/0.2088816813445404D-01/ DATA NULLW(15,3)/-0.1025551817987029D+00/ DATA NULLW(15,4)/-0.1456993480020755D+00/ DATA NULLW(15,5)/-0.8041833842953963D-01/ DATA NULLW(15,6)/0.4369891359834745D-01/ DATA NULLW(15,7)/0.1355713757017371D+00/ DATA NULLW(15,8)/0.1284341564046552D+00/ DATA NULLW(15,9)/0.2777799655524739D-01/ DATA NULLW(15,10)/-0.9304002613430708D-01/ DATA NULLW(15,11)/-0.1461812140165950D+00/ DATA NULLW(16,1)/0.1434250647895144D+00/ DATA NULLW(16,2)/0.5648832390790171D-01/ DATA NULLW(16,3)/-0.5370731005946019D-01/ DATA NULLW(16,4)/-0.1320553898020212D+00/ DATA NULLW(16,5)/-0.1399117916675364D+00/ DATA NULLW(16,6)/-0.7464504070837483D-01/ DATA NULLW(16,7)/0.2922259459390760D-01/ DATA NULLW(16,8)/0.1177993871727999D+00/ DATA NULLW(16,9)/0.1454239155303997D+00/ DATA NULLW(16,10)/0.9797588771824906D-01/ DATA NULLW(16,11)/0.0000000000000000D+00/ C C***FIRST EXECUTABLE STATEMENT SGAINT C C Define constants C CONST=SAFETY*RCRIT**(1-ALPHA) C C C Compute half-length and center of interval. C HLEN=(B-A)/2 CENTR=A+HLEN X(11)=CENTR C C Compute all abscissas C DO 10 I=1,10 ABSCIS=HLEN*G(I) X(I)=CENTR-ABSCIS X(22-I)=CENTR+ABSCIS 10 CONTINUE IF ((X(21).EQ.X(20)).OR.(X(1).EQ.X(2))) THEN FLAG=1 END IF C C Evaluate all NUMFUN functions in the 21 evaluation points. C CALL FUNSUB (X,NUMFUN,21,FUNVAL) C C Apply the basic rule: first center point C DO 20 J=1,NUMFUN BASVAL(J)=W(11)*FUNVAL(11,J) BASABS(J)=ABS(BASVAL(J)) 20 CONTINUE C C Apply all nullrules: first center point C DO 40 J=1,NUMFUN DO 30 I=1,NUMNUL NULL(I,J)=NULLW(I,11)*FUNVAL(11,J) 30 CONTINUE 40 CONTINUE C C Compute the basic rule contributions from the other points. C DO 60 J=1,NUMFUN DO 50 I=1,10 SUM=FUNVAL(I,J)+FUNVAL(22-I,J) ABSSUM=ABS(FUNVAL(I,J))+ABS(FUNVAL(22-I,J)) BASVAL(J)=W(I)*SUM+BASVAL(J) BASABS(J)=BASABS(J)+W(I)*ABSSUM 50 CONTINUE 60 CONTINUE C C Compute the null rule contributions from the other points. C Recall that one half of the nullrules is symmetric and the other C half is asymmetric. C DO 90 J=1,NUMFUN DO 80 K=1,10 SUM=FUNVAL(K,J)+FUNVAL(22-K,J) DIFF=FUNVAL(K,J)-FUNVAL(22-K,J) DO 70 I=1,NUMNUL/2 I2=2*I I1=I2-1 NULL(I1,J)=NULL(I1,J)+NULLW(I1,K)*SUM NULL(I2,J)=NULL(I2,J)+NULLW(I2,K)*DIFF 70 CONTINUE 80 CONTINUE 90 CONTINUE C C Scale the results with the length of the interval C DO 100 J=1,NUMFUN BASVAL(J)=HLEN*BASVAL(J) BASABS(J)=HLEN*BASABS(J) 100 CONTINUE C C Compute errors. C GREATE=0 DO 130 J=1,NUMFUN EMAX=0 DO 110 I=1,NUMNUL/2 I2=2*I I1=I2-1 E(I)=HLEN*SQRT(NULL(I1,J)**2+NULL(I2,J)**2) EMAX=MAX(EMAX,E(I)) 110 CONTINUE R=0 DO 120 I=1,NUMNUL/2-1 IF (E(I+1).NE.0) THEN RR(I)=E(I)/E(I+1) ELSE RR(I)=1 END IF R=MAX(R,RR(I)) 120 CONTINUE IF (R.GE.1) THEN RGNERR(J)=SAFETY*EMAX ELSE IF (R.GE.RCRIT) THEN RGNERR(J)=SAFETY*R*E(1) ELSE RGNERR(J)=CONST*(R**ALPHA)*E(1) END IF C C Check the noise level. C NOISE=50*EPMACH*BASABS(J) IF ((E(1).LT.NOISE).AND.(E(2).LT.NOISE)) THEN RGNERR(J)=NOISE ELSE RGNERR(J)=MAX(NOISE,RGNERR(J)) END IF C C Check if this is the greatest error so far among the remaining C problems. One exception: If the user wants to force the code C to do extra work (which seems unnecessary). C IF (.NOT.(MORE.AND.DONE(J)).AND.(RGNERR(J).GT.GREATE)) THEN GREATE=RGNERR(J) WORST=J END IF 130 CONTINUE C C Compute the new subdivision points. C C=A+(B-A)/3 D=B-(B-A)/3 RETURN C C***END SGAINT C END CUT HERE............ cat > sgainf.f <<'CUT HERE............' SUBROUTINE SGAINF (NUMFUN,F,B,PERIOD,GAMMA,NEVAL,IFAIL) C C***BEGIN PROLOGUE SGAINF C***KEYWORDS Linear extrapolation, C oscillating functions, asymptotic behavior. C***PURPOSE To compute an estimate of the exponent GAMMA which governs C the asymptotic behavior of an oscillating function: C F(X+PERIOD*K) appr. C(X)/(X+PERIOD*K)**GAMMA, for all X C and all integers K sufficiently great. C The function is assumed to C have asymptotical oscillating behavior and C F(X+PERIOD*(K+1/2)) approx. - F(X+PERIOD*K). C C***LAST MODIFICATION 93-01-31 C C***REFER TO SQAINF C C***DESCRIPTION Computation of the exponent GAMMA > 0 which governs C the asymptotic behavior of the functions F. All functions in the C vector F are assumed to be asymptotically periodic with the same C PERIOD and the same asymptotical behavior. Furthermore we assume C that for all functions in the vector we have C F(X+PERIOD*K) approx. C(X)/(X+PERIOD*K)**GAMMA, for all values of X C and K, such that (X+PERIOD*K) >> B. C may be zero for some values of C X. Assuming that C is non-zero then we have that C F(X+2*PERIOD*K)/F(X+PERIOD*K) approx. (1 + D)**GAMMA, with C D = PERIOD*K/(X+PERIOD*K). C C ON ENTRY C C NUMFUN Integer. C The number of components of the integral. C F Externally declared subroutine for computing C all components of the integrand at all NP C evaluation points. C It must have parameters (X,NUMFUN,NP,FUNVLS) C Input parameters: C X The NP x-coordinates of the evaluation points. C NUMFUN Integer that defines the number of C components of I. C NP Integer that gives the number of evaluation points C Output parameter: C FUNVLS Real array of dimension (NUMFUN,NP) C that defines NUMFUN components of the integrand fo C each evaluation point. C C B Real. C Asymptotic behavior assumed satisfied for all X >= B and C for all the integrands. C PERIOD Real. C All functions are assumed to be asymptotically periodic, C and with the same period. C C ON RETURN C C GAMMA Real C Estimated value of the rate of decay of the functions C when X tends to +infinity. C NEVAL Integer C The number of function evaluations. C IFAIL Integer. C IFAIL = 0 for normal exit. C IFAIL = 12 if either PERIOD is wrong or B(1) is too small. C IFAIL = 13 if unable to estimate GAMMA (In case KEY=2 only) C C***ROUTINES CALLED F,R1MACH C***END PROLOGUE SGAINF C C Global variables. C EXTERNAL F,R1MACH REAL PERIOD,GAMMA,B,R1MACH INTEGER NUMFUN,IFAIL,NEVAL C C Local variables. C REAL P,X,H,XB,FX,FXP,XR,FR,EPMACH,FOLD,FNEW,SIGN,XOLD, +XNEW INTEGER I,IMX,M PARAMETER (IMX=15) REAL A(0:IMX),DA,FACTOR(0:IMX),FI(0:IMX),DF,LAST C C***FIRST EXECUTABLE STATEMENT SGAINF C EPMACH=R1MACH(4) C C First: compute a starting point X >= B where the function value F(X) C is reasonably large in absolute value. Remark: We simply want to C avoid that F(X) becomes approximately 0. C XOLD=B CALL F (XOLD,1,1,FOLD) P=PERIOD/4 XNEW=XOLD DO 10 I=1,3,1 XNEW=XNEW+P CALL F (XNEW,1,1,FNEW) IF (ABS(FOLD).LT.ABS(FNEW)) THEN XOLD=XNEW IF ((FOLD*FNEW).GT.0) THEN P=-P END IF FOLD=FNEW ELSE IF ((FOLD*FNEW).GT.0) THEN XNEW=XOLD END IF END IF P=P/2 10 CONTINUE X=XOLD FX=FOLD C C Check that the function values do change sign C and that the absolute value of the function is decreasing C when evaluated at points which differ by PERIOD*integer/2. C CALL F (X+PERIOD/2,1,1,FXP) NEVAL=5 LAST=FXP SIGN=FX FI(0)=X IF (((SIGN*FXP).LT.0).AND.(ABS(FX).GT.ABS(FXP))) THEN FR=-FXP/FX XR=X/(X+PERIOD/2) A(0)=LOG(FR)/LOG(XR) FACTOR(0)=A(0)*(1/ABS(FR*LOG(FR))+1/ABS(XR*LOG(XR)))*EPMACH ELSE IFAIL=12 RETURN END IF C C Compute a sequence of estimates to GAMMA and extrapolate. C C FACTOR is meant to keep track of the increase in the error due to C the extrapolation. We assume that PERIOD is known correctly to full C precision. If X is an exact local extremum of F(X) then it is C not important to know PERIOD with such a high precision. C However, since X is just a point where ABS(F(X)) is not too small C the quality of PERIOD becomes essential when we estimate GAMMA. C H=PERIOD/2 XB=X DO 30 I=1,IMX H=2*H X=XB+H FI(I)=X CALL F (X+PERIOD/2,1,1,FXP) CALL F (X,1,1,FX) NEVAL=NEVAL+2 C C Check the computed function values. C IF (((SIGN*FX).GT.0).AND.((SIGN*FXP).LT.0).AND.(ABS(FX).GT. + ABS(FXP)).AND.(ABS(LAST).GT.ABS(FX))) THEN LAST=FXP FR=-FXP/FX XR=X/(X+PERIOD/2) A(I)=LOG(FR)/LOG(XR) FACTOR(I)=A(0)*(2**I/ABS(FR*LOG(FR))+1/ABS(XR*LOG(XR)))* + EPMACH ELSE C C Either PERIOD is wrong or B is too small. C IFAIL=12 RETURN END IF DO 20 M=I-1,0,-1 DA=(A(M+1)-A(M))/((FI(I)-FI(M))/FI(M)) DF=(FACTOR(M+1)+FACTOR(M))/((FI(I)-FI(M))/FI(M)) A(M)=A(M+1)+DA FACTOR(M)=FACTOR(M+1)+DF 20 CONTINUE IF (ABS(DA).LE.FACTOR(0)) THEN GAMMA=A(0) C C Normal return C RETURN END IF 30 CONTINUE C C We did not succeed in estimating GAMMA to sufficient precision. C GAMMA=A(0) IFAIL=13 RETURN C C***END SGAINF C END CUT HERE............ cat > sexinf.f <<'CUT HERE............' SUBROUTINE SEXINF (NUMFUN,GAMMA,C,KEY,N,T,CT,UPDATE,EMAX,U,E, +RESULT,ABSERR,EXTERR,BETA,MAXEER) C***BEGIN PROLOGUE SEXINF C***KEYWORDS Quasi-linear extrapolation, infinite integration, C oscillating functions, error estimation. C***PURPOSE To compute better estimates to a vector of approximations C to integrals of oscillating functions over an infinite C interval and to provide new and updated error estimates. C C***LAST MODIFICATION 94-03-11 C C***DESCRIPTION The routine uses linear extrapolation to compute better C approximations to each component in a vector of C integrals. All components are assumed to be integrals of C oscillating functions which asymptotically have the same C behavior at infinity. C A series approach is used, assuming C that the terms are given with estimates of the error in C each term. New error estimates are computed too. The C routine have two options: EITHER a new extrapolation term C is provided and we take a new extrapolation step, C OR we update one previously computed term in the series and C therefore have to update the extrapolation tableau. C C ON ENTRY C C NUMFUN Integer. C The number of components of the integral. C GAMMA Real. C Asymptotic decay of the functions; 1/x**(GAMMA), GAMMA>0. C C Real. C Reference point: User specified subdivision point/period. C KEY Integer C Choice of extrapolation method: C KEY = 0 : Euler's method. C KEY = 1 : Modified Euler. C KEY = 2 : Overholt's polynomial order 2-method. C This last method is the only one that needs C the value of GAMMA. C N Integer C The number of U-terms in the series: 0,1,2,...,N. C Makes it possible to perform N extrapolation steps. C T Real array of dimension (NUMFUN,0:EMAX) C Contains the last row in the extrapolation tableau for C each function in the vector. C CT Real array of dimension NUMFUN. C Contains the element to the left of the diagonal element C in the extrapolation tableau. C UPDATE Integer C If < 0 then this is a new extrapolation step. C If >= 1 then this is a step where we have to correct the C existing tableau. The value of UPDATE gives the index to C the U-element that has been modified. C EMAX Integer C The maximum allowed number of extrapolation steps. C U Real array of dimension (NUMFUN,0:N) C C E Real array of dimension (NUMFUN,0:N) C The estimated errors of all U-terms in the series. C ON RETURN C C T Real array of dimension (NUMFUN:0:EMAX) C Contains the diagonal element in the extrapolation tableau C for each function in the vector after the extrapolation. C CT Real array of dimension NUMFUN. C Contains the element to the left of the diagonal element C in the extrapolation tableau. C RESULT Real array of dimension NUMFUN C Contains the new approximations for the NUMFUN components C of the integral. C ABSERR Real array of dimension NUMFUN. C Returns the global errors for all components. C This includes both the pure extrapolation C error and the effect of not knowing the U-terms exactly. C EXTERR Real array of dimension NUMFUN. C These errors are the pure extrapolation errors. C BETA Real Array of dimension (EMAX +1)*(EMAX +1) C A table of coefficients to be used in the extrapolation C and the error estimation. C MAXEER Real. C The maximum extrapolation error. C***REFERENCES C C T.O.Espelid and K.J.Overholt, DQAINF: An Algorithm for Automatic C Integration of Infinite Oscillating Tails, C submitted for publication 1993. C C K.J.Overholt: The P-algorithms for extrapolation, C Department of Informatics, Report No. 36, 1989. C C***ROUTINES CALLED NONE C***END PROLOGUE SEXINF C C Global variables. C INTEGER N,EMAX,UPDATE,NUMFUN,KEY REAL GAMMA,T(NUMFUN),CT(NUMFUN),C,U(NUMFUN,0:N), +E(NUMFUN,0:N),BETA(0:EMAX,0:EMAX),MAXEER,RESULT(NUMFUN), +ABSERR(NUMFUN),EXTERR(NUMFUN) C C Local variables. C INTEGER I,J REAL SAVE1,SAVE2,CONST,BASE1,BASE2,P,Q PARAMETER (CONST=1.E0) C C CONST heuristic constant used in the error estimation. C C***FIRST EXECUTABLE STATEMENT SEXINF C C Initialize the beta-factors. C BETA(0,0)=1 IF (KEY.EQ.2) THEN BASE1=GAMMA BASE2=MAX(GAMMA-1,4*C) P=2.E0 ELSE BASE1=0.E0 BASE2=MAX(0.E0,4*C) P=1.E0 END IF C C Compute the new extrapolation coefficients. C IF (UPDATE.LT.0) THEN BETA(0,N)=1 BETA(N,0)=1 DO 20 I=1,N-1 SAVE1=1 DO 10 J=N-I,N-1 IF (KEY.EQ.0) THEN Q=0.5E0 ELSE Q=(1+(-BASE1-P*J)/(2*N+BASE2))/2.E0 END IF SAVE2=SAVE1-(SAVE1-BETA(I,J))*Q BETA(I,J)=SAVE1 10 SAVE1=SAVE2 BETA(I,N)=SAVE1 20 CONTINUE DO 30 J=1,N IF (KEY.EQ.0) THEN Q=0.5E0 ELSE Q=(1+(-BASE1-P*(J-1))/(2*N+BASE2))/2.E0 END IF 30 BETA(N,J)=(1-Q)*BETA(N,J-1) END IF C C Extrapolate; use the extrapolation coefficients. C DO 40 J=1,NUMFUN T(J)=0 CT(J)=0 DO 40 I=N,0,-1 CT(J)=CT(J)+BETA(I,N-1)*U(J,I) 40 T(J)=T(J)+BETA(I,N)*U(J,I) C C Then compute the error estimates. C First the extrapolation error and then the U-effect. C The maximum extrapolation error is computed too. C MAXEER=0 DO 60 J=1,NUMFUN EXTERR(J)=CONST*ABS(T(J)-CT(J))/Q MAXEER=MAX(EXTERR(J),MAXEER) C C Note: The last U-errors are effected by the extrapolation-process C Accumulate the total error in exterr. C DO 50 I=0,N 50 EXTERR(J)=EXTERR(J)+BETA(I,N)*E(J,I) 60 CONTINUE C C Define the results: update the results only if the error is improved C DO 70 J=1,NUMFUN IF ((EXTERR(J).LE.ABSERR(J)).OR.(ABSERR(J).EQ.0)) THEN RESULT(J)=T(J) ABSERR(J)=EXTERR(J) END IF 70 CONTINUE RETURN C C***END SEXINF C END CUT HERE............ cat > r1mach.f <<'CUT HERE............' REAL FUNCTION R1MACH(I) C C SINGLE-PRECISION MACHINE CONSTANTS C C R1MACH(1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. C C R1MACH(2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. C C R1MACH(3) = B**(-T), THE SMALLEST RELATIVE SPACING. C C R1MACH(4) = B**(1-T), THE LARGEST RELATIVE SPACING. C C R1MACH(5) = LOG10(B) C C TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT, C THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY C REMOVING THE C FROM COLUMN 1. C C FOR IEEE-ARITHMETIC MACHINES (BINARY STANDARD), THE FIRST C SET OF CONSTANTS BELOW SHOULD BE APPROPRIATE. C C WHERE POSSIBLE, DECIMAL, OCTAL OR HEXADECIMAL CONSTANTS ARE USED C TO SPECIFY THE CONSTANTS EXACTLY. SOMETIMES THIS REQUIRES USING C EQUIVALENT INTEGER ARRAYS. IF YOUR COMPILER USES HALF-WORD C INTEGERS BY DEFAULT (SOMETIMES CALLED INTEGER*2), YOU MAY NEED TO C CHANGE INTEGER TO INTEGER*4 OR OTHERWISE INSTRUCT YOUR COMPILER C TO USE FULL-WORD INTEGERS IN THE NEXT 5 DECLARATIONS. C C COMMENTS JUST BEFORE THE END STATEMENT (LINES STARTING WITH *) C GIVE C SOURCE FOR R1MACH. C INTEGER SMALL(2) INTEGER LARGE(2) INTEGER RIGHT(2) INTEGER DIVER(2) INTEGER LOG10(2) INTEGER SC C REAL RMACH(5) C EQUIVALENCE (RMACH(1),SMALL(1)) EQUIVALENCE (RMACH(2),LARGE(1)) EQUIVALENCE (RMACH(3),RIGHT(1)) EQUIVALENCE (RMACH(4),DIVER(1)) EQUIVALENCE (RMACH(5),LOG10(1)) C C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T C 3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T C PC 7300), AND 8087 BASED MICROS (E.G. IBM PC AND AT&T 6300). DATA SMALL(1) / 8388608 / DATA LARGE(1) / 2139095039 / DATA RIGHT(1) / 864026624 / DATA DIVER(1) / 872415232 / DATA LOG10(1) / 1050288283 /, SC/987/ C MACHINE CONSTANTS FOR AMDAHL MACHINES. C C DATA SMALL(1) / 1048576 / C DATA LARGE(1) / 2147483647 / C DATA RIGHT(1) / 990904320 / C DATA DIVER(1) / 1007681536 / C DATA LOG10(1) / 1091781651 /, SC/987/ C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA RMACH(1) / Z400800000 / C DATA RMACH(2) / Z5FFFFFFFF / C DATA RMACH(3) / Z4E9800000 / C DATA RMACH(4) / Z4EA800000 / C DATA RMACH(5) / Z500E730E8 /, SC/987/ C C MACHINE CONSTANTS FOR THE BURROUGHS 5700/6700/7700 SYSTEMS. C C DATA RMACH(1) / O1771000000000000 / C DATA RMACH(2) / O0777777777777777 / C DATA RMACH(3) / O1311000000000000 / C DATA RMACH(4) / O1301000000000000 / C DATA RMACH(5) / O1157163034761675 /, SC/987/ C C MACHINE CONSTANTS FOR FTN4 ON THE CDC 6000/7000 SERIES. C C DATA RMACH(1) / 00564000000000000000B / C DATA RMACH(2) / 37767777777777777776B / C DATA RMACH(3) / 16414000000000000000B / C DATA RMACH(4) / 16424000000000000000B / C DATA RMACH(5) / 17164642023241175720B /, SC/987/ C C MACHINE CONSTANTS FOR FTN5 ON THE CDC 6000/7000 SERIES. C C DATA RMACH(1) / O"00564000000000000000" / C DATA RMACH(2) / O"37767777777777777776" / C DATA RMACH(3) / O"16414000000000000000" / C DATA RMACH(4) / O"16424000000000000000" / C DATA RMACH(5) / O"17164642023241175720" /, SC/987/ C C MACHINE CONSTANTS FOR CONVEX C-1. C C DATA RMACH(1) / '00800000'X / C DATA RMACH(2) / '7FFFFFFF'X / C DATA RMACH(3) / '34800000'X / C DATA RMACH(4) / '35000000'X / C DATA RMACH(5) / '3F9A209B'X /, SC/987/ C C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. C C DATA RMACH(1) / 200034000000000000000B / C DATA RMACH(2) / 577767777777777777776B / C DATA RMACH(3) / 377224000000000000000B / C DATA RMACH(4) / 377234000000000000000B / C DATA RMACH(5) / 377774642023241175720B /, SC/987/ C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200. C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE - C STATIC RMACH(5) C C DATA SMALL/20K,0/,LARGE/77777K,177777K/ C DATA RIGHT/35420K,0/,DIVER/36020K,0/ C DATA LOG10/40423K,42023K/, SC/987/ C C MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7. C C DATA SMALL(1),SMALL(2) / '20000000, '00000201 / C DATA LARGE(1),LARGE(2) / '37777777, '00000177 / C DATA RIGHT(1),RIGHT(2) / '20000000, '00000352 / C DATA DIVER(1),DIVER(2) / '20000000, '00000353 / C DATA LOG10(1),LOG10(2) / '23210115, '00000377 /, SC/987/ C C MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. C C DATA RMACH(1) / O402400000000 / C DATA RMACH(2) / O376777777777 / C DATA RMACH(3) / O714400000000 / C DATA RMACH(4) / O716400000000 / C DATA RMACH(5) / O776464202324 /, SC/987/ C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86. C C DATA RMACH(1) / Z00100000 / C DATA RMACH(2) / Z7FFFFFFF / C DATA RMACH(3) / Z3B100000 / C DATA RMACH(4) / Z3C100000 / C DATA RMACH(5) / Z41134413 /, SC/987/ C C MACHINE CONSTANTS FOR THE INTERDATA 8/32 C WITH THE UNIX SYSTEM FORTRAN 77 COMPILER. C C FOR THE INTERDATA FORTRAN VII COMPILER REPLACE C THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S. C C DATA RMACH(1) / Z'00100000' / C DATA RMACH(2) / Z'7EFFFFFF' / C DATA RMACH(3) / Z'3B100000' / C DATA RMACH(4) / Z'3C100000' / C DATA RMACH(5) / Z'41134413' /, SC/987/ C C MACHINE CONSTANTS FOR THE PDP-10 (KA OR KI PROCESSOR). C C DATA RMACH(1) / "000400000000 / C DATA RMACH(2) / "377777777777 / C DATA RMACH(3) / "146400000000 / C DATA RMACH(4) / "147400000000 / C DATA RMACH(5) / "177464202324 /, SC/987/ C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1) / 8388608 / C DATA LARGE(1) / 2147483647 / C DATA RIGHT(1) / 880803840 / C DATA DIVER(1) / 889192448 / C DATA LOG10(1) / 1067065499 /, SC/987/ C C DATA RMACH(1) / O00040000000 / C DATA RMACH(2) / O17777777777 / C DATA RMACH(3) / O06440000000 / C DATA RMACH(4) / O06500000000 / C DATA RMACH(5) / O07746420233 /, SC/987/ C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1),SMALL(2) / 128, 0 / C DATA LARGE(1),LARGE(2) / 32767, -1 / C DATA RIGHT(1),RIGHT(2) / 13440, 0 / C DATA DIVER(1),DIVER(2) / 13568, 0 / C DATA LOG10(1),LOG10(2) / 16282, 8347 /, SC/987/ C C DATA SMALL(1),SMALL(2) / O000200, O000000 / C DATA LARGE(1),LARGE(2) / O077777, O177777 / C DATA RIGHT(1),RIGHT(2) / O032200, O000000 / C DATA DIVER(1),DIVER(2) / O032400, O000000 / C DATA LOG10(1),LOG10(2) / O037632, O020233 /, SC/987/ C C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000. C C DATA SMALL(1) / $00800000 / C DATA LARGE(1) / $7F7FFFFF / C DATA RIGHT(1) / $33800000 / C DATA DIVER(1) / $34000000 / C DATA LOG10(1) / $3E9A209B /, SC/987/ C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. C C DATA RMACH(1) / O000400000000 / C DATA RMACH(2) / O377777777777 / C DATA RMACH(3) / O146400000000 / C DATA RMACH(4) / O147400000000 / C DATA RMACH(5) / O177464202324 /, SC/987/ C C MACHINE CONSTANTS FOR THE VAX UNIX F77 COMPILER. C C DATA SMALL(1) / 128 / C DATA LARGE(1) / -32769 / C DATA RIGHT(1) / 13440 / C DATA DIVER(1) / 13568 / C DATA LOG10(1) / 547045274 /, SC/987/ C C MACHINE CONSTANTS FOR THE VAX-11 WITH C FORTRAN IV-PLUS COMPILER. C C DATA RMACH(1) / Z00000080 / C DATA RMACH(2) / ZFFFF7FFF / C DATA RMACH(3) / Z00003480 / C DATA RMACH(4) / Z00003500 / C DATA RMACH(5) / Z209B3F9A /, SC/987/ C C MACHINE CONSTANTS FOR VAX/VMS VERSION 2.2. C C DATA RMACH(1) / '80'X / C DATA RMACH(2) / 'FFFF7FFF'X / C DATA RMACH(3) / '3480'X / C DATA RMACH(4) / '3500'X / C DATA RMACH(5) / '209B3F9A'X /, SC/987/ C C *** ISSUE STOP 778 IF ALL DATA STATEMENTS ARE COMMENTED... IF (SC .NE. 987) STOP 778 IF (I .LT. 1 .OR. I .GT. 5) GOTO 999 R1MACH = RMACH(I) RETURN 999 WRITE(*,1999) I 1999 FORMAT(' R1MACH - I OUT OF BOUNDS',I10) STOP END CUT HERE............ #done