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 semiinfinite 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 semiinfinite 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 algebraicologarithmic 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 algebraicologarithmic 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 algebraicologarithmic 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 GaussKronrod quadrature rule. More...  
subroutine  dqk15i (f, boun, inf, a, b, result, abserr, resabs, resasc) 
DQK15I applies a 15 point GaussKronrod 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 GaussKronrod rule for a weighted integrand. More...  
subroutine  dqk21 (f, a, b, result, abserr, resabs, resasc) 
DQK21 carries out a 21 point GaussKronrod quadrature rule. More...  
subroutine  dqk31 (f, a, b, result, abserr, resabs, resasc) 
DQK31 carries out a 31 point GaussKronrod quadrature rule. More...  
subroutine  dqk41 (f, a, b, result, abserr, resabs, resasc) 
DQK41 carries out a 41 point GaussKronrod quadrature rule. More...  
subroutine  dqk51 (f, a, b, result, abserr, resabs, resasc) 
DQK51 carries out a 51 point GaussKronrod quadrature rule. More...  
subroutine  dqk61 (f, a, b, result, abserr, resabs, resasc) 
DQK61 carries out a 61 point GaussKronrod 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 nonadaptive 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: 9780898711721, 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:N1). 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 Kth 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(iresult)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.5d28), the routine will end with ier = 6. key  integer ( kind = 4 ) key for choice of local integration rule a gausskronrod 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(iresult)
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 specialpurpose 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.5d28)) 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+1last 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(ireslt).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.5d28), the routine will end with ier = 6.
key  integer ( kind = 4 ) key for choice of local integration rule a gausskronrod 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(iresult)
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 specialpurpose 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.5d28), 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+1last 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 semiinfinite 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(iresult).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 doublyinfinite)
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.5d28), 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(iresult)
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+1last 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 semiinfinite 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(iresult).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 doublyinfinite) 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.5d28), 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(iresult)
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+1last 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 (truevalue)
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 usersupplied 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 (npts22) 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.5d28), 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(iresult)
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 specialpurpose 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 underestimated. = 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.5d28)) 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 = (leniwnpts2)/2).
dimensioning parameters leniw  integer ( kind = 4 ) dimensioning parameter for iwork leniw determines limit = (leniwnpts2)/2, which is the maximum number of subintervals in the partition of the given integration interval (a,b), leniw.ge.(3*npts22). if leniw.lt.(3*npts22), the routine will end with ier = 6.
lenw  integer ( kind = 4 ) dimensioning parameter for work lenw must be at least leniw*2npts2. if lenw.lt.leniw*2npts2, 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+1last 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 userprovided break point or integration limit, then (aa,bb) has level l if abs(bbaa) = abs(p2p1)*2**(l), iwork(limit*2+1), ..., iwork(limit*2+npts2) have no significance for the user, note that limit = (leniwnpts2)/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 = (leniwnpts2)/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(iresult).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 usersupplied 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 (npts22) 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.5d28), 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(iresult)
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 specialpurpose 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 underestimated. = 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.5d28)) 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 userprovided break point or integration limit, then (aa,bb) has level l if abs(bbaa) = abs(p2p1)*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, ..., npts22, 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+1last 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 (truevalue)
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(iresult).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.5d28), 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(iresult)
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 specialpurpose 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 underestimated. = 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.5d28) 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+1last otherwise
work  real ( kind = 8 ) vector of dimension at least lenw on return work(1), ..., work(last) contain the left endpoints of the subintervals in the partition of (a,b), work(limit+1), ..., work(limit+last) contain the right endpoints, 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(iresult).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.5d28), 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(iresult)
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 specialpurpose 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 underestimated. = 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.5d28). 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+1last 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/((xc), c.ne.a, c.ne.b), hopefully satisfying following claim for accuracy abs(iresult).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.5d28), 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(iresult)
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.5d28)) 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+1last 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/(xc), (c.ne.a, c.ne.b), hopefully satisfying following claim for accuracy abs(iresult).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.5d28), 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(iresult)
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.5d28)) 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+1last 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(iresult).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(iresult)
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+(k1)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 (=(leniwlimlst) /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, (leniwlimlst)/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(ba)*2**(l), l = 0,1, ..., maxp12, 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(iresult).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(ba)*2**(l), l=0,1, ..., maxp12, 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(iresult)
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+(k1)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+(k1)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(iresult).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.5d28), 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(iresult)
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(ba)*2**(l), l=0,1, ..., maxp12, 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+1last 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(ba)*2**(1l).
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(iresult).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.5d28), 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 clenshawcurtis integration of degree 24) have been computed for intervals of lenghts (abs(ba))*2**(l), l=0,1,2,...momcom1. if icall.gt.1 this means that dqawoe has been called twice or more on intervals of the same length abs(ba). the chebyshev moments already computed are then reused 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(ba)*2**(l), l=0,1, ..., maxp12, 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(iresult)
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+1last 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(ba)*2**(1l)
on entry and return momcom  integer ( kind = 4 ) indicating that the chebyshev moments have been computed for intervals of lengths (abs(ba))*2**(l), l=0,1,2, ..., momcom1, 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 algebraicologarithmic 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(iresult).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 (xa)**alfa*(bx)**beta = 2 (xa)**alfa*(bx)**beta*log(xa) = 3 (xa)**alfa*(bx)**beta*log(bx) = 4 (xa)**alfa*(bx)**beta*log(xa)*log(bx) 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.5d28), 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(iresult)
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 algebraicologarithmic 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.5d28)) 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+1last 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 algebraicologarithmic 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(iresult).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 (xa)**alfa*(bx)**beta = 2 (xa)**alfa*(bx)**beta*log(xa) = 3 (xa)**alfa*(bx)**beta*log(bx) = 4 (xa)**alfa*(bx)**beta*log(xa)*log(bx) 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.5d28), 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(iresult)
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 algebraicologarithmic 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.5d28), 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+1last 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/(xc)
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 clenshawcurtis method if c lies within ten percent of the integration interval. in the other case the 15point kronrod rule obtained by optimal addition of abscissae to the 7point gauss rule, is applied.
abserr  real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(iresult)
krul  integer ( kind = 4 ) key which is decreased by 1 if the 15point gausskronrod 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  halflength 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 15point gausskronro rule is used. otherwise a generalized clenshawcurtis 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(bbaa)*2**(l), l = 0,1,2, ..., maxp12.
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(iresult)
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(fi/(ba))
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  halflength of the integration interval fval  value of the function f at the points (ba)*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*(ba)*omega*x)*f(0.5*(ba)*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 algebraicologarithmic 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 algebraicologarithmic 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 clenshawcurtis method (computed in routine dqmomo)
result  real ( kind = 8 ) approximation to the integral result is computed by using a generalized clenshawcurtis method if b1 = a or br = b. in all other cases the 15point kronrod rule is applied, obtained by optimal addition of abscissae to the 7point gauss rule.
abserr  real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(iresult)
resasc  real ( kind = 8 ) approximation to the integral of abs(f*wi/(ba))
integr  integer ( kind = 4 ) which determines the weight function = 1 w(x) = (xa)**alfa*(bx)**beta = 2 w(x) = (xa)**alfa*(bx)**beta*log(xa) = 3 w(x) = (xa)**alfa*(bx)**beta*log(bx) = 4 w(x) = (xa)**alfa*(bx)**beta*log(xa)* log(bx)
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 (brbl)*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  halflength 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(k1,x)), f(x) = sum(k=1,..,25) (cheb24(k)*t(k1,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+(ba)*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 righthand 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(e1e0)+abs(e2e1)+abs(newe2) 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 GaussKronrod 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 15point kronrod rule xgk(2), xgk(4), ... abscissae of the 7point gauss rule xgk(1), xgk(3), ... abscissae which are optimally added to the 7point gauss rule
wgk  weights of the 15point kronrod rule
wg  weights of the 7point 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 15point kronrod rule (resk) obtained by optimal addition of abscissae to the7point gauss rule(resg). abserr  real ( kind = 8 ) estimate of the modulus of the absolute error, which should not exceed abs(iresult) resabs  real ( kind = 8 ) approximation to the integral j resasc  real ( kind = 8 ) approximation to the integral of abs(fi/(ba)) over (a,b)
Local Parameters:
centr  mid point of the interval hlgth  halflength of the interval absc  abscissa fval*  function value resg  result of the 7point gauss formula resk  result of the 15point kronrod formula reskh  approximation to the mean value of f over (a,b), i.e. to i/(ba)
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 GaussKronrod 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 15point kronrod rule xgk(2), xgk(4), ... abscissae of the 7point gauss rule xgk(1), xgk(3), ... abscissae which are optimally added to the 7point gauss rule
wgk  weights of the 15point kronrod rule
wg  weights of the 7point 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 15point kronrod rule(resk) obtained by optimal addition of abscissae to the 7point gauss rule(resg). abserr  real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(iresult) resabs  real ( kind = 8 ) approximation to the integral j resasc  real ( kind = 8 ) approximation to the integral of abs((transformed integrand)i/(ba)) over (a,b)
Local Parameters:
centr  mid point of the interval hlgth  halflength of the interval absc*  abscissa tabsc*  transformed abscissa fval*  function value resg  result of the 7point gauss formula resk  result of the 15point kronrod formula reskh  approximation to the mean value of the transformed integrand over (a,b), i.e. to i/(ba)
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 GaussKronrod 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 15point kronrod rule (resk) obtained by optimal addition of abscissae to the 7point gauss rule (resg). abserr  real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(iresult) resabs  real ( kind = 8 ) approximation to the integral of abs(f) resasc  real ( kind = 8 ) approximation to the integral of abs(fi/(ba))
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 15point gausskronrod rule xgk(2), xgk(4), ... abscissae of the 7point gauss rule xgk(1), xgk(3), ... abscissae which are optimally added to the 7point gauss rule
wgk  weights of the 15point gausskronrod rule
wg  weights of the 7point gauss rule
centr  mid point of the interval hlgth  halflength of the interval absc*  abscissa fval*  function value resg  result of the 7point gauss formula resk  result of the 15point kronrod formula reskh  approximation to the mean value of f*w over (a,b), i.e. to i/(ba)
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 GaussKronrod 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 21point kronrod rule (resk) obtained by optimal addition of abscissae to the 10point gauss rule (resg). abserr  real ( kind = 8 ) estimate of the modulus of the absolute error, which should not exceed abs(iresult) resabs  real ( kind = 8 ) approximation to the integral j resasc  real ( kind = 8 ) approximation to the integral of abs(fi/(ba)) 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 21point kronrod rule xgk(2), xgk(4), ... abscissae of the 10point gauss rule xgk(1), xgk(3), ... abscissae which are optimally added to the 10point gauss rule wgk  weights of the 21point kronrod rule wg  weights of the 10point 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  halflength of the interval absc  abscissa fval*  function value resg  result of the 10point gauss formula resk  result of the 21point kronrod formula reskh  approximation to the mean value of f over (a,b), i.e. to i/(ba) 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 GaussKronrod 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 31point gausskronrod rule (resk), obtained by optimal addition of abscissae to the 15point gauss rule (resg). abserr  double precison estimate of the modulus of the modulus, which should not exceed abs(iresult) resabs  real ( kind = 8 ) approximation to the integral j resasc  real ( kind = 8 ) approximation to the integral of abs(fi/(ba)) 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 31point kronrod rule xgk(2), xgk(4), ... abscissae of the 15point gauss rule xgk(1), xgk(3), ... abscissae which are optimally added to the 15point gauss rule wgk  weights of the 31point kronrod rule wg  weights of the 15point 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  halflength of the interval absc  abscissa fval*  function value resg  result of the 15point gauss formula resk  result of the 31point kronrod formula reskh  approximation to the mean value of f over (a,b), i.e. to i/(ba) 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 GaussKronrod 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 41point gausskronrod rule (resk) obtained by optimal addition of abscissae to the 20point gauss rule (resg). abserr  real ( kind = 8 ) estimate of the modulus of the absolute error, which should not exceed abs(iresult) resabs  real ( kind = 8 ) approximation to the integral j resasc  real ( kind = 8 ) approximation to the integal of abs(fi/(ba)) 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 41point gausskronrod rule xgk(2), xgk(4), ... abscissae of the 20point gauss rule xgk(1), xgk(3), ... abscissae which are optimally added to the 20point gauss rule wgk  weights of the 41point gausskronrod rule wg  weights of the 20point 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  halflength of the interval absc  abscissa fval*  function value resg  result of the 20point gauss formula resk  result of the 41point kronrod formula reskh  approximation to mean value of f over (a,b), i.e. to i/(ba) 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 GaussKronrod 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 51point kronrod rule (resk) obtained by optimal addition of abscissae to the 25point gauss rule (resg). abserr  real ( kind = 8 ) estimate of the modulus of the absolute error, which should not exceed abs(iresult) resabs  real ( kind = 8 ) approximation to the integral j resasc  real ( kind = 8 ) approximation to the integral of abs(fi/(ba)) 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 51point kronrod rule xgk(2), xgk(4), ... abscissae of the 25point gauss rule xgk(1), xgk(3), ... abscissae which are optimally added to the 25point gauss rule
wgk  weights of the 51point kronrod rule
wg  weights of the 25point 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  halflength of the interval absc  abscissa fval*  function value resg  result of the 25point gauss formula resk  result of the 51point kronrod formula reskh  approximation to the mean value of f over (a,b), i.e. to i/(ba) 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 GaussKronrod 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 61point kronrod rule (resk) obtained by optimal addition of abscissae to the 30point gauss rule (resg).
abserr  real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs ( iresult)
resabs  real ( kind = 8 ) approximation to the integral j
resasc  real ( kind = 8 ) approximation to the integral of abs ( fi/(ba))
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 61point kronrod rule xgk(2), xgk(4) ... abscissae of the 30point gauss rule xgk(1), xgk(3) ... optimally added abscissae to the 30point gauss rule
wgk  weights of the 61point kronrod rule
wg  weigths of the 30point 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  halflength of the interval dabsc  abscissa fval*  function value resg  result of the 30point gauss rule resk  result of the 61point kronrod rule reskh  approximation to the mean value of f over (a,b), i.e. to i/(ba)
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 kth 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(k1,x), k = 1, ..., 25.
rj  real ( kind = 8 ) vector of dimension 25 rj(k) is the integral over (1,1) of (1x)**beta*t(k1,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(k1,x), k = 1, ..., 25.
rh  real ( kind = 8 ) vector of dimension 25 rh(k) is the integral over (1,1) of (1x)**beta*log((1x)/2)*t(k1,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 nonadaptive 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(iresult).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.5d28), the routine will end with ier = 6.
on return result  real ( kind = 8 ) approximation to the integral i result is obtained by applying the 21point gausskronrod rule (res21) obtained by optimal addition of abscissae to the 10point gauss rule (res10), or by applying the 43point rule (res43) obtained by optimal addition of abscissae to the 21point gausskronrod rule, or by applying the 87point rule (res87) obtained by optimal addition of abscissae to the 43point rule.
abserr  real ( kind = 8 ) estimate of the modulus of the absolute error, which should equal or exceed abs(iresult)
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.5d28). 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 87point rule x3 abscissae common to the 43 and 87point rule x4 abscissae of the 87point rule w10 weights of the 10point formula w21a weights of the 21point formula for abscissae x1 w21b weights of the 21point formula for abscissae x2 w43a weights of the 43point formula for abscissae x1, x3 w43b weights of the 43point formula for abscissae x3 w87a weights of the 87point formula for abscissae x1, x2, x3 w87b weights of the 87point formula for abscissae x4
gausskronrodpatterson 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  halflength 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  10point gauss result res21  21point kronrod result res43  43point result res87  87point result resabs  approximation to the integral of abs(f) resasc  approximation to the integral of abs(fi/(ba)) 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, topdown for the largest error estimate and bottomup 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 nrmaxth largest error estimate currently in the list ermax  real ( kind = 8 ) nrmaxth 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+1last 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.