#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