src/numerical/slatec/fortran/dqag.f

   1 *DECK DQAG
   2       SUBROUTINE DQAG (F, A, B, EPSABS, EPSREL, KEY, RESULT, ABSERR,
   3      +   NEVAL, IER, LIMIT, LENW, LAST, IWORK, WORK)
   4 C***BEGIN PROLOGUE  DQAG
   5 C***PURPOSE  The routine calculates an approximation result to a given
   6 C            definite integral I = integral of F over (A,B),
   7 C            hopefully satisfying following claim for accuracy
   8 C            ABS(I-RESULT)LE.MAX(EPSABS,EPSREL*ABS(I)).
   9 C***LIBRARY   SLATEC (QUADPACK)
  10 C***CATEGORY  H2A1A1
  11 C***TYPE      DOUBLE PRECISION (QAG-S, DQAG-D)
  12 C***KEYWORDS  AUTOMATIC INTEGRATOR, GAUSS-KRONROD RULES,
  13 C             GENERAL-PURPOSE, GLOBALLY ADAPTIVE, INTEGRAND EXAMINATOR,
  14 C             QUADPACK, QUADRATURE
  15 C***AUTHOR  Piessens, Robert
  16 C             Applied Mathematics and Programming Division
  17 C             K. U. Leuven
  18 C           de Doncker, Elise
  19 C             Applied Mathematics and Programming Division
  20 C             K. U. Leuven
  21 C***DESCRIPTION
  22 C
  23 C        Computation of a definite integral
  24 C        Standard fortran subroutine
  25 C        Double precision version
  26 C
  27 C            F      - Double precision
  28 C                     Function subprogram defining the integrand
  29 C                     Function F(X). The actual name for F needs to be
  30 C                     Declared E X T E R N A L in the driver program.
  31 C
  32 C            A      - Double precision
  33 C                     Lower limit of integration
  34 C
  35 C            B      - Double precision
  36 C                     Upper limit of integration
  37 C
  38 C            EPSABS - Double precision
  39 C                     Absolute accuracy requested
  40 C            EPSREL - Double precision
  41 C                     Relative accuracy requested
  42 C                     If  EPSABS.LE.0
  43 C                     And EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
  44 C                     The routine will end with IER = 6.
  45 C
  46 C            KEY    - Integer
  47 C                     Key for choice of local integration rule
  48 C                     A GAUSS-KRONROD PAIR is used with
  49 C                       7 - 15 POINTS If KEY.LT.2,
  50 C                      10 - 21 POINTS If KEY = 2,
  51 C                      15 - 31 POINTS If KEY = 3,
  52 C                      20 - 41 POINTS If KEY = 4,
  53 C                      25 - 51 POINTS If KEY = 5,
  54 C                      30 - 61 POINTS If KEY.GT.5.
  55 C
  56 C         ON RETURN
  57 C            RESULT - Double precision
  58 C                     Approximation to the integral
  59 C
  60 C            ABSERR - Double precision
  61 C                     Estimate of the modulus of the absolute error,
  62 C                     Which should EQUAL or EXCEED ABS(I-RESULT)
  63 C
  64 C            NEVAL  - Integer
  65 C                     Number of integrand evaluations
  66 C
  67 C            IER    - Integer
  68 C                     IER = 0 Normal and reliable termination of the
  69 C                             routine. It is assumed that the requested
  70 C                             accuracy has been achieved.
  71 C                     IER.GT.0 Abnormal termination of the routine
  72 C                             The estimates for RESULT and ERROR are
  73 C                             Less reliable. It is assumed that the
  74 C                             requested accuracy has not been achieved.
  75 C                      ERROR MESSAGES
  76 C                     IER = 1 Maximum number of subdivisions allowed
  77 C                             has been achieved. One can allow more
  78 C                             subdivisions by increasing the value of
  79 C                             LIMIT (and taking the according dimension
  80 C                             adjustments into account). HOWEVER, If
  81 C                             this yield no improvement it is advised
  82 C                             to analyze the integrand in order to
  83 C                             determine the integration difficulties.
  84 C                             If the position of a local difficulty can
  85 C                             be determined (I.E. SINGULARITY,
  86 C                             DISCONTINUITY WITHIN THE INTERVAL) One
  87 C                             will probably gain from splitting up the
  88 C                             interval at this point and calling the
  89 C                             INTEGRATOR on the SUBRANGES. If possible,
  90 C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
  91 C                             should be used which is designed for
  92 C                             handling the type of difficulty involved.
  93 C                         = 2 The occurrence of roundoff error is
  94 C                             detected, which prevents the requested
  95 C                             tolerance from being achieved.
  96 C                         = 3 Extremely bad integrand behaviour occurs
  97 C                             at some points of the integration
  98 C                             interval.
  99 C                         = 6 The input is invalid, because
 100 C                             (EPSABS.LE.0 AND
 101 C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))
 102 C                             OR LIMIT.LT.1 OR LENW.LT.LIMIT*4.
 103 C                             RESULT, ABSERR, NEVAL, LAST are set
 104 C                             to zero.
 105 C                             EXCEPT when LENW is invalid, IWORK(1),
 106 C                             WORK(LIMIT*2+1) and WORK(LIMIT*3+1) are
 107 C                             set to zero, WORK(1) is set to A and
 108 C                             WORK(LIMIT+1) to B.
 109 C
 110 C         DIMENSIONING PARAMETERS
 111 C            LIMIT - Integer
 112 C                    Dimensioning parameter for IWORK
 113 C                    Limit determines the maximum number of subintervals
 114 C                    in the partition of the given integration interval
 115 C                    (A,B), LIMIT.GE.1.
 116 C                    If LIMIT.LT.1, the routine will end with IER = 6.
 117 C
 118 C            LENW  - Integer
 119 C                    Dimensioning parameter for work
 120 C                    LENW must be at least LIMIT*4.
 121 C                    IF LENW.LT.LIMIT*4, the routine will end with
 122 C                    IER = 6.
 123 C
 124 C            LAST  - Integer
 125 C                    On return, LAST equals the number of subintervals
 126 C                    produced in the subdivision process, which
 127 C                    determines the number of significant elements
 128 C                    actually in the WORK ARRAYS.
 129 C
 130 C         WORK ARRAYS
 131 C            IWORK - Integer
 132 C                    Vector of dimension at least limit, the first K
 133 C                    elements of which contain pointers to the error
 134 C                    estimates over the subintervals, such that
 135 C                    WORK(LIMIT*3+IWORK(1)),... , WORK(LIMIT*3+IWORK(K))
 136 C                    form a decreasing sequence with K = LAST If
 137 C                    LAST.LE.(LIMIT/2+2), and K = LIMIT+1-LAST otherwise
 138 C
 139 C            WORK  - Double precision
 140 C                    Vector of dimension at least LENW
 141 C                    on return
 142 C                    WORK(1), ..., WORK(LAST) contain the left end
 143 C                    points of the subintervals in the partition of
 144 C                     (A,B),
 145 C                    WORK(LIMIT+1), ..., WORK(LIMIT+LAST) contain the
 146 C                     right end points,
 147 C                    WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) contain
 148 C                     the integral approximations over the subintervals,
 149 C                    WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST) contain
 150 C                     the error estimates.
 151 C
 152 C***REFERENCES  (NONE)
 153 C***ROUTINES CALLED  DQAGE, XERMSG
 154 C***REVISION HISTORY  (YYMMDD)
 155 C   800101  DATE WRITTEN
 156 C   890831  Modified array declarations.  (WRB)
 157 C   890831  REVISION DATE from Version 3.2
 158 C   891214  Prologue converted to Version 4.0 format.  (BAB)
 159 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
 160 C***END PROLOGUE  DQAG
 161       DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK
 162       INTEGER IER,IWORK,KEY,LAST,LENW,LIMIT,LVL,L1,L2,L3,NEVAL
 163 C
 164       DIMENSION IWORK(*),WORK(*)
 165 C
 166       EXTERNAL F
 167 C***FIRST EXECUTABLE STATEMENT  DQAG
 168       IER = 6
 169       NEVAL = 0
 170       LAST = 0
 171       RESULT = 0.0D+00
 172       ABSERR = 0.0D+00
 173       IF (LIMIT.GE.1 .AND. LENW.GE.LIMIT*4) THEN
 174 C
 175 C        PREPARE CALL FOR DQAGE.
 176 C
 177          L1 = LIMIT+1
 178          L2 = LIMIT+L1
 179          L3 = LIMIT+L2
 180 C
 181          CALL DQAGE(F,A,B,EPSABS,EPSREL,KEY,LIMIT,RESULT,ABSERR,NEVAL,
 182      1     IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST)
 183 C
 184 C        CALL ERROR HANDLER IF NECESSARY.
 185 C
 186          LVL = 0
 187       ENDIF
 188 C
 189       IF (IER.EQ.6) LVL = 1
 190       IF (IER .NE. 0) CALL XERMSG ('SLATEC', 'DQAG',
 191      +   'ABNORMAL RETURN', IER, LVL)
 192       RETURN
 193       END