src/numerical/slatec/fortran/dqagi.f

   1 *DECK DQAGI
   2       SUBROUTINE DQAGI (F, BOUND, INF, EPSABS, EPSREL, RESULT, ABSERR,
   3      +   NEVAL, IER, LIMIT, LENW, LAST, IWORK, WORK)
   4 C***BEGIN PROLOGUE  DQAGI
   5 C***PURPOSE  The routine calculates an approximation result to a given
   6 C            INTEGRAL   I = Integral of F over (BOUND,+INFINITY)
   7 C            OR I = Integral of F over (-INFINITY,BOUND)
   8 C            OR I = Integral of F over (-INFINITY,+INFINITY)
   9 C            Hopefully satisfying following claim for accuracy
  10 C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
  11 C***LIBRARY   SLATEC (QUADPACK)
  12 C***CATEGORY  H2A3A1, H2A4A1
  13 C***TYPE      DOUBLE PRECISION (QAGI-S, DQAGI-D)
  14 C***KEYWORDS  AUTOMATIC INTEGRATOR, EXTRAPOLATION, GENERAL-PURPOSE,
  15 C             GLOBALLY ADAPTIVE, INFINITE INTERVALS, QUADPACK,
  16 C             QUADRATURE, TRANSFORMATION
  17 C***AUTHOR  Piessens, Robert
  18 C             Applied Mathematics and Programming Division
  19 C             K. U. Leuven
  20 C           de Doncker, Elise
  21 C             Applied Mathematics and Programming Division
  22 C             K. U. Leuven
  23 C***DESCRIPTION
  24 C
  25 C        Integration over infinite intervals
  26 C        Standard fortran subroutine
  27 C
  28 C        PARAMETERS
  29 C         ON ENTRY
  30 C            F      - Double precision
  31 C                     Function subprogram defining the integrand
  32 C                     function F(X). The actual name for F needs to be
  33 C                     declared E X T E R N A L in the driver program.
  34 C
  35 C            BOUND  - Double precision
  36 C                     Finite bound of integration range
  37 C                     (has no meaning if interval is doubly-infinite)
  38 C
  39 C            INF    - Integer
  40 C                     indicating the kind of integration range involved
  41 C                     INF = 1 corresponds to  (BOUND,+INFINITY),
  42 C                     INF = -1            to  (-INFINITY,BOUND),
  43 C                     INF = 2             to (-INFINITY,+INFINITY).
  44 C
  45 C            EPSABS - Double precision
  46 C                     Absolute accuracy requested
  47 C            EPSREL - Double precision
  48 C                     Relative accuracy requested
  49 C                     If  EPSABS.LE.0
  50 C                     and EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
  51 C                     the routine will end with IER = 6.
  52 C
  53 C
  54 C         ON RETURN
  55 C            RESULT - Double precision
  56 C                     Approximation to the integral
  57 C
  58 C            ABSERR - Double precision
  59 C                     Estimate of the modulus of the absolute error,
  60 C                     which should equal or exceed ABS(I-RESULT)
  61 C
  62 C            NEVAL  - Integer
  63 C                     Number of integrand evaluations
  64 C
  65 C            IER    - Integer
  66 C                     IER = 0 normal and reliable termination of the
  67 C                             routine. It is assumed that the requested
  68 C                             accuracy has been achieved.
  69 C                   - IER.GT.0 abnormal termination of the routine. The
  70 C                             estimates for result and error are less
  71 C                             reliable. It is assumed that the requested
  72 C                             accuracy has not been achieved.
  73 C            ERROR MESSAGES
  74 C                     IER = 1 Maximum number of subdivisions allowed
  75 C                             has been achieved. One can allow more
  76 C                             subdivisions by increasing the value of
  77 C                             LIMIT (and taking the according dimension
  78 C                             adjustments into account). However, if
  79 C                             this yields no improvement it is advised
  80 C                             to analyze the integrand in order to
  81 C                             determine the integration difficulties. If
  82 C                             the position of a local difficulty can be
  83 C                             determined (e.g. SINGULARITY,
  84 C                             DISCONTINUITY within the interval) one
  85 C                             will probably gain from splitting up the
  86 C                             interval at this point and calling the
  87 C                             integrator on the subranges. If possible,
  88 C                             an appropriate special-purpose integrator
  89 C                             should be used, which is designed for
  90 C                             handling the type of difficulty involved.
  91 C                         = 2 The occurrence of roundoff error is
  92 C                             detected, which prevents the requested
  93 C                             tolerance from being achieved.
  94 C                             The error may be under-estimated.
  95 C                         = 3 Extremely bad integrand behaviour occurs
  96 C                             at some points of the integration
  97 C                             interval.
  98 C                         = 4 The algorithm does not converge.
  99 C                             Roundoff error is detected in the
 100 C                             extrapolation table.
 101 C                             It is assumed that the requested tolerance
 102 C                             cannot be achieved, and that the returned
 103 C                             RESULT is the best which can be obtained.
 104 C                         = 5 The integral is probably divergent, or
 105 C                             slowly convergent. It must be noted that
 106 C                             divergence can occur with any other value
 107 C                             of IER.
 108 C                         = 6 The input is invalid, because
 109 C                             (EPSABS.LE.0 and
 110 C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))
 111 C                              or LIMIT.LT.1 or LENIW.LT.LIMIT*4.
 112 C                             RESULT, ABSERR, NEVAL, LAST are set to
 113 C                             zero.  Except when LIMIT or LENIW is
 114 C                             invalid, IWORK(1), WORK(LIMIT*2+1) and
 115 C                             WORK(LIMIT*3+1) are set to ZERO, WORK(1)
 116 C                             is set to A and WORK(LIMIT+1) to B.
 117 C
 118 C         DIMENSIONING PARAMETERS
 119 C            LIMIT - Integer
 120 C                    Dimensioning parameter for IWORK
 121 C                    LIMIT determines the maximum number of subintervals
 122 C                    in the partition of the given integration interval
 123 C                    (A,B), LIMIT.GE.1.
 124 C                    If LIMIT.LT.1, the routine will end with IER = 6.
 125 C
 126 C            LENW  - Integer
 127 C                    Dimensioning parameter for WORK
 128 C                    LENW must be at least LIMIT*4.
 129 C                    If LENW.LT.LIMIT*4, the routine will end
 130 C                    with IER = 6.
 131 C
 132 C            LAST  - Integer
 133 C                    On return, LAST equals the number of subintervals
 134 C                    produced in the subdivision process, which
 135 C                    determines the number of significant elements
 136 C                    actually in the WORK ARRAYS.
 137 C
 138 C         WORK ARRAYS
 139 C            IWORK - Integer
 140 C                    Vector of dimension at least LIMIT, the first
 141 C                    K elements of which contain pointers
 142 C                    to the error estimates over the subintervals,
 143 C                    such that WORK(LIMIT*3+IWORK(1)),... ,
 144 C                    WORK(LIMIT*3+IWORK(K)) form a decreasing
 145 C                    sequence, with K = LAST if LAST.LE.(LIMIT/2+2), and
 146 C                    K = LIMIT+1-LAST otherwise
 147 C
 148 C            WORK  - Double precision
 149 C                    Vector of dimension at least LENW
 150 C                    on return
 151 C                    WORK(1), ..., WORK(LAST) contain the left
 152 C                     end points of the subintervals in the
 153 C                     partition of (A,B),
 154 C                    WORK(LIMIT+1), ..., WORK(LIMIT+LAST) Contain
 155 C                     the right end points,
 156 C                    WORK(LIMIT*2+1), ...,WORK(LIMIT*2+LAST) contain the
 157 C                     integral approximations over the subintervals,
 158 C                    WORK(LIMIT*3+1), ..., WORK(LIMIT*3)
 159 C                     contain the error estimates.
 160 C
 161 C***REFERENCES  (NONE)
 162 C***ROUTINES CALLED  DQAGIE, XERMSG
 163 C***REVISION HISTORY  (YYMMDD)
 164 C   800101  DATE WRITTEN
 165 C   890831  Modified array declarations.  (WRB)
 166 C   890831  REVISION DATE from Version 3.2
 167 C   891214  Prologue converted to Version 4.0 format.  (BAB)
 168 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
 169 C***END PROLOGUE  DQAGI
 170 C
 171       DOUBLE PRECISION ABSERR,BOUND,EPSABS,EPSREL,F,RESULT,WORK
 172       INTEGER IER,INF,IWORK,LAST,LENW,LIMIT,LVL,L1,L2,L3,NEVAL
 173 C
 174       DIMENSION IWORK(*),WORK(*)
 175 C
 176       EXTERNAL F
 177 C
 178 C         CHECK VALIDITY OF LIMIT AND LENW.
 179 C
 180 C***FIRST EXECUTABLE STATEMENT  DQAGI
 181       IER = 6
 182       NEVAL = 0
 183       LAST = 0
 184       RESULT = 0.0D+00
 185       ABSERR = 0.0D+00
 186       IF(LIMIT.LT.1.OR.LENW.LT.LIMIT*4) GO TO 10
 187 C
 188 C         PREPARE CALL FOR DQAGIE.
 189 C
 190       L1 = LIMIT+1
 191       L2 = LIMIT+L1
 192       L3 = LIMIT+L2
 193 C
 194       CALL DQAGIE(F,BOUND,INF,EPSABS,EPSREL,LIMIT,RESULT,ABSERR,
 195      1  NEVAL,IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST)
 196 C
 197 C         CALL ERROR HANDLER IF NECESSARY.
 198 C
 199        LVL = 0
 200 10    IF(IER.EQ.6) LVL = 1
 201       IF (IER .NE. 0) CALL XERMSG ('SLATEC', 'DQAGI',
 202      +   'ABNORMAL RETURN', IER, LVL)
 203       RETURN
 204       END