 |
EFTCAMB
Reference documentation for version 3.0
|
This module contains the relevant code for the double precision QUADPACK integration library. More...
Functions/Subroutines | |
| subroutine | xerror (xmess, nmess, nerr, level) |
| XERROR replaces the SLATEC XERROR routine. More... | |
| real(kind=8) function | dqwgts (x, a, b, alfa, beta, integr) |
| DQWGTS defines the weight functions used by DQC25S. More... | |
| real(kind=8) function | dqwgtc (x, c, p2, p3, p4, kp) |
| DQWGTC defines the weight function used by DQC25C. More... | |
| real(kind=8) function | dqwgtf (x, omega, p2, p3, p4, integr) |
| DQWGTF defines the weight functions used by DQC25F. More... | |
| subroutine | dgtsl (n, c, d, e, b, info) |
| DGTSL solves a general tridiagonal linear system. More... | |
| subroutine | dqage (f, a, b, epsabs, epsrel, key, limit, result, abserr, neval, ier, alist, blist, rlist, elist, iord, last) |
| DQAGE estimates a definite integral. More... | |
| subroutine | dqag (f, a, b, epsabs, epsrel, key, result, abserr, neval, ier, limit, lenw, last, iwork, work) |
| DQAG approximates an integral over a finite interval. More... | |
| subroutine | dqagie (f, bound, inf, epsabs, epsrel, limit, result, abserr, neval, ier, alist, blist, rlist, elist, iord, last) |
| DQAGIE estimates an integral over a semi-infinite or infinite interval. More... | |
| subroutine | dqagi (f, bound, inf, epsabs, epsrel, result, abserr, neval, ier, limit, lenw, last, iwork, work) |
| DQAGI estimates an integral over a semi-infinite or infinite interval. More... | |
| subroutine | dqagpe (f, a, b, npts2, points, epsabs, epsrel, limit, result, abserr, neval, ier, alist, blist, rlist, elist, pts, iord, level, ndin, last) |
| DQAGPE computes a definite integral. More... | |
| subroutine | dqagp (f, a, b, npts2, points, epsabs, epsrel, result, abserr, neval, ier, leniw, lenw, last, iwork, work) |
| DQAGP computes a definite integral. More... | |
| subroutine | dqagse (f, a, b, epsabs, epsrel, limit, result, abserr, neval, ier, alist, blist, rlist, elist, iord, last) |
| DQAGSE estimates the integral of a function. More... | |
| subroutine | dqags (f, a, b, epsabs, epsrel, result, abserr, neval, ier, limit, lenw, last, iwork, work) |
| DQAGS estimates the integral of a function. More... | |
| subroutine | dqawce (f, a, b, c, epsabs, epsrel, limit, result, abserr, neval, ier, alist, blist, rlist, elist, iord, last) |
| DQAWCE computes a Cauchy principal value. More... | |
| subroutine | dqawc (f, a, b, c, epsabs, epsrel, result, abserr, neval, ier, limit, lenw, last, iwork, work) |
| DQAWC computes a Cauchy principal value. More... | |
| subroutine | dqawfe (f, a, omega, integr, epsabs, limlst, limit, maxp1, result, abserr, neval, ier, rslst, erlst, ierlst, lst, alist, blist, rlist, elist, iord, nnlog, chebmo) |
| DQAWFE computes Fourier integrals. More... | |
| subroutine | dqawf (f, a, omega, integr, epsabs, result, abserr, neval, ier, limlst, lst, leniw, maxp1, lenw, iwork, work) |
| DQAWF computes Fourier integrals over the interval [ A, +Infinity ). More... | |
| subroutine | dqawoe (f, a, b, omega, integr, epsabs, epsrel, limit, icall, maxp1, result, abserr, neval, ier, last, alist, blist, rlist, elist, iord, nnlog, momcom, chebmo) |
| DQAWOE computes the integrals of oscillatory integrands. More... | |
| subroutine | dqawo (f, a, b, omega, integr, epsabs, epsrel, result, abserr, neval, ier, leniw, maxp1, lenw, last, iwork, work) |
| DQAWO computes the integrals of oscillatory integrands. More... | |
| subroutine | dqawse (f, a, b, alfa, beta, integr, epsabs, epsrel, limit, result, abserr, neval, ier, alist, blist, rlist, elist, iord, last) |
| DQAWSE estimates integrals with algebraico-logarithmic end singularities. More... | |
| subroutine | dqaws (f, a, b, alfa, beta, integr, epsabs, epsrel, result, abserr, neval, ier, limit, lenw, last, iwork, work) |
| DQAWS estimates integrals with algebraico-logarithmic endpoint singularities. More... | |
| subroutine | dqc25c (f, a, b, c, result, abserr, krul, neval) |
| DQC25C returns integration rules for Cauchy Principal Value integrals. More... | |
| subroutine | dqc25f (f, a, b, omega, integr, nrmom, maxp1, ksave, result, abserr, neval, resabs, resasc, momcom, chebmo) |
| DQC25F returns integration rules for functions with a COS or SIN factor. More... | |
| subroutine | dqc25s (f, a, b, bl, br, alfa, beta, ri, rj, rg, rh, result, abserr, resasc, integr, nev) |
| DQC25S returns rules for algebraico-logarithmic end point singularities. More... | |
| subroutine | dqcheb (x, fval, cheb12, cheb24) |
| DQCHEB computes the Chebyshev series expansion. More... | |
| subroutine | dqelg (n, epstab, result, abserr, res3la, nres) |
| DQELG carries out the Epsilon extrapolation algorithm. More... | |
| subroutine | dqk15 (f, a, b, result, abserr, resabs, resasc) |
| DQK15 carries out a 15 point Gauss-Kronrod quadrature rule. More... | |
| subroutine | dqk15i (f, boun, inf, a, b, result, abserr, resabs, resasc) |
| DQK15I applies a 15 point Gauss-Kronrod quadrature on an infinite interval. More... | |
| subroutine | dqk15w (f, w, p1, p2, p3, p4, kp, a, b, result, abserr, resabs, resasc) |
| DQK15W applies a 15 point Gauss-Kronrod rule for a weighted integrand. More... | |
| subroutine | dqk21 (f, a, b, result, abserr, resabs, resasc) |
| DQK21 carries out a 21 point Gauss-Kronrod quadrature rule. More... | |
| subroutine | dqk31 (f, a, b, result, abserr, resabs, resasc) |
| DQK31 carries out a 31 point Gauss-Kronrod quadrature rule. More... | |
| subroutine | dqk41 (f, a, b, result, abserr, resabs, resasc) |
| DQK41 carries out a 41 point Gauss-Kronrod quadrature rule. More... | |
| subroutine | dqk51 (f, a, b, result, abserr, resabs, resasc) |
| DQK51 carries out a 51 point Gauss-Kronrod quadrature rule. More... | |
| subroutine | dqk61 (f, a, b, result, abserr, resabs, resasc) |
| DQK61 carries out a 61 point Gauss-Kronrod quadrature rule. More... | |
| subroutine | dqmomo (alfa, beta, ri, rj, rg, rh, integr) |
| DQMOMO computes modified Chebyshev moments. More... | |
| subroutine | dqng (f, a, b, epsabs, epsrel, result, abserr, neval, ier) |
| DQNG estimates an integral, using non-adaptive integration. More... | |
| subroutine | dqpsrt (limit, last, maxerr, ermax, elist, iord, nrmax) |
| DQPSRT maintains the order of a list of local error estimates. More... | |
This module contains the relevant code for the double precision QUADPACK integration library.
| subroutine eftcamb_quadpack::dgtsl | ( | integer ( kind = 4 ) | n, |
| real ( kind = 8 ), dimension(n) | c, | ||
| real ( kind = 8 ), dimension(n) | d, | ||
| real ( kind = 8 ), dimension(n) | e, | ||
| real ( kind = 8 ), dimension(n) | b, | ||
| integer ( kind = 4 ) | info | ||
| ) |
DGTSL solves a general tridiagonal linear system.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
17 May 2005
Author:
FORTRAN90 version by John Burkardt.
Reference:
Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, LINPACK User's Guide, SIAM, 1979, ISBN13: 978-0-898711-72-1, LC: QA214.L56.
Parameters:
Input, integer ( kind = 4 ) N, the order of the tridiagonal matrix.
Input/output, real ( kind = 8 ) C(N), contains the subdiagonal of the tridiagonal matrix in entries C(2:N). On output, C is destroyed.
Input/output, real ( kind = 8 ) D(N). On input, the diagonal of the matrix. On output, D is destroyed.
Input/output, real ( kind = 8 ) E(N), contains the superdiagonal of the tridiagonal matrix in entries E(1:N-1). On output E is destroyed.
Input/output, real ( kind = 8 ) B(N). On input, the right hand side. On output, the solution.
Output, integer ( kind = 4 ) INFO, error flag. 0, normal value. K, the K-th element of the diagonal becomes exactly zero. The routine returns if this error condition is detected.
Definition at line 193 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqag | ( | real ( kind = 8 ), external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | epsabs, | ||
| real ( kind = 8 ) | epsrel, | ||
| integer ( kind = 4 ) | key, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| integer ( kind = 4 ) | neval, | ||
| integer ( kind = 4 ) | ier, | ||
| integer ( kind = 4 ) | limit, | ||
| integer ( kind = 4 ) | lenw, | ||
| integer ( kind = 4 ) | last, | ||
| integer ( kind = 4 ), dimension(limit) | iwork, | ||
| real ( kind = 8 ), dimension(lenw) | work | ||
| ) |
DQAG approximates an integral over a finite interval.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose the routine calculates an approximation result to a given definite integral i = integral of f over (a,b), hopefully satisfying following claim for accuracy abs(i-result)le.max(epsabs,epsrel*abs(i)).
Parameters:
f - real ( kind = 8 )
function subprogam defining the integrand
function f(x). the actual name for f needs to be
declared e x t e r n a l in the driver program.
a - real ( kind = 8 )
lower limit of integration
b - real ( kind = 8 )
upper limit of integration
epsabs - real ( kind = 8 )
absolute accoracy requested
epsrel - real ( kind = 8 )
relative accuracy requested
if epsabs.le.0
and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
the routine will end with ier = 6.
key - integer ( kind = 4 )
key for choice of local integration rule
a gauss-kronrod pair is used with
7 - 15 points if key.lt.2,
10 - 21 points if key = 2,
15 - 31 points if key = 3,
20 - 41 points if key = 4,
25 - 51 points if key = 5,
30 - 61 points if key.gt.5.
on return result - real ( kind = 8 ) approximation to the integral
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
neval - integer ( kind = 4 ) number of integrand evaluations
ier - integer ( kind = 4 ) ier = 0 normal and reliable termination of the routine. it is assumed that the requested accuracy has been achieved. ier.gt.0 abnormal termination of the routine the estimates for result and error are less reliable. it is assumed that the requested accuracy has not been achieved. error messages ier = 1 maximum number of subdivisions allowed has been achieved. one can allow more subdivisions by increasing the value of limit (and taking the according dimension adjustments into account). however, if this yield no improvement it is advised to analyze the integrand in order to determine the integration difficulaties. if the position of a local difficulty can be determined (i.e.singularity, discontinuity within the interval) one will probably gain from splitting up the interval at this point and calling the integrator on the subranges. if possible, an appropriate special-purpose integrator should be used which is designed for handling the type of difficulty involved. = 2 the occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. = 3 extremely bad integrand behaviour occurs at some points of the integration interval. = 6 the input is invalid, because (epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) or limit.lt.1 or lenw.lt.limit*4. result, abserr, neval, last are set to zero. except when lenw is invalid, iwork(1), work(limit*2+1) and work(limit*3+1) are set to zero, work(1) is set to a and work(limit+1) to b.
dimensioning parameters limit - integer ( kind = 4 ) dimensioning parameter for iwork limit determines the maximum number of subintervals in the partition of the given integration interval (a,b), limit.ge.1. if limit.lt.1, the routine will end with ier = 6.
lenw - integer ( kind = 4 ) dimensioning parameter for work lenw must be at least limit*4. if lenw.lt.limit*4, the routine will end with ier = 6.
last - integer ( kind = 4 ) on return, last equals the number of subintervals produced in the subdivision process, which determines the number of significant elements actually in the work arrays.
work arrays iwork - integer ( kind = 4 ) vector of dimension at least limit, the first k elements of which contain pointers to the error estimates over the subintervals, such that work(limit*3+iwork(1)),... , work(limit*3+iwork(k)) form a decreasing sequence with k = last if last.le.(limit/2+2), and k = limit+1-last otherwise
work - real ( kind = 8 ) vector of dimension at least lenw on return work(1), ..., work(last) contain the left end points of the subintervals in the partition of (a,b), work(limit+1), ..., work(limit+last) contain the right end points, work(limit*2+1), ..., work(limit*2+last) contain the integral approximations over the subintervals, work(limit*3+1), ..., work(limit*3+last) contain the error estimates.
Definition at line 806 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqage | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | epsabs, | ||
| real ( kind = 8 ) | epsrel, | ||
| integer ( kind = 4 ) | key, | ||
| integer ( kind = 4 ) | limit, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| integer ( kind = 4 ) | neval, | ||
| integer ( kind = 4 ) | ier, | ||
| real ( kind = 8 ), dimension(limit) | alist, | ||
| real ( kind = 8 ), dimension(limit) | blist, | ||
| real ( kind = 8 ), dimension(limit) | rlist, | ||
| real ( kind = 8 ), dimension(limit) | elist, | ||
| integer ( kind = 4 ), dimension(limit) | iord, | ||
| integer ( kind = 4 ) | last | ||
| ) |
DQAGE estimates a definite integral.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose the routine calculates an approximation result to a given definite integral i = integral of f over (a,b), hopefully satisfying following claim for accuracy abs(i-reslt).le.max(epsabs,epsrel*abs(i)).
Parameters:
on entry f - real ( kind = 8 ) function subprogram defining the integrand function f(x). the actual name for f needs to be declared e x t e r n a l in the driver program.
a - real ( kind = 8 ) lower limit of integration
b - real ( kind = 8 ) upper limit of integration
epsabs - real ( kind = 8 ) absolute accuracy requested epsrel - real ( kind = 8 ) relative accuracy requested if epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), the routine will end with ier = 6.
key - integer ( kind = 4 ) key for choice of local integration rule a gauss-kronrod pair is used with 7 - 15 points if key.lt.2, 10 - 21 points if key = 2, 15 - 31 points if key = 3, 20 - 41 points if key = 4, 25 - 51 points if key = 5, 30 - 61 points if key.gt.5.
limit - integer ( kind = 4 ) gives an upperbound on the number of subintervals in the partition of (a,b), limit.ge.1.
on return result - real ( kind = 8 ) approximation to the integral
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
neval - integer ( kind = 4 ) number of integrand evaluations
ier - integer ( kind = 4 ) ier = 0 normal and reliable termination of the routine. it is assumed that the requested accuracy has been achieved. ier.gt.0 abnormal termination of the routine the estimates for result and error are less reliable. it is assumed that the requested accuracy has not been achieved. error messages ier = 1 maximum number of subdivisions allowed has been achieved. one can allow more subdivisions by increasing the value of limit. however, if this yields no improvement it is rather advised to analyze the integrand in order to determine the integration difficulties. if the position of a local difficulty can be determined(e.g. singularity, discontinuity within the interval) one will probably gain from splitting up the interval at this point and calling the integrator on the subranges. if possible, an appropriate special-purpose integrator should be used which is designed for handling the type of difficulty involved. = 2 the occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. = 3 extremely bad integrand behaviour occurs at some points of the integration interval. = 6 the input is invalid, because (epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), result, abserr, neval, last, rlist(1) , elist(1) and iord(1) are set to zero. alist(1) and blist(1) are set to a and b respectively.
alist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the left end points of the subintervals in the partition of the given integration range (a,b)
blist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the right end points of the subintervals in the partition of the given integration range (a,b)
rlist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the integral approximations on the subintervals
elist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the moduli of the absolute error estimates on the subintervals
iord - integer ( kind = 4 ) vector of dimension at least limit, the first k elements of which are pointers to the error estimates over the subintervals, such that elist(iord(1)), ..., elist(iord(k)) form a decreasing sequence, with k = last if last.le.(limit/2+2), and k = limit+1-last otherwise
last - integer ( kind = 4 ) number of subintervals actually produced in the subdivision process
Local Parameters:
alist - list of left end points of all subintervals considered up to now blist - list of right end points of all subintervals considered up to now rlist(i) - approximation to the integral over (alist(i),blist(i)) elist(i) - error estimate applying to rlist(i) maxerr - pointer to the interval with largest error estimate errmax - elist(maxerr) area - sum of the integrals over the subintervals errsum - sum of the errors over the subintervals errbnd - requested accuracy max(epsabs,epsrel* abs(result)) *****1 - variable for the left subinterval *****2 - variable for the right subinterval last - index for subdivision
machine dependent constants epmach is the largest relative spacing. uflow is the smallest positive magnitude.
Definition at line 447 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqagi | ( | external | f, |
| real ( kind = 8 ) | bound, | ||
| integer ( kind = 4 ) | inf, | ||
| real ( kind = 8 ) | epsabs, | ||
| real ( kind = 8 ) | epsrel, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| integer ( kind = 4 ) | neval, | ||
| integer ( kind = 4 ) | ier, | ||
| integer ( kind = 4 ) | limit, | ||
| integer ( kind = 4 ) | lenw, | ||
| integer ( kind = 4 ) | last, | ||
| integer ( kind = 4 ), dimension(limit) | iwork, | ||
| real ( kind = 8 ), dimension(lenw) | work | ||
| ) |
DQAGI estimates an integral over a semi-infinite or infinite interval.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose the routine calculates an approximation result to a given integral i = integral of f over (bound,+infinity) or i = integral of f over (-infinity,bound) or i = integral of f over (-infinity,+infinity) hopefully satisfying following claim for accuracy abs(i-result).le.max(epsabs,epsrel*abs(i)).
Parameters:
on entry f - real ( kind = 8 ) function subprogram defining the integrand function f(x). the actual name for f needs to be declared e x t e r n a l in the driver program.
bound - real ( kind = 8 ) finite bound of integration range (has no meaning if interval is doubly-infinite)
inf - integer ( kind = 4 ) indicating the kind of integration range involved inf = 1 corresponds to (bound,+infinity), inf = -1 to (-infinity,bound), inf = 2 to (-infinity,+infinity).
epsabs - real ( kind = 8 ) absolute accuracy requested epsrel - real ( kind = 8 ) relative accuracy requested if epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), the routine will end with ier = 6.
on return result - real ( kind = 8 ) approximation to the integral
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
neval - integer ( kind = 4 ) number of integrand evaluations
ier - integer ( kind = 4 ) ier = 0 normal and reliable termination of the routine. it is assumed that the requested accuracy has been achieved.
dimensioning parameters limit - integer ( kind = 4 ) dimensioning parameter for iwork limit determines the maximum number of subintervals in the partition of the given integration interval (a,b), limit.ge.1. if limit.lt.1, the routine will end with ier = 6.
lenw - integer ( kind = 4 ) dimensioning parameter for work lenw must be at least limit*4. if lenw.lt.limit*4, the routine will end with ier = 6.
last - integer ( kind = 4 ) on return, last equals the number of subintervals produced in the subdivision process, which determines the number of significant elements actually in the work arrays.
work arrays iwork - integer ( kind = 4 ) vector of dimension at least limit, the first k elements of which contain pointers to the error estimates over the subintervals, such that work(limit*3+iwork(1)),... , work(limit*3+iwork(k)) form a decreasing sequence, with k = last if last.le.(limit/2+2), and k = limit+1-last otherwise
work - real ( kind = 8 ) vector of dimension at least lenw on return work(1), ..., work(last) contain the left end points of the subintervals in the partition of (a,b), work(limit+1), ..., work(limit+last) contain the right end points, work(limit*2+1), ...,work(limit*2+last) contain the integral approximations over the subintervals, work(limit*3+1), ..., work(limit*3) contain the error estimates.
Definition at line 1478 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqagie | ( | external | f, |
| real ( kind = 8 ) | bound, | ||
| integer ( kind = 4 ) | inf, | ||
| real ( kind = 8 ) | epsabs, | ||
| real ( kind = 8 ) | epsrel, | ||
| integer ( kind = 4 ) | limit, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| integer ( kind = 4 ) | neval, | ||
| integer ( kind = 4 ) | ier, | ||
| real ( kind = 8 ), dimension(limit) | alist, | ||
| real ( kind = 8 ), dimension(limit) | blist, | ||
| real ( kind = 8 ), dimension(limit) | rlist, | ||
| real ( kind = 8 ), dimension(limit) | elist, | ||
| integer ( kind = 4 ), dimension(limit) | iord, | ||
| integer ( kind = 4 ) | last | ||
| ) |
DQAGIE estimates an integral over a semi-infinite or infinite interval.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose the routine calculates an approximation result to a given integral i = integral of f over (bound,+infinity) or i = integral of f over (-infinity,bound) or i = integral of f over (-infinity,+infinity), hopefully satisfying following claim for accuracy abs(i-result).le.max(epsabs,epsrel*abs(i))
Parameters:
f - real ( kind = 8 )
function subprogram defining the integrand
function f(x). the actual name for f needs to be
declared e x t e r n a l in the driver program.
bound - real ( kind = 8 )
finite bound of integration range
(has no meaning if interval is doubly-infinite)
inf - real ( kind = 8 )
indicating the kind of integration range involved
inf = 1 corresponds to (bound,+infinity),
inf = -1 to (-infinity,bound),
inf = 2 to (-infinity,+infinity).
epsabs - real ( kind = 8 )
absolute accuracy requested
epsrel - real ( kind = 8 )
relative accuracy requested
if epsabs.le.0
and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
the routine will end with ier = 6.
limit - integer ( kind = 4 )
gives an upper bound on the number of subintervals
in the partition of (a,b), limit.ge.1
on return result - real ( kind = 8 ) approximation to the integral
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
neval - integer ( kind = 4 ) number of integrand evaluations
ier - integer ( kind = 4 ) ier = 0 normal and reliable termination of the routine. it is assumed that the requested accuracy has been achieved.
alist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the left end points of the subintervals in the partition of the transformed integration range (0,1).
blist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the right end points of the subintervals in the partition of the transformed integration range (0,1).
rlist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the integral approximations on the subintervals
elist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the moduli of the absolute error estimates on the subintervals
iord - integer ( kind = 4 ) vector of dimension limit, the first k elements of which are pointers to the error estimates over the subintervals, such that elist(iord(1)), ..., elist(iord(k)) form a decreasing sequence, with k = last if last.le.(limit/2+2), and k = limit+1-last otherwise
last - integer ( kind = 4 ) number of subintervals actually produced in the subdivision process
Local Parameters:
the dimension of rlist2 is determined by the value of limexp in routine dqelg.
alist - list of left end points of all subintervals considered up to now blist - list of right end points of all subintervals considered up to now rlist(i) - approximation to the integral over (alist(i),blist(i)) rlist2 - array of dimension at least (limexp+2), containing the part of the epsilon table wich is still needed for further computations elist(i) - error estimate applying to rlist(i) maxerr - pointer to the interval with largest error estimate errmax - elist(maxerr) erlast - error on the interval currently subdivided (before that subdivision has taken place) area - sum of the integrals over the subintervals errsum - sum of the errors over the subintervals errbnd - requested accuracy max(epsabs,epsrel* abs(result)) *****1 - variable for the left subinterval *****2 - variable for the right subinterval last - index for subdivision nres - number of calls to the extrapolation routine numrl2 - number of elements currently in rlist2. if an appropriate approximation to the compounded integral has been obtained, it is put in rlist2(numrl2) after numrl2 has been increased by one. small - length of the smallest interval considered up to now, multiplied by 1.5 erlarg - sum of the errors over the intervals larger than the smallest interval considered up to now extrap - logical variable denoting that the routine is attempting to perform extrapolation. i.e. before subdividing the smallest interval we try to decrease the value of erlarg. noext - logical variable denoting that extrapolation is no longer allowed (true-value)
machine dependent constants
epmach is the largest relative spacing. uflow is the smallest positive magnitude. oflow is the largest positive magnitude.
Definition at line 1055 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqagp | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| integer ( kind = 4 ) | npts2, | ||
| real ( kind = 8 ), dimension(npts2) | points, | ||
| real ( kind = 8 ) | epsabs, | ||
| real ( kind = 8 ) | epsrel, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| integer ( kind = 4 ) | neval, | ||
| integer ( kind = 4 ) | ier, | ||
| integer ( kind = 4 ) | leniw, | ||
| integer ( kind = 4 ) | lenw, | ||
| integer ( kind = 4 ) | last, | ||
| integer ( kind = 4 ), dimension(leniw) | iwork, | ||
| real ( kind = 8 ), dimension(lenw) | work | ||
| ) |
DQAGP computes a definite integral.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose the routine calculates an approximation result to a given definite integral i = integral of f over (a,b), hopefully satisfying following claim for accuracy break points of the integration interval, where local difficulties of the integrand may occur (e.g. singularities, discontinuities), are provided by the user.
Parameters:
on entry f - real ( kind = 8 ) function subprogram defining the integrand function f(x). the actual name for f needs to be declared e x t e r n a l in the driver program.
a - real ( kind = 8 ) lower limit of integration
b - real ( kind = 8 ) upper limit of integration
npts2 - integer ( kind = 4 ) number equal to two more than the number of user-supplied break points within the integration range, npts.ge.2. if npts2.lt.2, the routine will end with ier = 6.
points - real ( kind = 8 ) vector of dimension npts2, the first (npts2-2) elements of which are the user provided break points. if these points do not constitute an ascending sequence there will be an automatic sorting.
epsabs - real ( kind = 8 ) absolute accuracy requested epsrel - real ( kind = 8 ) relative accuracy requested if epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), the routine will end with ier = 6.
on return result - real ( kind = 8 ) approximation to the integral
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
neval - integer ( kind = 4 ) number of integrand evaluations
ier - integer ( kind = 4 ) ier = 0 normal and reliable termination of the routine. it is assumed that the requested accuracy has been achieved. ier.gt.0 abnormal termination of the routine. the estimates for integral and error are less reliable. it is assumed that the requested accuracy has not been achieved. error messages ier = 1 maximum number of subdivisions allowed has been achieved. one can allow more subdivisions by increasing the value of limit (and taking the according dimension adjustments into account). however, if this yields no improvement it is advised to analyze the integrand in order to determine the integration difficulties. if the position of a local difficulty can be determined (i.e. singularity, discontinuity within the interval), it should be supplied to the routine as an element of the vector points. if necessary an appropriate special-purpose integrator must be used, which is designed for handling the type of difficulty involved. = 2 the occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. the error may be under-estimated. = 3 extremely bad integrand behaviour occurs at some points of the integration interval. = 4 the algorithm does not converge. roundoff error is detected in the extrapolation table. it is presumed that the requested tolerance cannot be achieved, and that the returned result is the best which can be obtained. = 5 the integral is probably divergent, or slowly convergent. it must be noted that divergence can occur with any other value of ier.gt.0. = 6 the input is invalid because npts2.lt.2 or break points are specified outside the integration range or (epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) result, abserr, neval, last are set to zero. exept when leniw or lenw or npts2 is invalid, iwork(1), iwork(limit+1), work(limit*2+1) and work(limit*3+1) are set to zero. work(1) is set to a and work(limit+1) to b (where limit = (leniw-npts2)/2).
dimensioning parameters leniw - integer ( kind = 4 ) dimensioning parameter for iwork leniw determines limit = (leniw-npts2)/2, which is the maximum number of subintervals in the partition of the given integration interval (a,b), leniw.ge.(3*npts2-2). if leniw.lt.(3*npts2-2), the routine will end with ier = 6.
lenw - integer ( kind = 4 ) dimensioning parameter for work lenw must be at least leniw*2-npts2. if lenw.lt.leniw*2-npts2, the routine will end with ier = 6.
last - integer ( kind = 4 ) on return, last equals the number of subintervals produced in the subdivision process, which determines the number of significant elements actually in the work arrays.
work arrays iwork - integer ( kind = 4 ) vector of dimension at least leniw. on return, the first k elements of which contain pointers to the error estimates over the subintervals, such that work(limit*3+iwork(1)),..., work(limit*3+iwork(k)) form a decreasing sequence, with k = last if last.le.(limit/2+2), and k = limit+1-last otherwise iwork(limit+1), ...,iwork(limit+last) contain the subdivision levels of the subintervals, i.e. if (aa,bb) is a subinterval of (p1,p2) where p1 as well as p2 is a user-provided break point or integration limit, then (aa,bb) has level l if abs(bb-aa) = abs(p2-p1)*2**(-l), iwork(limit*2+1), ..., iwork(limit*2+npts2) have no significance for the user, note that limit = (leniw-npts2)/2.
work - real ( kind = 8 ) vector of dimension at least lenw on return work(1), ..., work(last) contain the left end points of the subintervals in the partition of (a,b), work(limit+1), ..., work(limit+last) contain the right end points, work(limit*2+1), ..., work(limit*2+last) contain the integral approximations over the subintervals, work(limit*3+1), ..., work(limit*3+last) contain the corresponding error estimates, work(limit*4+1), ..., work(limit*4+npts2) contain the integration limits and the break points sorted in an ascending sequence. note that limit = (leniw-npts2)/2.
Definition at line 2252 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqagpe | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| integer ( kind = 4 ) | npts2, | ||
| real ( kind = 8 ), dimension(npts2) | points, | ||
| real ( kind = 8 ) | epsabs, | ||
| real ( kind = 8 ) | epsrel, | ||
| integer ( kind = 4 ) | limit, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| integer ( kind = 4 ) | neval, | ||
| integer ( kind = 4 ) | ier, | ||
| real ( kind = 8 ), dimension(limit) | alist, | ||
| real ( kind = 8 ), dimension(limit) | blist, | ||
| real ( kind = 8 ), dimension(limit) | rlist, | ||
| real ( kind = 8 ), dimension(limit) | elist, | ||
| real ( kind = 8 ), dimension(npts2) | pts, | ||
| integer ( kind = 4 ), dimension(limit) | iord, | ||
| integer ( kind = 4 ), dimension(limit) | level, | ||
| integer ( kind = 4 ), dimension(npts2) | ndin, | ||
| integer ( kind = 4 ) | last | ||
| ) |
DQAGPE computes a definite integral.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose the routine calculates an approximation result to a given definite integral i = integral of f over (a,b), hopefully satisfying following claim for accuracy abs(i-result).le. max(epsabs,epsrel*abs(i)). break points of the integration interval, where local difficulties of the integrand may occur(e.g. singularities,discontinuities),provided by user.
Parameters:
on entry f - real ( kind = 8 ) function subprogram defining the integrand function f(x). the actual name for f needs to be declared e x t e r n a l in the driver program.
a - real ( kind = 8 ) lower limit of integration
b - real ( kind = 8 ) upper limit of integration
npts2 - integer ( kind = 4 ) number equal to two more than the number of user-supplied break points within the integration range, npts2.ge.2. if npts2.lt.2, the routine will end with ier = 6.
points - real ( kind = 8 ) vector of dimension npts2, the first (npts2-2) elements of which are the user provided break points. if these points do not constitute an ascending sequence there will be an automatic sorting.
epsabs - real ( kind = 8 ) absolute accuracy requested epsrel - real ( kind = 8 ) relative accuracy requested if epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), the routine will end with ier = 6.
limit - integer ( kind = 4 ) gives an upper bound on the number of subintervals in the partition of (a,b), limit.ge.npts2 if limit.lt.npts2, the routine will end with ier = 6.
on return result - real ( kind = 8 ) approximation to the integral
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
neval - integer ( kind = 4 ) number of integrand evaluations
ier - integer ( kind = 4 ) ier = 0 normal and reliable termination of the routine. it is assumed that the requested accuracy has been achieved. ier.gt.0 abnormal termination of the routine. the estimates for integral and error are less reliable. it is assumed that the requested accuracy has not been achieved. error messages ier = 1 maximum number of subdivisions allowed has been achieved. one can allow more subdivisions by increasing the value of limit (and taking the according dimension adjustments into account). however, if this yields no improvement it is advised to analyze the integrand in order to determine the integration difficulties. if the position of a local difficulty can be determined (i.e. singularity, discontinuity within the interval), it should be supplied to the routine as an element of the vector points. if necessary an appropriate special-purpose integrator must be used, which is designed for handling the type of difficulty involved. = 2 the occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. the error may be under-estimated. = 3 extremely bad integrand behaviour occurs at some points of the integration interval. = 4 the algorithm does not converge. roundoff error is detected in the extrapolation table. it is presumed that the requested tolerance cannot be achieved, and that the returned result is the best which can be obtained. = 5 the integral is probably divergent, or slowly convergent. it must be noted that divergence can occur with any other value of ier.gt.0. = 6 the input is invalid because npts2.lt.2 or break points are specified outside the integration range or (epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) or limit.lt.npts2. result, abserr, neval, last, rlist(1), and elist(1) are set to zero. alist(1) and blist(1) are set to a and b respectively.
alist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the left end points of the subintervals in the partition of the given integration range (a,b)
blist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the right end points of the subintervals in the partition of the given integration range (a,b)
rlist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the integral approximations on the subintervals
elist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the moduli of the absolute error estimates on the subintervals
pts - real ( kind = 8 ) vector of dimension at least npts2, containing the integration limits and the break points of the interval in ascending sequence.
level - integer ( kind = 4 ) vector of dimension at least limit, containing the subdivision levels of the subinterval, i.e. if (aa,bb) is a subinterval of (p1,p2) where p1 as well as p2 is a user-provided break point or integration limit, then (aa,bb) has level l if abs(bb-aa) = abs(p2-p1)*2**(-l).
ndin - integer ( kind = 4 ) vector of dimension at least npts2, after first integration over the intervals (pts(i)),pts(i+1), i = 0,1, ..., npts2-2, the error estimates over some of the intervals may have been increased artificially, in order to put their subdivision forward. if this happens for the subinterval numbered k, ndin(k) is put to 1, otherwise ndin(k) = 0.
iord - integer ( kind = 4 ) vector of dimension at least limit, the first k elements of which are pointers to the error estimates over the subintervals, such that elist(iord(1)), ..., elist(iord(k)) form a decreasing sequence, with k = last if last.le.(limit/2+2), and k = limit+1-last otherwise
last - integer ( kind = 4 ) number of subintervals actually produced in the subdivisions process
Local Parameters:
the dimension of rlist2 is determined by the value of limexp in routine epsalg (rlist2 should be of dimension (limexp+2) at least).
alist - list of left end points of all subintervals considered up to now blist - list of right end points of all subintervals considered up to now rlist(i) - approximation to the integral over (alist(i),blist(i)) rlist2 - array of dimension at least limexp+2 containing the part of the epsilon table which is still needed for further computations elist(i) - error estimate applying to rlist(i) maxerr - pointer to the interval with largest error estimate errmax - elist(maxerr) erlast - error on the interval currently subdivided (before that subdivision has taken place) area - sum of the integrals over the subintervals errsum - sum of the errors over the subintervals errbnd - requested accuracy max(epsabs,epsrel* abs(result)) *****1 - variable for the left subinterval *****2 - variable for the right subinterval last - index for subdivision nres - number of calls to the extrapolation routine numrl2 - number of elements in rlist2. if an appropriate approximation to the compounded integral has been obtained, it is put in rlist2(numrl2) after numrl2 has been increased by one. erlarg - sum of the errors over the intervals larger than the smallest interval considered up to now extrap - logical variable denoting that the routine is attempting to perform extrapolation. i.e. before subdividing the smallest interval we try to decrease the value of erlarg. noext - logical variable denoting that extrapolation is no longer allowed (true-value)
machine dependent constants
epmach is the largest relative spacing. uflow is the smallest positive magnitude. oflow is the largest positive magnitude.
Definition at line 1750 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqags | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | epsabs, | ||
| real ( kind = 8 ) | epsrel, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| integer ( kind = 4 ) | neval, | ||
| integer ( kind = 4 ) | ier, | ||
| integer ( kind = 4 ) | limit, | ||
| integer ( kind = 4 ) | lenw, | ||
| integer ( kind = 4 ) | last, | ||
| integer ( kind = 4 ), dimension(limit) | iwork, | ||
| real ( kind = 8 ), dimension(lenw) | work | ||
| ) |
DQAGS estimates the integral of a function.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose the routine calculates an approximation result to a given definite integral i = integral of f over (a,b), hopefully satisfying following claim for accuracy abs(i-result).le.max(epsabs,epsrel*abs(i)).
Parameters:
on entry f - real ( kind = 8 ) function subprogram defining the integrand function f(x). the actual name for f needs to be declared e x t e r n a l in the driver program.
a - real ( kind = 8 ) lower limit of integration
b - real ( kind = 8 ) upper limit of integration
epsabs - real ( kind = 8 ) absolute accuracy requested epsrel - real ( kind = 8 ) relative accuracy requested if epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), the routine will end with ier = 6.
on return result - real ( kind = 8 ) approximation to the integral
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
neval - integer ( kind = 4 ) number of integrand evaluations
ier - integer ( kind = 4 ) ier = 0 normal and reliable termination of the routine. it is assumed that the requested accuracy has been achieved. ier.gt.0 abnormal termination of the routine the estimates for integral and error are less reliable. it is assumed that the requested accuracy has not been achieved. error messages ier = 1 maximum number of subdivisions allowed has been achieved. one can allow more sub- divisions by increasing the value of limit (and taking the according dimension adjustments into account. however, if this yields no improvement it is advised to analyze the integrand in order to determine the integration difficulties. if the position of a local difficulty can be determined (e.g. singularity, discontinuity within the interval) one will probably gain from splitting up the interval at this point and calling the integrator on the subranges. if possible, an appropriate special-purpose integrator should be used, which is designed for handling the type of difficulty involved. = 2 the occurrence of roundoff error is detec- ted, which prevents the requested tolerance from being achieved. the error may be under-estimated. = 3 extremely bad integrand behaviour occurs at some points of the integration interval. = 4 the algorithm does not converge. roundoff error is detected in the extrapolation table. it is presumed that the requested tolerance cannot be achieved, and that the returned result is the best which can be obtained. = 5 the integral is probably divergent, or slowly convergent. it must be noted that divergence can occur with any other value of ier. = 6 the input is invalid, because (epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28) or limit.lt.1 or lenw.lt.limit*4. result, abserr, neval, last are set to zero.except when limit or lenw is invalid, iwork(1), work(limit*2+1) and work(limit*3+1) are set to zero, work(1) is set to a and work(limit+1) to b.
dimensioning parameters limit - integer ( kind = 4 ) dimensioning parameter for iwork limit determines the maximum number of subintervals in the partition of the given integration interval (a,b), limit.ge.1. if limit.lt.1, the routine will end with ier = 6.
lenw - integer ( kind = 4 ) dimensioning parameter for work lenw must be at least limit*4. if lenw.lt.limit*4, the routine will end with ier = 6.
last - integer ( kind = 4 ) on return, last equals the number of subintervals produced in the subdivision process, detemines the number of significant elements actually in the work arrays.
work arrays iwork - integer ( kind = 4 ) vector of dimension at least limit, the first k elements of which contain pointers to the error estimates over the subintervals such that work(limit*3+iwork(1)),... , work(limit*3+iwork(k)) form a decreasing sequence, with k = last if last.le.(limit/2+2), and k = limit+1-last otherwise
work - real ( kind = 8 ) vector of dimension at least lenw on return work(1), ..., work(last) contain the left end-points of the subintervals in the partition of (a,b), work(limit+1), ..., work(limit+last) contain the right end-points, work(limit*2+1), ..., work(limit*2+last) contain the integral approximations over the subintervals, work(limit*3+1), ..., work(limit*3+last) contain the error estimates.
Definition at line 2878 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqagse | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | epsabs, | ||
| real ( kind = 8 ) | epsrel, | ||
| integer ( kind = 4 ) | limit, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| integer ( kind = 4 ) | neval, | ||
| integer ( kind = 4 ) | ier, | ||
| real ( kind = 8 ), dimension(limit) | alist, | ||
| real ( kind = 8 ), dimension(limit) | blist, | ||
| real ( kind = 8 ), dimension(limit) | rlist, | ||
| real ( kind = 8 ), dimension(limit) | elist, | ||
| integer ( kind = 4 ), dimension(limit) | iord, | ||
| integer ( kind = 4 ) | last | ||
| ) |
DQAGSE estimates the integral of a function.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose the routine calculates an approximation result to a given definite integral i = integral of f over (a,b), hopefully satisfying following claim for accuracy abs(i-result).le.max(epsabs,epsrel*abs(i)).
Parameters:
on entry f - real ( kind = 8 ) function subprogram defining the integrand function f(x). the actual name for f needs to be declared e x t e r n a l in the driver program.
a - real ( kind = 8 ) lower limit of integration
b - real ( kind = 8 ) upper limit of integration
epsabs - real ( kind = 8 ) absolute accuracy requested epsrel - real ( kind = 8 ) relative accuracy requested if epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), the routine will end with ier = 6.
limit - integer ( kind = 4 ) gives an upperbound on the number of subintervals in the partition of (a,b)
on return result - real ( kind = 8 ) approximation to the integral
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
neval - integer ( kind = 4 ) number of integrand evaluations
ier - integer ( kind = 4 ) ier = 0 normal and reliable termination of the routine. it is assumed that the requested accuracy has been achieved. ier.gt.0 abnormal termination of the routine the estimates for integral and error are less reliable. it is assumed that the requested accuracy has not been achieved. error messages = 1 maximum number of subdivisions allowed has been achieved. one can allow more sub- divisions by increasing the value of limit (and taking the according dimension adjustments into account). however, if this yields no improvement it is advised to analyze the integrand in order to determine the integration difficulties. if the position of a local difficulty can be determined (e.g. singularity, discontinuity within the interval) one will probably gain from splitting up the interval at this point and calling the integrator on the subranges. if possible, an appropriate special-purpose integrator should be used, which is designed for handling the type of difficulty involved. = 2 the occurrence of roundoff error is detec- ted, which prevents the requested tolerance from being achieved. the error may be under-estimated. = 3 extremely bad integrand behaviour occurs at some points of the integration interval. = 4 the algorithm does not converge. roundoff error is detected in the extrapolation table. it is presumed that the requested tolerance cannot be achieved, and that the returned result is the best which can be obtained. = 5 the integral is probably divergent, or slowly convergent. it must be noted that divergence can occur with any other value of ier. = 6 the input is invalid, because epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28). result, abserr, neval, last, rlist(1), iord(1) and elist(1) are set to zero. alist(1) and blist(1) are set to a and b respectively.
alist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the left end points of the subintervals in the partition of the given integration range (a,b)
blist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the right end points of the subintervals in the partition of the given integration range (a,b)
rlist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the integral approximations on the subintervals
elist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the moduli of the absolute error estimates on the subintervals
iord - integer ( kind = 4 ) vector of dimension at least limit, the first k elements of which are pointers to the error estimates over the subintervals, such that elist(iord(1)), ..., elist(iord(k)) form a decreasing sequence, with k = last if last.le.(limit/2+2), and k = limit+1-last otherwise
last - integer ( kind = 4 ) number of subintervals actually produced in the subdivision process
Local parameters:
the dimension of rlist2 is determined by the value of limexp in routine dqelg (rlist2 should be of dimension (limexp+2) at least). list of major variables
alist - list of left end points of all subintervals considered up to now blist - list of right end points of all subintervals considered up to now rlist(i) - approximation to the integral over (alist(i),blist(i)) rlist2 - array of dimension at least limexp+2 containing the part of the epsilon table which is still needed for further computations elist(i) - error estimate applying to rlist(i) maxerr - pointer to the interval with largest error estimate errmax - elist(maxerr) erlast - error on the interval currently subdivided (before that subdivision has taken place) area - sum of the integrals over the subintervals errsum - sum of the errors over the subintervals errbnd - requested accuracy max(epsabs,epsrel* abs(result)) *****1 - variable for the left interval *****2 - variable for the right interval last - index for subdivision nres - number of calls to the extrapolation routine numrl2 - number of elements currently in rlist2. if an appropriate approximation to the compounded integral has been obtained it is put in rlist2(numrl2) after numrl2 has been increased by one. small - length of the smallest interval considered up to now, multiplied by 1.5 erlarg - sum of the errors over the intervals larger than the smallest interval considered up to now extrap - logical variable denoting that the routine is attempting to perform extrapolation i.e. before subdividing the smallest interval we try to decrease the value of erlarg. noext - logical variable denoting that extrapolation is no longer allowed (true value)
machine dependent constants
epmach is the largest relative spacing. uflow is the smallest positive magnitude. oflow is the largest positive magnitude.
Definition at line 2492 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqawc | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | c, | ||
| real ( kind = 8 ) | epsabs, | ||
| real ( kind = 8 ) | epsrel, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| integer ( kind = 4 ) | neval, | ||
| integer ( kind = 4 ) | ier, | ||
| integer ( kind = 4 ) | limit, | ||
| integer ( kind = 4 ) | lenw, | ||
| integer ( kind = 4 ) | last, | ||
| integer ( kind = 4 ), dimension(limit) | iwork, | ||
| real ( kind = 8 ), dimension(lenw) | work | ||
| ) |
DQAWC computes a Cauchy principal value.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose the routine calculates an approximation result to a cauchy principal value i = integral of f*w over (a,b) (w(x) = 1/((x-c), c.ne.a, c.ne.b), hopefully satisfying following claim for accuracy abs(i-result).le.max(epsabe,epsrel*abs(i)).
Parameters:
on entry f - real ( kind = 8 ) function subprogram defining the integrand function f(x). the actual name for f needs to be declared e x t e r n a l in the driver program.
a - real ( kind = 8 ) under limit of integration
b - real ( kind = 8 ) upper limit of integration
c - parameter in the weight function, c.ne.a, c.ne.b. if c = a or c = b, the routine will end with ier = 6 .
epsabs - real ( kind = 8 ) absolute accuracy requested epsrel - real ( kind = 8 ) relative accuracy requested if epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), the routine will end with ier = 6.
on return result - real ( kind = 8 ) approximation to the integral
abserr - real ( kind = 8 ) estimate or the modulus of the absolute error, which should equal or exceed abs(i-result)
neval - integer ( kind = 4 ) number of integrand evaluations
ier - integer ( kind = 4 ) ier = 0 normal and reliable termination of the routine. it is assumed that the requested accuracy has been achieved. ier.gt.0 abnormal termination of the routine the estimates for integral and error are less reliable. it is assumed that the requested accuracy has not been achieved. error messages ier = 1 maximum number of subdivisions allowed has been achieved. one can allow more sub- divisions by increasing the value of limit (and taking the according dimension adjustments into account). however, if this yields no improvement it is advised to analyze the integrand in order to determine the integration difficulties. if the position of a local difficulty can be determined (e.g. singularity, discontinuity within the interval) one will probably gain from splitting up the interval at this point and calling appropriate integrators on the subranges. = 2 the occurrence of roundoff error is detec- ted, which prevents the requested tolerance from being achieved. = 3 extremely bad integrand behaviour occurs at some points of the integration interval. = 6 the input is invalid, because c = a or c = b or (epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) or limit.lt.1 or lenw.lt.limit*4. result, abserr, neval, last are set to zero. exept when lenw or limit is invalid, iwork(1), work(limit*2+1) and work(limit*3+1) are set to zero, work(1) is set to a and work(limit+1) to b.
dimensioning parameters limit - integer ( kind = 4 ) dimensioning parameter for iwork limit determines the maximum number of subintervals in the partition of the given integration interval (a,b), limit.ge.1. if limit.lt.1, the routine will end with ier = 6.
lenw - integer ( kind = 4 ) dimensioning parameter for work lenw must be at least limit*4. if lenw.lt.limit*4, the routine will end with ier = 6.
last - integer ( kind = 4 ) on return, last equals the number of subintervals produced in the subdivision process, which determines the number of significant elements actually in the work arrays.
work arrays iwork - integer ( kind = 4 ) vector of dimension at least limit, the first k elements of which contain pointers to the error estimates over the subintervals, such that work(limit*3+iwork(1)), ... , work(limit*3+iwork(k)) form a decreasing sequence, with k = last if last.le.(limit/2+2), and k = limit+1-last otherwise
work - real ( kind = 8 ) vector of dimension at least lenw on return work(1), ..., work(last) contain the left end points of the subintervals in the partition of (a,b), work(limit+1), ..., work(limit+last) contain the right end points, work(limit*2+1), ..., work(limit*2+last) contain the integral approximations over the subintervals, work(limit*3+1), ..., work(limit*3+last) contain the error estimates.
Definition at line 3384 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqawce | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | c, | ||
| real ( kind = 8 ) | epsabs, | ||
| real ( kind = 8 ) | epsrel, | ||
| integer ( kind = 4 ) | limit, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| integer ( kind = 4 ) | neval, | ||
| integer ( kind = 4 ) | ier, | ||
| real ( kind = 8 ), dimension(limit) | alist, | ||
| real ( kind = 8 ), dimension(limit) | blist, | ||
| real ( kind = 8 ), dimension(limit) | rlist, | ||
| real ( kind = 8 ), dimension(limit) | elist, | ||
| integer ( kind = 4 ), dimension(limit) | iord, | ||
| integer ( kind = 4 ) | last | ||
| ) |
DQAWCE computes a Cauchy principal value.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
*** purpose the routine calculates an approximation result to a cauchy principal value i = integral of f*w over (a,b) (w(x) = 1/(x-c), (c.ne.a, c.ne.b), hopefully satisfying following claim for accuracy abs(i-result).le.max(epsabs,epsrel*abs(i))
Parameters:
on entry f - real ( kind = 8 ) function subprogram defining the integrand function f(x). the actual name for f needs to be declared e x t e r n a l in the driver program.
a - real ( kind = 8 ) lower limit of integration
b - real ( kind = 8 ) upper limit of integration
c - real ( kind = 8 ) parameter in the weight function, c.ne.a, c.ne.b if c = a or c = b, the routine will end with ier = 6.
epsabs - real ( kind = 8 ) absolute accuracy requested epsrel - real ( kind = 8 ) relative accuracy requested if epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), the routine will end with ier = 6.
limit - integer ( kind = 4 ) gives an upper bound on the number of subintervals in the partition of (a,b), limit.ge.1
on return result - real ( kind = 8 ) approximation to the integral
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
neval - integer ( kind = 4 ) number of integrand evaluations
ier - integer ( kind = 4 ) ier = 0 normal and reliable termination of the routine. it is assumed that the requested accuracy has been achieved. ier.gt.0 abnormal termination of the routine the estimates for integral and error are less reliable. it is assumed that the requested accuracy has not been achieved. error messages ier = 1 maximum number of subdivisions allowed has been achieved. one can allow more sub- divisions by increasing the value of limit. however, if this yields no improvement it is advised to analyze the the integrand, in order to determine the the integration difficulties. if the position of a local difficulty can be determined (e.g. singularity, discontinuity within the interval) one will probably gain from splitting up the interval at this point and calling appropriate integrators on the subranges. = 2 the occurrence of roundoff error is detec- ted, which prevents the requested tolerance from being achieved. = 3 extremely bad integrand behaviour occurs at some interior points of the integration interval. = 6 the input is invalid, because c = a or c = b or (epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) or limit.lt.1. result, abserr, neval, rlist(1), elist(1), iord(1) and last are set to zero. alist(1) and blist(1) are set to a and b respectively.
alist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the left end points of the subintervals in the partition of the given integration range (a,b)
blist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the right end points of the subintervals in the partition of the given integration range (a,b)
rlist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the integral approximations on the subintervals
elist - real ( kind = 8 ) vector of dimension limit, the first last elements of which are the moduli of the absolute error estimates on the subintervals
iord - integer ( kind = 4 ) vector of dimension at least limit, the first k elements of which are pointers to the error estimates over the subintervals, so that elist(iord(1)), ..., elist(iord(k)) with k = last if last.le.(limit/2+2), and k = limit+1-last otherwise, form a decreasing sequence
last - integer ( kind = 4 ) number of subintervals actually produced in the subdivision process
Local Parameters:
alist - list of left end points of all subintervals considered up to now blist - list of right end points of all subintervals considered up to now rlist(i) - approximation to the integral over (alist(i),blist(i)) elist(i) - error estimate applying to rlist(i) maxerr - pointer to the interval with largest error estimate errmax - elist(maxerr) area - sum of the integrals over the subintervals errsum - sum of the errors over the subintervals errbnd - requested accuracy max(epsabs,epsrel* abs(result)) *****1 - variable for the left subinterval *****2 - variable for the right subinterval last - index for subdivision
machine dependent constants epmach is the largest relative spacing. uflow is the smallest positive magnitude.
Definition at line 3073 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqawf | ( | real ( kind = 8 ), external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | omega, | ||
| integer ( kind = 4 ) | integr, | ||
| real ( kind = 8 ) | epsabs, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| integer ( kind = 4 ) | neval, | ||
| integer ( kind = 4 ) | ier, | ||
| integer ( kind = 4 ) | limlst, | ||
| integer ( kind = 4 ) | lst, | ||
| integer ( kind = 4 ) | leniw, | ||
| integer ( kind = 4 ) | maxp1, | ||
| integer ( kind = 4 ) | lenw, | ||
| integer ( kind = 4 ), dimension(leniw) | iwork, | ||
| real ( kind = 8 ), dimension(lenw) | work | ||
| ) |
DQAWF computes Fourier integrals over the interval [ A, +Infinity ).
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose the routine calculates an approximation result to a given fourier integral i=integral of f(x)*w(x) over (a,infinity) where w(x) = cos(omega*x) or w(x) = sin(omega*x). hopefully satisfying following claim for accuracy abs(i-result).le.epsabs.
Parameters:
on entry f - real ( kind = 8 ) function subprogram defining the integrand function f(x). the actual name for f needs to be declared e x t e r n a l in the driver program.
a - real ( kind = 8 ) lower limit of integration
omega - real ( kind = 8 ) parameter in the integrand weight function
integr - integer ( kind = 4 ) indicates which of the weight functions is used integr = 1 w(x) = cos(omega*x) integr = 2 w(x) = sin(omega*x) if integr.ne.1.and.integr.ne.2, the routine will end with ier = 6.
epsabs - real ( kind = 8 ) absolute accuracy requested, epsabs.gt.0. if epsabs.le.0, the routine will end with ier = 6.
on return result - real ( kind = 8 ) approximation to the integral
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
neval - integer ( kind = 4 ) number of integrand evaluations
ier - integer ( kind = 4 ) ier = 0 normal and reliable termination of the routine. it is assumed that the requested accuracy has been achieved. ier.gt.0 abnormal termination of the routine. the estimates for integral and error are less reliable. it is assumed that the requested accuracy has not been achieved. error messages if omega.ne.0 ier = 1 maximum number of cycles allowed has been achieved, i.e. of subintervals (a+(k-1)c,a+kc) where c = (2*int(abs(omega))+1)*pi/abs(omega), for k = 1, 2, ..., lst. one can allow more cycles by increasing the value of limlst (and taking the according dimension adjustments into account). examine the array iwork which contains the error flags on the cycles, in order to look for eventual local integration difficulties. if the position of a local difficulty can be determined (e.g. singularity, discontinuity within the interval) one will probably gain from splitting up the interval at this point and calling appropriate integrators on the subranges. = 4 the extrapolation table constructed for convergence accelaration of the series formed by the integral contributions over the cycles, does not converge to within the requested accuracy. as in the case of ier = 1, it is advised to examine the array iwork which contains the error flags on the cycles. = 6 the input is invalid because (integr.ne.1 and integr.ne.2) or epsabs.le.0 or limlst.lt.1 or leniw.lt.(limlst+2) or maxp1.lt.1 or lenw.lt.(leniw*2+maxp1*25). result, abserr, neval, lst are set to zero. = 7 bad integrand behaviour occurs within one or more of the cycles. location and type of the difficulty involved can be determined from the first lst elements of vector iwork. here lst is the number of cycles actually needed (see below). iwork(k) = 1 the maximum number of subdivisions (=(leniw-limlst) /2) has been achieved on the k th cycle. = 2 occurrence of roundoff error is detected and prevents the tolerance imposed on the k th cycle, from being achieved on this cycle. = 3 extremely bad integrand behaviour occurs at some points of the k th cycle. = 4 the integration procedure over the k th cycle does not converge (to within the required accuracy) due to roundoff in the extrapolation procedure invoked on this cycle. it is assumed that the result on this interval is the best which can be obtained. = 5 the integral over the k th cycle is probably divergent or slowly convergent. it must be noted that divergence can occur with any other value of iwork(k). if omega = 0 and integr = 1, the integral is calculated by means of dqagie, and ier = iwork(1) (with meaning as described for iwork(k),k = 1).
dimensioning parameters limlst - integer ( kind = 4 ) limlst gives an upper bound on the number of cycles, limlst.ge.3. if limlst.lt.3, the routine will end with ier = 6.
lst - integer ( kind = 4 ) on return, lst indicates the number of cycles actually needed for the integration. if omega = 0, then lst is set to 1.
leniw - integer ( kind = 4 ) dimensioning parameter for iwork. on entry, (leniw-limlst)/2 equals the maximum number of subintervals allowed in the partition of each cycle, leniw.ge.(limlst+2). if leniw.lt.(limlst+2), the routine will end with ier = 6.
maxp1 - integer ( kind = 4 ) maxp1 gives an upper bound on the number of chebyshev moments which can be stored, i.e. for the intervals of lengths abs(b-a)*2**(-l), l = 0,1, ..., maxp1-2, maxp1.ge.1. if maxp1.lt.1, the routine will end with ier = 6. lenw - integer ( kind = 4 ) dimensioning parameter for work lenw must be at least leniw*2+maxp1*25. if lenw.lt.(leniw*2+maxp1*25), the routine will end with ier = 6.
work arrays iwork - integer ( kind = 4 ) vector of dimension at least leniw on return, iwork(k) for k = 1, 2, ..., lst contain the error flags on the cycles.
work - real ( kind = 8 ) vector of dimension at least lenw on return, work(1), ..., work(lst) contain the integral approximations over the cycles, work(limlst+1), ..., work(limlst+lst) contain the error extimates over the cycles. further elements of work have no specific meaning for the user.
Definition at line 3969 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqawfe | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | omega, | ||
| integer ( kind = 4 ) | integr, | ||
| real ( kind = 8 ) | epsabs, | ||
| integer ( kind = 4 ) | limlst, | ||
| integer ( kind = 4 ) | limit, | ||
| integer ( kind = 4 ) | maxp1, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| integer ( kind = 4 ) | neval, | ||
| integer ( kind = 4 ) | ier, | ||
| real ( kind = 8 ), dimension(limlst) | rslst, | ||
| real ( kind = 8 ), dimension(limlst) | erlst, | ||
| integer ( kind = 4 ), dimension(limlst) | ierlst, | ||
| integer ( kind = 4 ) | lst, | ||
| real ( kind = 8 ), dimension(limit) | alist, | ||
| real ( kind = 8 ), dimension(limit) | blist, | ||
| real ( kind = 8 ), dimension(limit) | rlist, | ||
| real ( kind = 8 ), dimension(limit) | elist, | ||
| integer ( kind = 4 ), dimension(limit) | iord, | ||
| integer ( kind = 4 ), dimension(limit) | nnlog, | ||
| real ( kind = 8 ), dimension(maxp1,25) | chebmo | ||
| ) |
DQAWFE computes Fourier integrals.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose the routine calculates an approximation result to a given fourier integal i = integral of f(x)*w(x) over (a,infinity) where w(x)=cos(omega*x) or w(x)=sin(omega*x), hopefully satisfying following claim for accuracy abs(i-result).le.epsabs.
Parameters:
on entry f - real ( kind = 8 ) function subprogram defining the integrand function f(x). the actual name for f needs to be declared e x t e r n a l in the driver program.
a - real ( kind = 8 ) lower limit of integration
omega - real ( kind = 8 ) parameter in the weight function
integr - integer ( kind = 4 ) indicates which weight function is used integr = 1 w(x) = cos(omega*x) integr = 2 w(x) = sin(omega*x) if integr.ne.1.and.integr.ne.2, the routine will end with ier = 6.
epsabs - real ( kind = 8 ) absolute accuracy requested, epsabs.gt.0 if epsabs.le.0, the routine will end with ier = 6.
limlst - integer ( kind = 4 ) limlst gives an upper bound on the number of cycles, limlst.ge.1. if limlst.lt.3, the routine will end with ier = 6.
limit - integer ( kind = 4 ) gives an upper bound on the number of subintervals allowed in the partition of each cycle, limit.ge.1 each cycle, limit.ge.1.
maxp1 - integer ( kind = 4 ) gives an upper bound on the number of chebyshev moments which can be stored, i.e. for the intervals of lengths abs(b-a)*2**(-l), l=0,1, ..., maxp1-2, maxp1.ge.1
on return result - real ( kind = 8 ) approximation to the integral x
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
neval - integer ( kind = 4 ) number of integrand evaluations
ier - ier = 0 normal and reliable termination of the routine. it is assumed that the requested accuracy has been achieved. ier.gt.0 abnormal termination of the routine. the estimates for integral and error are less reliable. it is assumed that the requested accuracy has not been achieved. error messages if omega.ne.0 ier = 1 maximum number of cycles allowed has been achieved., i.e. of subintervals (a+(k-1)c,a+kc) where c = (2*int(abs(omega))+1)*pi/abs(omega), for k = 1, 2, ..., lst. one can allow more cycles by increasing the value of limlst (and taking the according dimension adjustments into account). examine the array iwork which contains the error flags on the cycles, in order to look for eventual local integration difficulties. if the position of a local difficulty can be determined (e.g. singularity, discontinuity within the interval) one will probably gain from splitting up the interval at this point and calling appropriate integrators on the subranges. = 4 the extrapolation table constructed for convergence acceleration of the series formed by the integral contributions over the cycles, does not converge to within the requested accuracy. as in the case of ier = 1, it is advised to examine the array iwork which contains the error flags on the cycles. = 6 the input is invalid because (integr.ne.1 and integr.ne.2) or epsabs.le.0 or limlst.lt.3. result, abserr, neval, lst are set to zero. = 7 bad integrand behaviour occurs within one or more of the cycles. location and type of the difficulty involved can be determined from the vector ierlst. here lst is the number of cycles actually needed (see below). ierlst(k) = 1 the maximum number of subdivisions (= limit) has been achieved on the k th cycle. = 2 occurrence of roundoff error is detected and prevents the tolerance imposed on the k th cycle, from being achieved. = 3 extremely bad integrand behaviour occurs at some points of the k th cycle. = 4 the integration procedure over the k th cycle does not converge (to within the required accuracy) due to roundoff in the extrapolation procedure invoked on this cycle. it is assumed that the result on this interval is the best which can be obtained. = 5 the integral over the k th cycle is probably divergent or slowly convergent. it must be noted that divergence can occur with any other value of ierlst(k). if omega = 0 and integr = 1, the integral is calculated by means of dqagie and ier = ierlst(1) (with meaning as described for ierlst(k), k = 1).
rslst - real ( kind = 8 ) vector of dimension at least limlst rslst(k) contains the integral contribution over the interval (a+(k-1)c,a+kc) where c = (2*int(abs(omega))+1)*pi/abs(omega), k = 1, 2, ..., lst. note that, if omega = 0, rslst(1) contains the value of the integral over (a,infinity).
erlst - real ( kind = 8 ) vector of dimension at least limlst erlst(k) contains the error estimate corresponding with rslst(k).
ierlst - integer ( kind = 4 ) vector of dimension at least limlst ierlst(k) contains the error flag corresponding with rslst(k). for the meaning of the local error flags see description of output parameter ier.
lst - integer ( kind = 4 ) number of subintervals needed for the integration if omega = 0 then lst is set to 1.
alist, blist, rlist, elist - real ( kind = 8 ) vector of dimension at least limit,
iord, nnlog - integer ( kind = 4 ) vector of dimension at least limit, providing space for the quantities needed in the subdivision process of each cycle
chebmo - real ( kind = 8 ) array of dimension at least (maxp1,25), providing space for the chebyshev moments needed within the cycles
Local Parameters:
the dimension of psum is determined by the value of limexp in routine dqelg (psum must be of dimension (limexp+2) at least).
c1, c2 - end points of subinterval (of length cycle) cycle - (2*int(abs(omega))+1)*pi/abs(omega) psum - vector of dimension at least (limexp+2) (see routine dqelg) psum contains the part of the epsilon table which is still needed for further computations. each element of psum is a partial sum of the series which should sum to the value of the integral. errsum - sum of error estimates over the subintervals, calculated cumulatively epsa - absolute tolerance requested over current subinterval chebmo - array containing the modified chebyshev moments (see also routine dqc25f)
Definition at line 3635 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqawo | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | omega, | ||
| integer ( kind = 4 ) | integr, | ||
| real ( kind = 8 ) | epsabs, | ||
| real ( kind = 8 ) | epsrel, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| integer ( kind = 4 ) | neval, | ||
| integer ( kind = 4 ) | ier, | ||
| integer ( kind = 4 ) | leniw, | ||
| integer ( kind = 4 ) | maxp1, | ||
| integer ( kind = 4 ) | lenw, | ||
| integer ( kind = 4 ) | last, | ||
| integer ( kind = 4 ), dimension(leniw) | iwork, | ||
| real ( kind = 8 ), dimension(lenw) | work | ||
| ) |
DQAWO computes the integrals of oscillatory integrands.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose the routine calculates an approximation result to a given definite integral i=integral of f(x)*w(x) over (a,b) where w(x) = cos(omega*x) or w(x) = sin(omega*x), hopefully satisfying following claim for accuracy abs(i-result).le.max(epsabs,epsrel*abs(i)).
Parameters:
on entry f - real ( kind = 8 ) function subprogram defining the function f(x). the actual name for f needs to be declared e x t e r n a l in the driver program.
a - real ( kind = 8 ) lower limit of integration
b - real ( kind = 8 ) upper limit of integration
omega - real ( kind = 8 ) parameter in the integrand weight function
integr - integer ( kind = 4 ) indicates which of the weight functions is used integr = 1 w(x) = cos(omega*x) integr = 2 w(x) = sin(omega*x) if integr.ne.1.and.integr.ne.2, the routine will end with ier = 6.
epsabs - real ( kind = 8 ) absolute accuracy requested epsrel - real ( kind = 8 ) relative accuracy requested if epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), the routine will end with ier = 6.
on return result - real ( kind = 8 ) approximation to the integral
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
neval - integer ( kind = 4 ) number of integrand evaluations
ier - integer ( kind = 4 ) ier = 0 normal and reliable termination of the routine. it is assumed that the requested accuracy has been achieved.
dimensioning parameters leniw - integer ( kind = 4 ) dimensioning parameter for iwork. leniw/2 equals the maximum number of subintervals allowed in the partition of the given integration interval (a,b), leniw.ge.2. if leniw.lt.2, the routine will end with ier = 6.
maxp1 - integer ( kind = 4 ) gives an upper bound on the number of chebyshev moments which can be stored, i.e. for the intervals of lengths abs(b-a)*2**(-l), l=0,1, ..., maxp1-2, maxp1.ge.1 if maxp1.lt.1, the routine will end with ier = 6.
lenw - integer ( kind = 4 ) dimensioning parameter for work lenw must be at least leniw*2+maxp1*25. if lenw.lt.(leniw*2+maxp1*25), the routine will end with ier = 6.
last - integer ( kind = 4 ) on return, last equals the number of subintervals produced in the subdivision process, which determines the number of significant elements actually in the work arrays.
work arrays iwork - integer ( kind = 4 ) vector of dimension at least leniw on return, the first k elements of which contain pointers to the error estimates over the subintervals, such that work(limit*3+iwork(1)), .. work(limit*3+iwork(k)) form a decreasing sequence, with limit = lenw/2 , and k = last if last.le.(limit/2+2), and k = limit+1-last otherwise. furthermore, iwork(limit+1), ..., iwork(limit+ last) indicate the subdivision levels of the subintervals, such that iwork(limit+i) = l means that the subinterval numbered i is of length abs(b-a)*2**(1-l).
work - real ( kind = 8 ) vector of dimension at least lenw on return work(1), ..., work(last) contain the left end points of the subintervals in the partition of (a,b), work(limit+1), ..., work(limit+last) contain the right end points, work(limit*2+1), ..., work(limit*2+last) contain the integral approximations over the subintervals, work(limit*3+1), ..., work(limit*3+last) contain the error estimates. work(limit*4+1), ..., work(limit*4+maxp1*25) provide space for storing the chebyshev moments. note that limit = lenw/2.
Definition at line 4738 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqawoe | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | omega, | ||
| integer ( kind = 4 ) | integr, | ||
| real ( kind = 8 ) | epsabs, | ||
| real ( kind = 8 ) | epsrel, | ||
| integer ( kind = 4 ) | limit, | ||
| integer ( kind = 4 ) | icall, | ||
| integer ( kind = 4 ) | maxp1, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| integer ( kind = 4 ) | neval, | ||
| integer ( kind = 4 ) | ier, | ||
| integer ( kind = 4 ) | last, | ||
| real ( kind = 8 ), dimension(limit) | alist, | ||
| real ( kind = 8 ), dimension(limit) | blist, | ||
| real ( kind = 8 ), dimension(limit) | rlist, | ||
| real ( kind = 8 ), dimension(limit) | elist, | ||
| integer ( kind = 4 ), dimension(limit) | iord, | ||
| integer ( kind = 4 ), dimension(limit) | nnlog, | ||
| integer ( kind = 4 ) | momcom, | ||
| real ( kind = 8 ), dimension(maxp1,25) | chebmo | ||
| ) |
DQAWOE computes the integrals of oscillatory integrands.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose the routine calculates an approximation result to a given definite integral i = integral of f(x)*w(x) over (a,b) where w(x) = cos(omega*x) or w(x)=sin(omega*x), hopefully satisfying following claim for accuracy abs(i-result).le.max(epsabs,epsrel*abs(i)).
Parameters:
on entry f - real ( kind = 8 ) function subprogram defining the integrand function f(x). the actual name for f needs to be declared e x t e r n a l in the driver program.
a - real ( kind = 8 ) lower limit of integration
b - real ( kind = 8 ) upper limit of integration
omega - real ( kind = 8 ) parameter in the integrand weight function
integr - integer ( kind = 4 ) indicates which of the weight functions is to be used integr = 1 w(x) = cos(omega*x) integr = 2 w(x) = sin(omega*x) if integr.ne.1 and integr.ne.2, the routine will end with ier = 6.
epsabs - real ( kind = 8 ) absolute accuracy requested epsrel - real ( kind = 8 ) relative accuracy requested if epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), the routine will end with ier = 6.
limit - integer ( kind = 4 ) gives an upper bound on the number of subdivisions in the partition of (a,b), limit.ge.1.
icall - integer ( kind = 4 ) if dqawoe is to be used only once, icall must be set to 1. assume that during this call, the chebyshev moments (for clenshaw-curtis integration of degree 24) have been computed for intervals of lenghts (abs(b-a))*2**(-l), l=0,1,2,...momcom-1. if icall.gt.1 this means that dqawoe has been called twice or more on intervals of the same length abs(b-a). the chebyshev moments already computed are then re-used in subsequent calls. if icall.lt.1, the routine will end with ier = 6.
maxp1 - integer ( kind = 4 ) gives an upper bound on the number of chebyshev moments which can be stored, i.e. for the intervals of lenghts abs(b-a)*2**(-l), l=0,1, ..., maxp1-2, maxp1.ge.1. if maxp1.lt.1, the routine will end with ier = 6.
on return result - real ( kind = 8 ) approximation to the integral
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
neval - integer ( kind = 4 ) number of integrand evaluations
ier - integer ( kind = 4 ) ier = 0 normal and reliable termination of the routine. it is assumed that the requested accuracy has been achieved.
last - integer ( kind = 4 ) on return, last equals the number of subintervals produces in the subdivision process, which determines the number of significant elements actually in the work arrays. alist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the left end points of the subintervals in the partition of the given integration range (a,b)
blist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the right end points of the subintervals in the partition of the given integration range (a,b)
rlist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the integral approximations on the subintervals
elist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the moduli of the absolute error estimates on the subintervals
iord - integer ( kind = 4 ) vector of dimension at least limit, the first k elements of which are pointers to the error estimates over the subintervals, such that elist(iord(1)), ..., elist(iord(k)) form a decreasing sequence, with k = last if last.le.(limit/2+2), and k = limit+1-last otherwise.
nnlog - integer ( kind = 4 ) vector of dimension at least limit, containing the subdivision levels of the subintervals, i.e. iwork(i) = l means that the subinterval numbered i is of length abs(b-a)*2**(1-l)
on entry and return momcom - integer ( kind = 4 ) indicating that the chebyshev moments have been computed for intervals of lengths (abs(b-a))*2**(-l), l=0,1,2, ..., momcom-1, momcom.lt.maxp1
chebmo - real ( kind = 8 ) array of dimension (maxp1,25) containing the chebyshev moments
Local Parameters:
the dimension of rlist2 is determined by the value of limexp in routine dqelg (rlist2 should be of dimension (limexp+2) at least). list of major variables
alist - list of left end points of all subintervals considered up to now blist - list of right end points of all subintervals considered up to now rlist(i) - approximation to the integral over (alist(i),blist(i)) rlist2 - array of dimension at least limexp+2 containing the part of the epsilon table which is still needed for further computations elist(i) - error estimate applying to rlist(i) maxerr - pointer to the interval with largest error estimate errmax - elist(maxerr) erlast - error on the interval currently subdivided area - sum of the integrals over the subintervals errsum - sum of the errors over the subintervals errbnd - requested accuracy max(epsabs,epsrel* abs(result)) *****1 - variable for the left subinterval *****2 - variable for the right subinterval last - index for subdivision nres - number of calls to the extrapolation routine numrl2 - number of elements in rlist2. if an appropriate approximation to the compounded integral has been obtained it is put in rlist2(numrl2) after numrl2 has been increased by one small - length of the smallest interval considered up to now, multiplied by 1.5 erlarg - sum of the errors over the intervals larger than the smallest interval considered up to now extrap - logical variable denoting that the routine is attempting to perform extrapolation, i.e. before subdividing the smallest interval we try to decrease the value of erlarg noext - logical variable denoting that extrapolation is no longer allowed (true value)
machine dependent constants
epmach is the largest relative spacing. uflow is the smallest positive magnitude. oflow is the largest positive magnitude.
Definition at line 4283 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqaws | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | alfa, | ||
| real ( kind = 8 ) | beta, | ||
| integer ( kind = 4 ) | integr, | ||
| real ( kind = 8 ) | epsabs, | ||
| real ( kind = 8 ) | epsrel, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| integer ( kind = 4 ) | neval, | ||
| integer ( kind = 4 ) | ier, | ||
| integer ( kind = 4 ) | limit, | ||
| integer ( kind = 4 ) | lenw, | ||
| integer ( kind = 4 ) | last, | ||
| integer ( kind = 4 ), dimension(limit) | iwork, | ||
| real ( kind = 8 ), dimension(lenw) | work | ||
| ) |
DQAWS estimates integrals with algebraico-logarithmic endpoint singularities.
Modified:
12 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose the routine calculates an approximation result to a given definite integral i = integral of f*w over (a,b), (where w shows a singular behaviour at the end points see parameter integr). hopefully satisfying following claim for accuracy abs(i-result).le.max(epsabs,epsrel*abs(i)).
Parameters:
on entry f - real ( kind = 8 ) function subprogram defining the integrand function f(x). the actual name for f needs to be declared e x t e r n a l in the driver program.
a - real ( kind = 8 ) lower limit of integration
b - real ( kind = 8 ) upper limit of integration, b.gt.a if b.le.a, the routine will end with ier = 6.
alfa - real ( kind = 8 ) parameter in the integrand function, alfa.gt.(-1) if alfa.le.(-1), the routine will end with ier = 6.
beta - real ( kind = 8 ) parameter in the integrand function, beta.gt.(-1) if beta.le.(-1), the routine will end with ier = 6.
integr - integer ( kind = 4 ) indicates which weight function is to be used = 1 (x-a)**alfa*(b-x)**beta = 2 (x-a)**alfa*(b-x)**beta*log(x-a) = 3 (x-a)**alfa*(b-x)**beta*log(b-x) = 4 (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x) if integr.lt.1 or integr.gt.4, the routine will end with ier = 6.
epsabs - real ( kind = 8 ) absolute accuracy requested epsrel - real ( kind = 8 ) relative accuracy requested if epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), the routine will end with ier = 6.
on return result - real ( kind = 8 ) approximation to the integral
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
neval - integer ( kind = 4 ) number of integrand evaluations
ier - integer ( kind = 4 ) ier = 0 normal and reliable termination of the routine. it is assumed that the requested accuracy has been achieved. ier.gt.0 abnormal termination of the routine the estimates for the integral and error are less reliable. it is assumed that the requested accuracy has not been achieved. error messages ier = 1 maximum number of subdivisions allowed has been achieved. one can allow more subdivisions by increasing the value of limit (and taking the according dimension adjustments into account). however, if this yields no improvement it is advised to analyze the integrand, in order to determine the integration difficulties which prevent the requested tolerance from being achieved. in case of a jump discontinuity or a local singularity of algebraico-logarithmic type at one or more interior points of the integration range, one should proceed by splitting up the interval at these points and calling the integrator on the subranges. = 2 the occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. = 3 extremely bad integrand behaviour occurs at some points of the integration interval. = 6 the input is invalid, because b.le.a or alfa.le.(-1) or beta.le.(-1) or or integr.lt.1 or integr.gt.4 or (epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) or limit.lt.2 or lenw.lt.limit*4. result, abserr, neval, last are set to zero. except when lenw or limit is invalid iwork(1), work(limit*2+1) and work(limit*3+1) are set to zero, work(1) is set to a and work(limit+1) to b.
dimensioning parameters limit - integer ( kind = 4 ) dimensioning parameter for iwork limit determines the maximum number of subintervals in the partition of the given integration interval (a,b), limit.ge.2. if limit.lt.2, the routine will end with ier = 6.
lenw - integer ( kind = 4 ) dimensioning parameter for work lenw must be at least limit*4. if lenw.lt.limit*4, the routine will end with ier = 6.
last - integer ( kind = 4 ) on return, last equals the number of subintervals produced in the subdivision process, which determines the significant number of elements actually in the work arrays.
work arrays iwork - integer ( kind = 4 ) vector of dimension limit, the first k elements of which contain pointers to the error estimates over the subintervals, such that work(limit*3+iwork(1)), ..., work(limit*3+iwork(k)) form a decreasing sequence with k = last if last.le.(limit/2+2), and k = limit+1-last otherwise
work - real ( kind = 8 ) vector of dimension lenw on return work(1), ..., work(last) contain the left end points of the subintervals in the partition of (a,b), work(limit+1), ..., work(limit+last) contain the right end points, work(limit*2+1), ..., work(limit*2+last) contain the integral approximations over the subintervals, work(limit*3+1), ..., work(limit*3+last) contain the error estimates.
Definition at line 5318 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqawse | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | alfa, | ||
| real ( kind = 8 ) | beta, | ||
| integer ( kind = 4 ) | integr, | ||
| real ( kind = 8 ) | epsabs, | ||
| real ( kind = 8 ) | epsrel, | ||
| integer ( kind = 4 ), dimension( kind = 4 ) | limit, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| integer ( kind = 4 ) | neval, | ||
| integer ( kind = 4 ) | ier, | ||
| real ( kind = 8 ), dimension(limit) | alist, | ||
| real ( kind = 8 ), dimension(limit) | blist, | ||
| real ( kind = 8 ), dimension(limit) | rlist, | ||
| real ( kind = 8 ), dimension(limit) | elist, | ||
| integer ( kind = 4 ), dimension(limit) | iord, | ||
| integer ( kind = 4 ) | last | ||
| ) |
DQAWSE estimates integrals with algebraico-logarithmic end singularities.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose the routine calculates an approximation result to a given definite integral i = integral of f*w over (a,b), (where w shows a singular behaviour at the end points, see parameter integr). hopefully satisfying following claim for accuracy abs(i-result).le.max(epsabs,epsrel*abs(i)).
Parameters:
on entry f - real ( kind = 8 ) function subprogram defining the integrand function f(x). the actual name for f needs to be declared e x t e r n a l in the driver program.
a - real ( kind = 8 ) lower limit of integration
b - real ( kind = 8 ) upper limit of integration, b.gt.a if b.le.a, the routine will end with ier = 6.
alfa - real ( kind = 8 ) parameter in the weight function, alfa.gt.(-1) if alfa.le.(-1), the routine will end with ier = 6.
beta - real ( kind = 8 ) parameter in the weight function, beta.gt.(-1) if beta.le.(-1), the routine will end with ier = 6.
integr - integer ( kind = 4 ) indicates which weight function is to be used = 1 (x-a)**alfa*(b-x)**beta = 2 (x-a)**alfa*(b-x)**beta*log(x-a) = 3 (x-a)**alfa*(b-x)**beta*log(b-x) = 4 (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x) if integr.lt.1 or integr.gt.4, the routine will end with ier = 6.
epsabs - real ( kind = 8 ) absolute accuracy requested epsrel - real ( kind = 8 ) relative accuracy requested if epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), the routine will end with ier = 6.
limit - integer ( kind = 4 ) gives an upper bound on the number of subintervals in the partition of (a,b), limit.ge.2 if limit.lt.2, the routine will end with ier = 6.
on return result - real ( kind = 8 ) approximation to the integral
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
neval - integer ( kind = 4 ) number of integrand evaluations
ier - integer ( kind = 4 ) ier = 0 normal and reliable termination of the routine. it is assumed that the requested accuracy has been achieved. ier.gt.0 abnormal termination of the routine the estimates for the integral and error are less reliable. it is assumed that the requested accuracy has not been achieved. error messages = 1 maximum number of subdivisions allowed has been achieved. one can allow more subdivisions by increasing the value of limit. however, if this yields no improvement, it is advised to analyze the integrand in order to determine the integration difficulties which prevent the requested tolerance from being achieved. in case of a jump discontinuity or a local singularity of algebraico-logarithmic type at one or more interior points of the integration range, one should proceed by splitting up the interval at these points and calling the integrator on the subranges. = 2 the occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. = 3 extremely bad integrand behaviour occurs at some points of the integration interval. = 6 the input is invalid, because b.le.a or alfa.le.(-1) or beta.le.(-1), or integr.lt.1 or integr.gt.4, or (epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), or limit.lt.2. result, abserr, neval, rlist(1), elist(1), iord(1) and last are set to zero. alist(1) and blist(1) are set to a and b respectively.
alist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the left end points of the subintervals in the partition of the given integration range (a,b)
blist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the right end points of the subintervals in the partition of the given integration range (a,b)
rlist - real ( kind = 8 ) vector of dimension at least limit,the first last elements of which are the integral approximations on the subintervals
elist - real ( kind = 8 ) vector of dimension at least limit, the first last elements of which are the moduli of the absolute error estimates on the subintervals
iord - integer ( kind = 4 ) vector of dimension at least limit, the first k of which are pointers to the error estimates over the subintervals, so that elist(iord(1)), ..., elist(iord(k)) with k = last if last.le.(limit/2+2), and k = limit+1-last otherwise form a decreasing sequence
last - integer ( kind = 4 ) number of subintervals actually produced in the subdivision process
Local parameters:
alist - list of left end points of all subintervals considered up to now blist - list of right end points of all subintervals considered up to now rlist(i) - approximation to the integral over (alist(i),blist(i)) elist(i) - error estimate applying to rlist(i) maxerr - pointer to the interval with largest error estimate errmax - elist(maxerr) area - sum of the integrals over the subintervals errsum - sum of the errors over the subintervals errbnd - requested accuracy max(epsabs,epsrel* abs(result)) *****1 - variable for the left subinterval *****2 - variable for the right subinterval last - index for subdivision
machine dependent constants
epmach is the largest relative spacing. uflow is the smallest positive magnitude.
Definition at line 4960 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqc25c | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | c, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| integer ( kind = 4 ) | krul, | ||
| integer ( kind = 4 ) | neval | ||
| ) |
DQC25C returns integration rules for Cauchy Principal Value integrals.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose to compute i = integral of f*w over (a,b) with error estimate, where w(x) = 1/(x-c)
Parameters:
f - real ( kind = 8 ) function subprogram defining the integrand function f(x). the actual name for f needs to be declared e x t e r n a l in the driver program.
a - real ( kind = 8 ) left end point of the integration interval
b - real ( kind = 8 ) right end point of the integration interval, b.gt.a
c - real ( kind = 8 ) parameter in the weight function
result - real ( kind = 8 ) approximation to the integral result is computed by using a generalized clenshaw-curtis method if c lies within ten percent of the integration interval. in the other case the 15-point kronrod rule obtained by optimal addition of abscissae to the 7-point gauss rule, is applied.
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
krul - integer ( kind = 4 ) key which is decreased by 1 if the 15-point gauss-kronrod scheme has been used
neval - integer ( kind = 4 ) number of integrand evaluations
Local Parameters:
fval - value of the function f at the points cos(k*pi/24), k = 0, ..., 24 cheb12 - chebyshev series expansion coefficients, for the function f, of degree 12 cheb24 - chebyshev series expansion coefficients, for the function f, of degree 24 res12 - approximation to the integral corresponding to the use of cheb12 res24 - approximation to the integral corresponding to the use of cheb24 dqwgtc - external function subprogram defining the weight function hlgth - half-length of the interval centr - mid point of the interval
the vector x contains the values cos(k*pi/24), k = 1, ..., 11, to be used for the chebyshev series expansion of f
Definition at line 5426 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqc25f | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | omega, | ||
| integer ( kind = 4 ) | integr, | ||
| integer ( kind = 4 ) | nrmom, | ||
| integer ( kind = 4 ) | maxp1, | ||
| integer ( kind = 4 ) | ksave, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| integer ( kind = 4 ) | neval, | ||
| real ( kind = 8 ) | resabs, | ||
| real ( kind = 8 ) | resasc, | ||
| integer ( kind = 4 ) | momcom, | ||
| real ( kind = 8 ), dimension(maxp1,25) | chebmo | ||
| ) |
DQC25F returns integration rules for functions with a COS or SIN factor.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose to compute the integral i=integral of f(x) over (a,b) where w(x) = cos(omega*x) or w(x)=sin(omega*x) and to compute j = integral of abs(f) over (a,b). for small value of omega or small intervals (a,b) the 15-point gauss-kronro rule is used. otherwise a generalized clenshaw-curtis method is used.
Parameters:
on entry f - real ( kind = 8 ) function subprogram defining the integrand function f(x). the actual name for f needs to be declared e x t e r n a l in the calling program.
a - real ( kind = 8 ) lower limit of integration
b - real ( kind = 8 ) upper limit of integration
omega - real ( kind = 8 ) parameter in the weight function
integr - integer ( kind = 4 ) indicates which weight function is to be used integr = 1 w(x) = cos(omega*x) integr = 2 w(x) = sin(omega*x)
nrmom - integer ( kind = 4 ) the length of interval (a,b) is equal to the length of the original integration interval divided by 2**nrmom (we suppose that the routine is used in an adaptive integration process, otherwise set nrmom = 0). nrmom must be zero at the first call.
maxp1 - integer ( kind = 4 ) gives an upper bound on the number of chebyshev moments which can be stored, i.e. for the intervals of lengths abs(bb-aa)*2**(-l), l = 0,1,2, ..., maxp1-2.
ksave - integer ( kind = 4 ) key which is one when the moments for the current interval have been computed
on return result - real ( kind = 8 ) approximation to the integral i
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
neval - integer ( kind = 4 ) number of integrand evaluations
resabs - real ( kind = 8 ) approximation to the integral j
resasc - real ( kind = 8 ) approximation to the integral of abs(f-i/(b-a))
on entry and return momcom - integer ( kind = 4 ) for each interval length we need to compute the chebyshev moments. momcom counts the number of intervals for which these moments have already been computed. if nrmom.lt.momcom or ksave = 1, the chebyshev moments for the interval (a,b) have already been computed and stored, otherwise we compute them and we increase momcom.
chebmo - real ( kind = 8 ) array of dimension at least (maxp1,25) containing the modified chebyshev moments for the first momcom momcom interval lengths
Local Parameters:
the vector x contains the values cos(k*pi/24) k = 1, ...,11, to be used for the chebyshev expansion of f
centr - mid point of the integration interval hlgth - half-length of the integration interval fval - value of the function f at the points (b-a)*0.5*cos(k*pi/12) + (b+a)*0.5, k = 0, ..., 24 cheb12 - coefficients of the chebyshev series expansion of degree 12, for the function f, in the interval (a,b) cheb24 - coefficients of the chebyshev series expansion of degree 24, for the function f, in the interval (a,b) resc12 - approximation to the integral of cos(0.5*(b-a)*omega*x)*f(0.5*(b-a)*x+0.5*(b+a)) over (-1,+1), using the chebyshev series expansion of degree 12 resc24 - approximation to the same integral, using the chebyshev series expansion of degree 24 ress12 - the analogue of resc12 for the sine ress24 - the analogue of resc24 for the sine
machine dependent constant oflow is the largest positive magnitude.
Definition at line 5637 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqc25s | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | bl, | ||
| real ( kind = 8 ) | br, | ||
| real ( kind = 8 ) | alfa, | ||
| real ( kind = 8 ) | beta, | ||
| real ( kind = 8 ), dimension(25) | ri, | ||
| real ( kind = 8 ), dimension(25) | rj, | ||
| real ( kind = 8 ), dimension(25) | rg, | ||
| real ( kind = 8 ), dimension(25) | rh, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| real ( kind = 8 ) | resasc, | ||
| integer ( kind = 4 ) | integr, | ||
| integer ( kind = 4 ) | nev | ||
| ) |
DQC25S returns rules for algebraico-logarithmic end point singularities.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose to compute i = integral of f*w over (bl,br), with error estimate, where the weight function w has a singular behaviour of algebraico-logarithmic type at the points a and/or b. (bl,br) is a part of (a,b).
Parameters:
f - real ( kind = 8 ) function subprogram defining the integrand f(x). the actual name for f needs to be declared e x t e r n a l in the driver program.
a - real ( kind = 8 ) left end point of the original interval
b - real ( kind = 8 ) right end point of the original interval, b.gt.a
bl - real ( kind = 8 ) lower limit of integration, bl.ge.a
br - real ( kind = 8 ) upper limit of integration, br.le.b
alfa - real ( kind = 8 ) parameter in the weight function
beta - real ( kind = 8 ) parameter in the weight function
ri,rj,rg,rh - real ( kind = 8 ) modified chebyshev moments for the application of the generalized clenshaw-curtis method (computed in routine dqmomo)
result - real ( kind = 8 ) approximation to the integral result is computed by using a generalized clenshaw-curtis method if b1 = a or br = b. in all other cases the 15-point kronrod rule is applied, obtained by optimal addition of abscissae to the 7-point gauss rule.
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
resasc - real ( kind = 8 ) approximation to the integral of abs(f*w-i/(b-a))
integr - integer ( kind = 4 ) which determines the weight function = 1 w(x) = (x-a)**alfa*(b-x)**beta = 2 w(x) = (x-a)**alfa*(b-x)**beta*log(x-a) = 3 w(x) = (x-a)**alfa*(b-x)**beta*log(b-x) = 4 w(x) = (x-a)**alfa*(b-x)**beta*log(x-a)* log(b-x)
nev - integer ( kind = 4 ) number of integrand evaluations
Local Parameters:
the vector x contains the values cos(k*pi/24) k = 1, ..., 11, to be used for the computation of the chebyshev series expansion of f.
fval - value of the function f at the points (br-bl)*0.5*cos(k*pi/24)+(br+bl)*0.5 k = 0, ..., 24 cheb12 - coefficients of the chebyshev series expansion of degree 12, for the function f, in the interval (bl,br) cheb24 - coefficients of the chebyshev series expansion of degree 24, for the function f, in the interval (bl,br) res12 - approximation to the integral obtained from cheb12 res24 - approximation to the integral obtained from cheb24 dqwgts - external function subprogram defining the four possible weight functions hlgth - half-length of the interval (bl,br) centr - mid point of the interval (bl,br)
Definition at line 5966 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqcheb | ( | real ( kind = 8 ), dimension(11) | x, |
| real ( kind = 8 ), dimension(25) | fval, | ||
| real ( kind = 8 ), dimension(13) | cheb12, | ||
| real ( kind = 8 ), dimension(25) | cheb24 | ||
| ) |
DQCHEB computes the Chebyshev series expansion.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose this routine computes the chebyshev series expansion of degrees 12 and 24 of a function using a fast fourier transform method f(x) = sum(k=1,..,13) (cheb12(k)*t(k-1,x)), f(x) = sum(k=1,..,25) (cheb24(k)*t(k-1,x)), where t(k,x) is the chebyshev polynomial of degree k.
Parameters:
on entry x - real ( kind = 8 ) vector of dimension 11 containing the values cos(k*pi/24), k = 1, ..., 11
fval - real ( kind = 8 ) vector of dimension 25 containing the function values at the points (b+a+(b-a)*cos(k*pi/24))/2, k = 0, ...,24, where (a,b) is the approximation interval. fval(1) and fval(25) are divided by two (these values are destroyed at output).
on return cheb12 - real ( kind = 8 ) vector of dimension 13 containing the chebyshev coefficients for degree 12
cheb24 - real ( kind = 8 ) vector of dimension 25 containing the chebyshev coefficients for degree 24
Definition at line 6256 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqelg | ( | integer ( kind = 4 ) | n, |
| real ( kind = 8 ), dimension(52) | epstab, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| real ( kind = 8 ), dimension(3) | res3la, | ||
| integer ( kind = 4 ) | nres | ||
| ) |
DQELG carries out the Epsilon extrapolation algorithm.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose the routine determines the limit of a given sequence of approximations, by means of the epsilon algorithm of p.wynn. an estimate of the absolute error is also given. the condensed epsilon table is computed. only those elements needed for the computation of the next diagonal are preserved.
Parameters:
n - integer ( kind = 4 )
epstab(n) contains the new element in the
first column of the epsilon table.
epstab - real ( kind = 8 )
vector of dimension 52 containing the elements
of the two lower diagonals of the triangular
epsilon table. the elements are numbered
starting at the right-hand corner of the
triangle.
result - real ( kind = 8 )
resulting approximation to the integral
abserr - real ( kind = 8 )
estimate of the absolute error computed from
result and the 3 previous results
res3la - real ( kind = 8 )
vector of dimension 3 containing the last 3
results
nres - integer ( kind = 4 )
number of calls to the routine
(should be zero at first call)
Local Parameters:
e0 - the 4 elements on which the computation of a new e1 element in the epsilon table is based e2 e3 e0 e3 e1 new e2 newelm - number of elements to be computed in the new diagonal error - error = abs(e1-e0)+abs(e2-e1)+abs(new-e2) result - the element in the new diagonal with least value of error
machine dependent constants
epmach is the largest relative spacing. oflow is the largest positive magnitude. limexp is the maximum number of elements the epsilon table can contain. if this number is reached, the upper diagonal of the epsilon table is deleted.
Definition at line 6441 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqk15 | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| real ( kind = 8 ) | resabs, | ||
| real ( kind = 8 ) | resasc | ||
| ) |
DQK15 carries out a 15 point Gauss-Kronrod quadrature rule.
the abscissae and weights are given for the interval (-1,1). because of symmetry only the positive abscissae and their corresponding weights are given.
xgk - abscissae of the 15-point kronrod rule xgk(2), xgk(4), ... abscissae of the 7-point gauss rule xgk(1), xgk(3), ... abscissae which are optimally added to the 7-point gauss rule
wgk - weights of the 15-point kronrod rule
wg - weights of the 7-point gauss rule
gauss quadrature weights and kronron quadrature abscissae and weights as evaluated with 80 decimal digit arithmetic by l. w. fullerton, bell labs, nov. 1981.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose to compute i = integral of f over (a,b), with error estimate j = integral of abs(f) over (a,b) Parameters:
on entry
f - real ( kind = 8 )
function subprogram defining the integrand
function f(x). the actual name for f needs to be
declared e x t e r n a l in the calling program.
a - real ( kind = 8 )
lower limit of integration
b - real ( kind = 8 )
upper limit of integration
on return
result - real ( kind = 8 )
approximation to the integral i
result is computed by applying the 15-point
kronrod rule (resk) obtained by optimal addition
of abscissae to the7-point gauss rule(resg).
abserr - real ( kind = 8 )
estimate of the modulus of the absolute error,
which should not exceed abs(i-result)
resabs - real ( kind = 8 )
approximation to the integral j
resasc - real ( kind = 8 )
approximation to the integral of abs(f-i/(b-a))
over (a,b)
Local Parameters:
centr - mid point of the interval hlgth - half-length of the interval absc - abscissa fval* - function value resg - result of the 7-point gauss formula resk - result of the 15-point kronrod formula reskh - approximation to the mean value of f over (a,b), i.e. to i/(b-a)
machine dependent constants
epmach is the largest relative spacing. uflow is the smallest positive magnitude.
Definition at line 6641 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqk15i | ( | external | f, |
| real ( kind = 8 ) | boun, | ||
| integer ( kind = 4 ) | inf, | ||
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| real ( kind = 8 ) | resabs, | ||
| real ( kind = 8 ) | resasc | ||
| ) |
DQK15I applies a 15 point Gauss-Kronrod quadrature on an infinite interval.
the abscissae and weights are supplied for the interval (-1,1). because of symmetry only the positive abscissae and their corresponding weights are given.
xgk - abscissae of the 15-point kronrod rule xgk(2), xgk(4), ... abscissae of the 7-point gauss rule xgk(1), xgk(3), ... abscissae which are optimally added to the 7-point gauss rule
wgk - weights of the 15-point kronrod rule
wg - weights of the 7-point gauss rule, corresponding to the abscissae xgk(2), xgk(4), ... wg(1), wg(3), ... are set to zero.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose the original (infinite integration range is mapped onto the interval (0,1) and (a,b) is a part of (0,1). it is the purpose to compute i = integral of transformed integrand over (a,b), j = integral of abs(transformed integrand) over (a,b).
Parameters:
on entry
f - real ( kind = 8 )
fuction subprogram defining the integrand
function f(x). the actual name for f needs to be
declared e x t e r n a l in the calling program.
boun - real ( kind = 8 )
finite bound of original integration
range (set to zero if inf = +2)
inf - integer ( kind = 4 )
if inf = -1, the original interval is
(-infinity,bound),
if inf = +1, the original interval is
(bound,+infinity),
if inf = +2, the original interval is
(-infinity,+infinity) and
the integral is computed as the sum of two
integrals, one over (-infinity,0) and one over
(0,+infinity).
a - real ( kind = 8 )
lower limit for integration over subrange
of (0,1)
b - real ( kind = 8 )
upper limit for integration over subrange
of (0,1)
on return
result - real ( kind = 8 )
approximation to the integral i
result is computed by applying the 15-point
kronrod rule(resk) obtained by optimal addition
of abscissae to the 7-point gauss rule(resg).
abserr - real ( kind = 8 )
estimate of the modulus of the absolute error,
which should equal or exceed abs(i-result)
resabs - real ( kind = 8 )
approximation to the integral j
resasc - real ( kind = 8 )
approximation to the integral of
abs((transformed integrand)-i/(b-a)) over (a,b)
Local Parameters:
centr - mid point of the interval hlgth - half-length of the interval absc* - abscissa tabsc* - transformed abscissa fval* - function value resg - result of the 7-point gauss formula resk - result of the 15-point kronrod formula reskh - approximation to the mean value of the transformed integrand over (a,b), i.e. to i/(b-a)
machine dependent constants
epmach is the largest relative spacing. uflow is the smallest positive magnitude.
Definition at line 6832 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqk15w | ( | external | f, |
| external | w, | ||
| real ( kind = 8 ) | p1, | ||
| real ( kind = 8 ) | p2, | ||
| real ( kind = 8 ) | p3, | ||
| real ( kind = 8 ) | p4, | ||
| integer ( kind = 4 ) | kp, | ||
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| real ( kind = 8 ) | resabs, | ||
| real ( kind = 8 ) | resasc | ||
| ) |
DQK15W applies a 15 point Gauss-Kronrod rule for a weighted integrand.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose to compute i = integral of f*w over (a,b), with error estimate j = integral of abs(f*w) over (a,b)
Parameters:
on entry
f - real ( kind = 8 )
function subprogram defining the integrand
function f(x). the actual name for f needs to be
declared e x t e r n a l in the driver program.
w - real ( kind = 8 )
function subprogram defining the integrand
weight function w(x). the actual name for w
needs to be declared e x t e r n a l in the
calling program.
p1, p2, p3, p4 - real ( kind = 8 )
parameters in the weight function
kp - integer ( kind = 4 )
key for indicating the type of weight function
a - real ( kind = 8 )
lower limit of integration
b - real ( kind = 8 )
upper limit of integration
on return
result - real ( kind = 8 )
approximation to the integral i
result is computed by applying the 15-point
kronrod rule (resk) obtained by optimal addition
of abscissae to the 7-point gauss rule (resg).
abserr - real ( kind = 8 )
estimate of the modulus of the absolute error,
which should equal or exceed abs(i-result)
resabs - real ( kind = 8 )
approximation to the integral of abs(f)
resasc - real ( kind = 8 )
approximation to the integral of abs(f-i/(b-a))
Local Parameters:
the abscissae and weights are given for the interval (-1,1). because of symmetry only the positive abscissae and their corresponding weights are given.
xgk - abscissae of the 15-point gauss-kronrod rule xgk(2), xgk(4), ... abscissae of the 7-point gauss rule xgk(1), xgk(3), ... abscissae which are optimally added to the 7-point gauss rule
wgk - weights of the 15-point gauss-kronrod rule
wg - weights of the 7-point gauss rule
centr - mid point of the interval hlgth - half-length of the interval absc* - abscissa fval* - function value resg - result of the 7-point gauss formula resk - result of the 15-point kronrod formula reskh - approximation to the mean value of f*w over (a,b), i.e. to i/(b-a)
machine dependent constants
epmach is the largest relative spacing. uflow is the smallest positive magnitude.
Definition at line 7015 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqk21 | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| real ( kind = 8 ) | resabs, | ||
| real ( kind = 8 ) | resasc | ||
| ) |
DQK21 carries out a 21 point Gauss-Kronrod quadrature rule.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose to compute i = integral of f over (a,b), with error estimate j = integral of abs(f) over (a,b)
Parameters:
on entry
f - real ( kind = 8 )
function subprogram defining the integrand
function f(x). the actual name for f needs to be
declared e x t e r n a l in the driver program.
a - real ( kind = 8 )
lower limit of integration
b - real ( kind = 8 )
upper limit of integration
on return
result - real ( kind = 8 )
approximation to the integral i
result is computed by applying the 21-point
kronrod rule (resk) obtained by optimal addition
of abscissae to the 10-point gauss rule (resg).
abserr - real ( kind = 8 )
estimate of the modulus of the absolute error,
which should not exceed abs(i-result)
resabs - real ( kind = 8 )
approximation to the integral j
resasc - real ( kind = 8 )
approximation to the integral of abs(f-i/(b-a))
over (a,b)
Local Parameters:
the abscissae and weights are given for the interval (-1,1).
because of symmetry only the positive abscissae and their
corresponding weights are given.
xgk - abscissae of the 21-point kronrod rule
xgk(2), xgk(4), ... abscissae of the 10-point
gauss rule
xgk(1), xgk(3), ... abscissae which are optimally
added to the 10-point gauss rule
wgk - weights of the 21-point kronrod rule
wg - weights of the 10-point gauss rule
gauss quadrature weights and kronron quadrature abscissae and weights as evaluated with 80 decimal digit arithmetic by l. w. fullerton, bell labs, nov. 1981.
centr - mid point of the interval
hlgth - half-length of the interval
absc - abscissa
fval* - function value
resg - result of the 10-point gauss formula
resk - result of the 21-point kronrod formula
reskh - approximation to the mean value of f over (a,b),
i.e. to i/(b-a)
machine dependent constants
epmach is the largest relative spacing.
uflow is the smallest positive magnitude.
Definition at line 7190 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqk31 | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| real ( kind = 8 ) | resabs, | ||
| real ( kind = 8 ) | resasc | ||
| ) |
DQK31 carries out a 31 point Gauss-Kronrod quadrature rule.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose to compute i = integral of f over (a,b) with error estimate j = integral of abs(f) over (a,b)
Parameters:
on entry
f - real ( kind = 8 )
function subprogram defining the integrand
function f(x). the actual name for f needs to be
declared e x t e r n a l in the calling program.
a - real ( kind = 8 )
lower limit of integration
b - real ( kind = 8 )
upper limit of integration
on return
result - real ( kind = 8 )
approximation to the integral i
result is computed by applying the 31-point
gauss-kronrod rule (resk), obtained by optimal
addition of abscissae to the 15-point gauss
rule (resg).
abserr - double precison
estimate of the modulus of the modulus,
which should not exceed abs(i-result)
resabs - real ( kind = 8 )
approximation to the integral j
resasc - real ( kind = 8 )
approximation to the integral of abs(f-i/(b-a))
over (a,b)
Local Parameters:
the abscissae and weights are given for the interval (-1,1).
because of symmetry only the positive abscissae and their
corresponding weights are given.
xgk - abscissae of the 31-point kronrod rule
xgk(2), xgk(4), ... abscissae of the 15-point
gauss rule
xgk(1), xgk(3), ... abscissae which are optimally
added to the 15-point gauss rule
wgk - weights of the 31-point kronrod rule
wg - weights of the 15-point gauss rule
gauss quadrature weights and kronron quadrature abscissae and weights as evaluated with 80 decimal digit arithmetic by l. w. fullerton, bell labs, nov. 1981.
centr - mid point of the interval
hlgth - half-length of the interval
absc - abscissa
fval* - function value
resg - result of the 15-point gauss formula
resk - result of the 31-point kronrod formula
reskh - approximation to the mean value of f over (a,b),
i.e. to i/(b-a)
machine dependent constants
epmach is the largest relative spacing.
uflow is the smallest positive magnitude.
Definition at line 7372 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqk41 | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| real ( kind = 8 ) | resabs, | ||
| real ( kind = 8 ) | resasc | ||
| ) |
DQK41 carries out a 41 point Gauss-Kronrod quadrature rule.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose to compute i = integral of f over (a,b), with error estimate j = integral of abs(f) over (a,b)
Parameters:
on entry
f - real ( kind = 8 )
function subprogram defining the integrand
function f(x). the actual name for f needs to be
declared e x t e r n a l in the calling program.
a - real ( kind = 8 )
lower limit of integration
b - real ( kind = 8 )
upper limit of integration
on return
result - real ( kind = 8 )
approximation to the integral i
result is computed by applying the 41-point
gauss-kronrod rule (resk) obtained by optimal
addition of abscissae to the 20-point gauss
rule (resg).
abserr - real ( kind = 8 )
estimate of the modulus of the absolute error,
which should not exceed abs(i-result)
resabs - real ( kind = 8 )
approximation to the integral j
resasc - real ( kind = 8 )
approximation to the integal of abs(f-i/(b-a))
over (a,b)
Local Parameters:
the abscissae and weights are given for the interval (-1,1).
because of symmetry only the positive abscissae and their
corresponding weights are given.
xgk - abscissae of the 41-point gauss-kronrod rule
xgk(2), xgk(4), ... abscissae of the 20-point
gauss rule
xgk(1), xgk(3), ... abscissae which are optimally
added to the 20-point gauss rule
wgk - weights of the 41-point gauss-kronrod rule
wg - weights of the 20-point gauss rule
gauss quadrature weights and kronron quadrature abscissae and weights as evaluated with 80 decimal digit arithmetic by l. w. fullerton, bell labs, nov. 1981.
centr - mid point of the interval
hlgth - half-length of the interval
absc - abscissa
fval* - function value
resg - result of the 20-point gauss formula
resk - result of the 41-point kronrod formula
reskh - approximation to mean value of f over (a,b), i.e.
to i/(b-a)
machine dependent constants
epmach is the largest relative spacing.
uflow is the smallest positive magnitude.
Definition at line 7569 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqk51 | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| real ( kind = 8 ) | resabs, | ||
| real ( kind = 8 ) | resasc | ||
| ) |
DQK51 carries out a 51 point Gauss-Kronrod quadrature rule.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose to compute i = integral of f over (a,b) with error estimate j = integral of abs(f) over (a,b)
Parameters:
on entry
f - real ( kind = 8 )
function defining the integrand
function f(x). the actual name for f needs to be
declared e x t e r n a l in the calling program.
a - real ( kind = 8 )
lower limit of integration
b - real ( kind = 8 )
upper limit of integration
on return
result - real ( kind = 8 )
approximation to the integral i
result is computed by applying the 51-point
kronrod rule (resk) obtained by optimal addition
of abscissae to the 25-point gauss rule (resg).
abserr - real ( kind = 8 )
estimate of the modulus of the absolute error,
which should not exceed abs(i-result)
resabs - real ( kind = 8 )
approximation to the integral j
resasc - real ( kind = 8 )
approximation to the integral of abs(f-i/(b-a))
over (a,b)
Local Parameters:
the abscissae and weights are given for the interval (-1,1). because of symmetry only the positive abscissae and their corresponding weights are given.
xgk - abscissae of the 51-point kronrod rule xgk(2), xgk(4), ... abscissae of the 25-point gauss rule xgk(1), xgk(3), ... abscissae which are optimally added to the 25-point gauss rule
wgk - weights of the 51-point kronrod rule
wg - weights of the 25-point gauss rule
gauss quadrature weights and kronron quadrature abscissae and weights as evaluated with 80 decimal digit arithmetic by l. w. fullerton, bell labs, nov. 1981.
centr - mid point of the interval
hlgth - half-length of the interval
absc - abscissa
fval* - function value
resg - result of the 25-point gauss formula
resk - result of the 51-point kronrod formula
reskh - approximation to the mean value of f over (a,b),
i.e. to i/(b-a)
machine dependent constants
epmach is the largest relative spacing.
uflow is the smallest positive magnitude.
Definition at line 7775 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqk61 | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| real ( kind = 8 ) | resabs, | ||
| real ( kind = 8 ) | resasc | ||
| ) |
DQK61 carries out a 61 point Gauss-Kronrod quadrature rule.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose to compute i = integral of f over (a,b) with error estimate j = integral of abs ( f) over (a,b)
Parameters:
on entry f - real ( kind = 8 ) function subprogram defining the integrand function f(x). the actual name for f needs to be declared e x t e r n a l in the calling program.
a - real ( kind = 8 ) lower limit of integration
b - real ( kind = 8 ) upper limit of integration
on return result - real ( kind = 8 ) approximation to the integral i result is computed by applying the 61-point kronrod rule (resk) obtained by optimal addition of abscissae to the 30-point gauss rule (resg).
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs ( i-result)
resabs - real ( kind = 8 ) approximation to the integral j
resasc - real ( kind = 8 ) approximation to the integral of abs ( f-i/(b-a))
Local Parameters:
the abscissae and weights are given for the interval (-1,1). because of symmetry only the positive abscissae and their corresponding weights are given.
xgk - abscissae of the 61-point kronrod rule xgk(2), xgk(4) ... abscissae of the 30-point gauss rule xgk(1), xgk(3) ... optimally added abscissae to the 30-point gauss rule
wgk - weights of the 61-point kronrod rule
wg - weigths of the 30-point gauss rule
gauss quadrature weights and kronron quadrature abscissae and weights as evaluated with 80 decimal digit arithmetic by l. w. fullerton, bell labs, nov. 1981.
centr - mid point of the interval hlgth - half-length of the interval dabsc - abscissa fval* - function value resg - result of the 30-point gauss rule resk - result of the 61-point kronrod rule reskh - approximation to the mean value of f over (a,b), i.e. to i/(b-a)
machine dependent constants
epmach is the largest relative spacing. uflow is the smallest positive magnitude.
Definition at line 7994 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqmomo | ( | real ( kind = 8 ) | alfa, |
| real ( kind = 8 ) | beta, | ||
| real ( kind = 8 ), dimension(25) | ri, | ||
| real ( kind = 8 ), dimension(25) | rj, | ||
| real ( kind = 8 ), dimension(25) | rg, | ||
| real ( kind = 8 ), dimension(25) | rh, | ||
| integer ( kind = 4 ) | integr | ||
| ) |
DQMOMO computes modified Chebyshev moments.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose this routine computes modified chebsyshev moments. the k-th modified chebyshev moment is defined as the integral over (-1,1) of w(x)*t(k,x), where t(k,x) is the chebyshev polynomial of degree k.
Parameters:
alfa - real ( kind = 8 ) parameter in the weight function w(x), alfa.gt.(-1)
beta - real ( kind = 8 ) parameter in the weight function w(x), beta.gt.(-1)
ri - real ( kind = 8 ) vector of dimension 25 ri(k) is the integral over (-1,1) of (1+x)**alfa*t(k-1,x), k = 1, ..., 25.
rj - real ( kind = 8 ) vector of dimension 25 rj(k) is the integral over (-1,1) of (1-x)**beta*t(k-1,x), k = 1, ..., 25.
rg - real ( kind = 8 ) vector of dimension 25 rg(k) is the integral over (-1,1) of (1+x)**alfa*log((1+x)/2)*t(k-1,x), k = 1, ..., 25.
rh - real ( kind = 8 ) vector of dimension 25 rh(k) is the integral over (-1,1) of (1-x)**beta*log((1-x)/2)*t(k-1,x), k = 1, ..., 25.
integr - integer ( kind = 4 ) input parameter indicating the modified moments to be computed integr = 1 compute ri, rj = 2 compute ri, rj, rg = 3 compute ri, rj, rh = 4 compute ri, rj, rg, rh
Definition at line 8196 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqng | ( | external | f, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | epsabs, | ||
| real ( kind = 8 ) | epsrel, | ||
| real ( kind = 8 ) | result, | ||
| real ( kind = 8 ) | abserr, | ||
| integer ( kind = 4 ) | neval, | ||
| integer ( kind = 4 ) | ier | ||
| ) |
DQNG estimates an integral, using non-adaptive integration.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose the routine calculates an approximation result to a given definite integral i = integral of f over (a,b), hopefully satisfying following claim for accuracy abs(i-result).le.max(epsabs,epsrel*abs(i)).
Parameters:
f - real ( kind = 8 ) function subprogram defining the integrand function f(x). the actual name for f needs to be declared e x t e r n a l in the driver program.
a - real ( kind = 8 ) lower limit of integration
b - real ( kind = 8 ) upper limit of integration
epsabs - real ( kind = 8 ) absolute accuracy requested epsrel - real ( kind = 8 ) relative accuracy requested if epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), the routine will end with ier = 6.
on return result - real ( kind = 8 ) approximation to the integral i result is obtained by applying the 21-point gauss-kronrod rule (res21) obtained by optimal addition of abscissae to the 10-point gauss rule (res10), or by applying the 43-point rule (res43) obtained by optimal addition of abscissae to the 21-point gauss-kronrod rule, or by applying the 87-point rule (res87) obtained by optimal addition of abscissae to the 43-point rule.
abserr - real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
neval - integer ( kind = 4 ) number of integrand evaluations
ier - ier = 0 normal and reliable termination of the routine. it is assumed that the requested accuracy has been achieved. ier.gt.0 abnormal termination of the routine. it is assumed that the requested accuracy has not been achieved. error messages ier = 1 the maximum number of steps has been executed. the integral is probably too difficult to be calculated by dqng. = 6 the input is invalid, because epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28). result, abserr and neval are set to zero.
Local Parameters:
the data statements contain the abscissae and weights of the integration rules used.
x1 abscissae common to the 10-, 21-, 43- and 87- point rule x2 abscissae common to the 21-, 43- and 87-point rule x3 abscissae common to the 43- and 87-point rule x4 abscissae of the 87-point rule w10 weights of the 10-point formula w21a weights of the 21-point formula for abscissae x1 w21b weights of the 21-point formula for abscissae x2 w43a weights of the 43-point formula for abscissae x1, x3 w43b weights of the 43-point formula for abscissae x3 w87a weights of the 87-point formula for abscissae x1, x2, x3 w87b weights of the 87-point formula for abscissae x4
gauss-kronrod-patterson quadrature coefficients for use in quadpack routine qng. these coefficients were calculated with 101 decimal digit arithmetic by l. w. fullerton, bell labs, nov 1981.
centr - mid point of the integration interval
hlgth - half-length of the integration interval
fcentr - function value at mid point
absc - abscissa
fval - function value
savfun - array of function values which have already been
computed
res10 - 10-point gauss result
res21 - 21-point kronrod result
res43 - 43-point result
res87 - 87-point result
resabs - approximation to the integral of abs(f)
resasc - approximation to the integral of abs(f-i/(b-a))
machine dependent constants
epmach is the largest relative spacing.
uflow is the smallest positive magnitude.
Definition at line 8394 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::dqpsrt | ( | integer ( kind = 4 ) | limit, |
| integer ( kind = 4 ) | last, | ||
| integer ( kind = 4 ) | maxerr, | ||
| real ( kind = 8 ) | ermax, | ||
| real ( kind = 8 ), dimension(last) | elist, | ||
| integer ( kind = 4 ), dimension(last) | iord, | ||
| integer ( kind = 4 ) | nrmax | ||
| ) |
DQPSRT maintains the order of a list of local error estimates.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
***purpose this routine maintains the descending ordering in the list of the local error estimated resulting from the interval subdivision process. at each call two error estimates are inserted using the sequential search method, top-down for the largest error estimate and bottom-up for the smallest error estimate.
Parameters:
limit - integer ( kind = 4 )
maximum number of error estimates the list
can contain
last - integer ( kind = 4 )
number of error estimates currently in the list
maxerr - integer ( kind = 4 )
maxerr points to the nrmax-th largest error
estimate currently in the list
ermax - real ( kind = 8 )
nrmax-th largest error estimate
ermax = elist(maxerr)
elist - real ( kind = 8 )
vector of dimension last containing
the error estimates
iord - integer ( kind = 4 )
vector of dimension last, the first k elements
of which contain pointers to the error
estimates, such that
elist(iord(1)),..., elist(iord(k))
form a decreasing sequence, with
k = last if last.le.(limit/2+2), and
k = limit+1-last otherwise
nrmax - integer ( kind = 4 )
maxerr = iord(nrmax)
Definition at line 8723 of file 02_quadpack_double.f90.
| real ( kind = 8 ) function eftcamb_quadpack::dqwgtc | ( | real ( kind = 8 ) | x, |
| real ( kind = 8 ) | c, | ||
| real ( kind = 8 ) | p2, | ||
| real ( kind = 8 ) | p3, | ||
| real ( kind = 8 ) | p4, | ||
| integer ( kind = 4 ) | kp | ||
| ) |
DQWGTC defines the weight function used by DQC25C.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
Definition at line 104 of file 02_quadpack_double.f90.
| real ( kind = 8 ) function eftcamb_quadpack::dqwgtf | ( | real ( kind = 8 ) | x, |
| real ( kind = 8 ) | omega, | ||
| real ( kind = 8 ) | p2, | ||
| real ( kind = 8 ) | p3, | ||
| real ( kind = 8 ) | p4, | ||
| integer ( kind = 4 ) | integr | ||
| ) |
DQWGTF defines the weight functions used by DQC25F.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
Definition at line 129 of file 02_quadpack_double.f90.
| real ( kind = 8 ) function eftcamb_quadpack::dqwgts | ( | real ( kind = 8 ) | x, |
| real ( kind = 8 ) | a, | ||
| real ( kind = 8 ) | b, | ||
| real ( kind = 8 ) | alfa, | ||
| real ( kind = 8 ) | beta, | ||
| integer ( kind = 4 ) | integr | ||
| ) |
DQWGTS defines the weight functions used by DQC25S.
Modified:
11 September 2015
Author:
Robert Piessens, Elise de Doncker
Definition at line 71 of file 02_quadpack_double.f90.
| subroutine eftcamb_quadpack::xerror | ( | character ( len = * ) | xmess, |
| integer ( kind = 4 ) | nmess, | ||
| integer ( kind = 4 ) | nerr, | ||
| integer ( kind = 4 ) | level | ||
| ) |
XERROR replaces the SLATEC XERROR routine.
Modified:
12 September 2015
Definition at line 41 of file 02_quadpack_double.f90.
1.8.11