      PROGRAM tglf6
************************************************************************
*     (Test Gauss-Laguerre Quadrature with Function values)
*
*     This is a driver program for subroutine glqf() which returns the
*     nodes and weights necessary for the ordinary Laguerre quadrature.
*
*     The test function is sin(sigma*x) with weight
*
*             x^alpha * exp(-x),      (alpha > -1)
*
*     integrated from 0 to infinity.
*
*     Normal input is on stdin, and normal output on stdout.
*
*     However, there is also an auxiliary file created called simply
*     "output.dat"; it contains data suitable for comparing against
*     similar files from other programs in an integral-independent
*     way.
*
*     Lines in this file are either comment lines (beginning with a
*     sharp sign), or blank or empty lines, or whitespace-separated
*     integral data lines of the form
*
*     np p(1) p(2) ... p(np) neval result opt-relerr opt-abserr
*
*     Here, np is the number of parameters that follow in the p(*)
*     values, neval is the number of function evaluations (0 if
*     unknown or unavailable), result is the floating-point computed
*     value of the integral, opt-relerr is an optional relative error
*     estimate, and opt-abserr is an optional absolute error estimate.
*     If opt-relerr is omitted, then opt-abserr cannot be specified.
*     A zero value for either implies that the value is unknown.
*
*     Normally, comment lines will document what integral is
*     evaluated, and relate any parameters in the integrand to the
*     array elements p(*).  It is acceptable for np to be 0: no p(*)
*     elements are then provided.
*
*     The availability of data files in this standard format makes it
*     relatively easy to compare results of different integrators.  In
*     particular, high-precision values can be computed in symbolic
*     algebra systems, such as Maple, Mathematica, Axiom, Reduce,
*     muPAD, ..., and used to evaluate the accuracy of results from
*     other integrators.
*
*     (13-May-2000)
************************************************************************
*
*     External functions
*
      EXTERNAL            deps,        derbit,      dfloat,      dgamma
      EXTERNAL            dnan,        dvsum
*
      DOUBLE PRECISION    dabs,        datan,       deps,        derbit
      DOUBLE PRECISION    dfloat,      dgamma,      dsin,        dnan
      DOUBLE PRECISION    dvsum
*
*     Statement functions
*
      DOUBLE PRECISION    f
*
*     Parameter variables
*
      CHARACTER*(*)       METHOD
      PARAMETER           (METHOD = 'glqf()')
*
      DOUBLE PRECISION    HALF
      PARAMETER           (HALF = 0.5d+00)
*
      DOUBLE PRECISION    ONE
      PARAMETER           (ONE = 1.0d+00)
*
      DOUBLE PRECISION    ZERO
      PARAMETER           (ZERO = 0.0d+00)
*
      INCLUDE 'maxpts.inc'
      INCLUDE 'stdio.inc'
*
*     Local variables
*
      DOUBLE PRECISION    alfmax,      alfmin,      alpha,       exact
      DOUBLE PRECISION    onepa,       oneps2,      relavg,      relerr
      DOUBLE PRECISION    relmax,      relmin,      relulp,      result
      DOUBLE PRECISION    sigarg,      sigma,       sigmax,      sigmin
      DOUBLE PRECISION    theta,       ulp,         v(MAXPTS)
      DOUBLE PRECISION    w(MAXPTS),   wxm1(MAXPTS)
      DOUBLE PRECISION    x(MAXPTS),   xarg,        y(MAXPTS)
      DOUBLE PRECISION    z(MAXPTS),   zinit(1)
*
      INTEGER             i,           ierr,        neval,       nqmax
      INTEGER             nqmin,       nquad,       nruns
*
      f(xarg,sigarg) = dsin(sigarg*xarg)
*     fprime(xarg,sigarg) = sigarg * dcos(sigarg*xarg)
*
*     Sun floating-point error trapping:
*
*     INTEGER             ieee_handler
*     EXTERNAL            trapit
*     ierr = ieee_handler ( 'set', 'invalid', trapit )
*     if (ierr .ne. 0) print *,'failed to set trap for invalid'
*
*     Initialize all local floating-point variables to NaN, or at
*     least an approximation thereto.
*
      zinit(1) = dnan()
      CALL dcopy (MAXPTS, zinit, 0, w, 1)
      CALL dcopy (MAXPTS, zinit, 0, wxm1, 1)
      CALL dcopy (MAXPTS, zinit, 0, x, 1)
      CALL dcopy (MAXPTS, zinit, 0, y, 1)
      CALL dcopy (MAXPTS, zinit, 0, z, 1)
      alfmax = zinit(1)
      alfmin = zinit(1)
      alpha = zinit(1)
      exact = zinit(1)
      onepa = zinit(1)
      oneps2 = zinit(1)
      relavg = zinit(1)
      relerr = zinit(1)
      relmax = zinit(1)
      relmin = zinit(1)
      relulp = zinit(1)
      result = zinit(1)
      sigarg = zinit(1)
      sigma = zinit(1)
      sigmax = zinit(1)
      sigmin = zinit(1)
      theta = zinit(1)
      ulp = zinit(1)
      xarg = zinit(1)
*
      OPEN (UNIT=stddat, FILE='output.dat', STATUS='unknown',
     X    FORM='formatted')
*
      WRITE (stddat,'(A)') '### Numerical integration with ' // METHOD
      WRITE (stddat,'(A)') '###'
      WRITE (stddat,'(A)') '### int(sin(sigma*x) * x^alpha * ' //
     X     'exp(-x), x = 0..infinity)'
      WRITE (stddat,'(A)') '###'
      WRITE (stddat,'(A)')
     X    '### Line format: 2 sigma alpha neval result relerr abserr'
      WRITE (stddat,'(A)') '###'
*
      ulp = deps(ONE)
      sigmin = 1.0d+75
      sigmax = ZERO
      nqmax = 0
      nqmin = 0
      nruns = 0
      relavg = ZERO
      relmax = ZERO
      relmin = 1.0d+75
*
*     Write a standard output header identifying the integral, and the
*     host precision
*
      CALL prthdr(stdout,
     X     'int(sin(sigma*x) * x^alpha * exp(-x), x = 0..infinity)')
      WRITE (stdout, 10000)
*
*     Loop reading input data:
*
  100 READ (stdin, *, END=300, ERR=300) nquad, sigma, alpha
      CALL glqf(x, w, wxm1, y, z, alpha, nquad, ierr)
      IF (ierr .NE. 0) THEN
          WRITE (stderr, 60000) METHOD, ierr
          GO TO 100
      END IF
*
*     Evaluate test integral as quadrature sum:
*
      DO 200 i = 1, nquad
          v(i) = w(i) * f(x(i),sigma)
  200 CONTINUE
      neval = nquad
      result = dvsum(v, nquad)
*
*     Analytic calculation of test integral.  This is a known integral
*     (Gradshteyn and Ryzhik, 4th edition, 3.944.5, p. 490):
*
*     S = \int_{0}^{\infty} dx\, x^\alpha e^{-x} \sin(\sigma x)
*       = \Gamma(1 + \alpha) (1 + \sigma^2)^{-(1 + \alpha)/2}
*         \sin((1 + \alpha)\theta)
*
*     where
*
*     \theta = arctan(\sigma)
*
      onepa = ONE + alpha
      oneps2 = ONE + sigma**2
      theta = datan(dabs(sigma))
      exact = dgamma(onepa) * dsin(onepa * theta) /
     X    oneps2**(HALF * onepa)
      IF (sigma .LT. ZERO) exact = -exact
      relerr = (exact - result)/exact
      relulp = relerr/ulp
*
*     For low quadrature order, the relative error in ulps can be
*     quite large, so we switch from Fw.d to Ew.d output format.
*
      IF (dabs(relulp) .LT. 10000.0d+00) THEN
          WRITE (stdout, 20000) nquad, sigma, alpha, result, exact,
     X         relerr, relulp, derbit(relerr,ulp)
      ELSE
          WRITE (stdout, 30000) nquad, sigma, alpha, result, exact,
     X         relerr, relulp, derbit(relerr,ulp)
      END IF
      WRITE (stddat, 50000) 2, sigma, alpha, neval, result, relerr,
     X        (exact - result)
      if (relerr .eq. relerr) THEN
*
*         Track the average, minimum, and maximum relative error as
*         long as it is not a NaN.
*
          relavg = relavg + dabs(relerr)
          nruns = nruns + 1
          IF (dabs(relerr) .GT. dabs(relmax)) THEN
              relmax = relerr
              nqmax = nquad
              sigmax = sigma
              alfmax = alpha
          END IF
          IF (dabs(relerr) .LT. dabs(relmin)) THEN
              relmin = relerr
              nqmin = nquad
              sigmin = sigma
              alfmin = alpha
          END IF
      END IF
      GO TO 100
*
  300 IF (nruns .GT. 0) THEN
          relavg = relavg / dfloat(nruns)
          WRITE (stdout, '()')
          WRITE (stdout, 70000) 'Maximum', relmax, relmax/ulp,
     X        derbit(relmax,ulp), 'at nquad =', nqmax,
     X        'sigma =', sigmax, 'alpha =', alfmax
          WRITE (stdout, 70000) 'Minimum', relmin, relmin/ulp,
     X        derbit(relmin,ulp), 'at nquad =', nqmin,
     X        'sigma =', sigmin, 'alpha =', alfmin
          WRITE (stdout, 70000) 'Average', relavg, relavg/ulp,
     X        derbit(relavg,ulp)
      END IF
      WRITE (stdout, 40000)
      CLOSE (UNIT=stddat)
*
10000 FORMAT (/,2X, 1X, 'nquad', 3X, 2X, 'sigma', 3X, 2X, 'alpha', 3X,
     X    2X, 'Quadrature Result', 5X,
     X    3X, 'Exact Integral', 4X,
     X    1X, 'Rel. Error', 1X, 'RelE (ULPs)', 2X, 'Err (bits)')
20000 FORMAT (2X, I6, 1X, F9.4, 1X, F9.4, 1P, 1X, E23.15, 1X, E23.15,
     X    1X, E10.2, 2X, 0P, 1X, F9.2, 1X, F9.2)
30000 FORMAT (2X, I6, 1X, F9.4, 1X, F9.4, 1P, 1X, E23.15, 1X, E23.15,
     X    1X, E10.2, 2X, E10.2, 0P, 1X, F9.2)
40000 FORMAT (/, 2X, 'Done')
50000 FORMAT (I1, 1X, 1X, 1P, 1X, E26.18, 1X, E26.18, 1X, I10, 1X,
     X    E27.18, 1X, E9.2, 1X, E9.2)
60000 FORMAT (/, 2X, 'ERROR: ', A, ' returns ierr = ', I10)
70000 FORMAT (2X, A, ' relative error = ', 1P, E10.2, 1X,
     X    1P, E10.2, ' ulps', 1X, 0P, F10.2, ' bits', 1X, A, I5,
     X    1X, A, F10.4, 1X, A, F10.4)
*
      END
