src/numerical/slatec/fortran/dqagp.f

   1 *DECK DQAGP
   2       SUBROUTINE DQAGP (F, A, B, NPTS2, POINTS, EPSABS, EPSREL, RESULT,
   3      +   ABSERR, NEVAL, IER, LENIW, LENW, LAST, IWORK, WORK)
   4 C***BEGIN PROLOGUE  DQAGP
   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            break points of the integration interval, where local
   9 C            difficulties of the integrand may occur (e.g.
  10 C            SINGULARITIES, DISCONTINUITIES), are provided by the user.
  11 C***LIBRARY   SLATEC (QUADPACK)
  12 C***CATEGORY  H2A2A1
  13 C***TYPE      DOUBLE PRECISION (QAGP-S, DQAGP-D)
  14 C***KEYWORDS  AUTOMATIC INTEGRATOR, EXTRAPOLATION, GENERAL-PURPOSE,
  15 C             GLOBALLY ADAPTIVE, QUADPACK, QUADRATURE,
  16 C             SINGULARITIES AT USER SPECIFIED POINTS
  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        Computation of a definite integral
  26 C        Standard fortran subroutine
  27 C        Double precision version
  28 C
  29 C        PARAMETERS
  30 C         ON ENTRY
  31 C            F      - Double precision
  32 C                     Function subprogram defining the integrand
  33 C                     Function F(X). The actual name for F needs to be
  34 C                     declared E X T E R N A L in the driver program.
  35 C
  36 C            A      - Double precision
  37 C                     Lower limit of integration
  38 C
  39 C            B      - Double precision
  40 C                     Upper limit of integration
  41 C
  42 C            NPTS2  - Integer
  43 C                     Number equal to two more than the number of
  44 C                     user-supplied break points within the integration
  45 C                     range, NPTS.GE.2.
  46 C                     If NPTS2.LT.2, The routine will end with IER = 6.
  47 C
  48 C            POINTS - Double precision
  49 C                     Vector of dimension NPTS2, the first (NPTS2-2)
  50 C                     elements of which are the user provided break
  51 C                     points. If these points do not constitute an
  52 C                     ascending sequence there will be an automatic
  53 C                     sorting.
  54 C
  55 C            EPSABS - Double precision
  56 C                     Absolute accuracy requested
  57 C            EPSREL - Double precision
  58 C                     Relative accuracy requested
  59 C                     If  EPSABS.LE.0
  60 C                     And EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
  61 C                     The routine will end with IER = 6.
  62 C
  63 C         ON RETURN
  64 C            RESULT - Double precision
  65 C                     Approximation to the integral
  66 C
  67 C            ABSERR - Double precision
  68 C                     Estimate of the modulus of the absolute error,
  69 C                     which should equal or exceed ABS(I-RESULT)
  70 C
  71 C            NEVAL  - Integer
  72 C                     Number of integrand evaluations
  73 C
  74 C            IER    - Integer
  75 C                     IER = 0 Normal and reliable termination of the
  76 C                             routine. It is assumed that the requested
  77 C                             accuracy has been achieved.
  78 C                     IER.GT.0 Abnormal termination of the routine.
  79 C                             The estimates for integral and error are
  80 C                             less reliable. it is assumed that the
  81 C                             requested accuracy has not been achieved.
  82 C            ERROR MESSAGES
  83 C                     IER = 1 Maximum number of subdivisions allowed
  84 C                             has been achieved. one can allow more
  85 C                             subdivisions by increasing the value of
  86 C                             LIMIT (and taking the according dimension
  87 C                             adjustments into account). However, if
  88 C                             this yields no improvement it is advised
  89 C                             to analyze the integrand in order to
  90 C                             determine the integration difficulties. If
  91 C                             the position of a local difficulty can be
  92 C                             determined (i.e. SINGULARITY,
  93 C                             DISCONTINUITY within the interval), it
  94 C                             should be supplied to the routine as an
  95 C                             element of the vector points. If necessary
  96 C                             an appropriate special-purpose integrator
  97 C                             must be used, which is designed for
  98 C                             handling the type of difficulty involved.
  99 C                         = 2 The occurrence of roundoff error is
 100 C                             detected, which prevents the requested
 101 C                             tolerance from being achieved.
 102 C                             The error may be under-estimated.
 103 C                         = 3 Extremely bad integrand behaviour occurs
 104 C                             at some points of the integration
 105 C                             interval.
 106 C                         = 4 The algorithm does not converge.
 107 C                             roundoff error is detected in the
 108 C                             extrapolation table.
 109 C                             It is presumed that the requested
 110 C                             tolerance cannot be achieved, and that
 111 C                             the returned RESULT is the best which
 112 C                             can be obtained.
 113 C                         = 5 The integral is probably divergent, or
 114 C                             slowly convergent. it must be noted that
 115 C                             divergence can occur with any other value
 116 C                             of IER.GT.0.
 117 C                         = 6 The input is invalid because
 118 C                             NPTS2.LT.2 or
 119 C                             break points are specified outside
 120 C                             the integration range or
 121 C                             (EPSABS.LE.0 and
 122 C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))
 123 C                             RESULT, ABSERR, NEVAL, LAST are set to
 124 C                             zero.  Except when LENIW or LENW or NPTS2
 125 C                             is invalid, IWORK(1), IWORK(LIMIT+1),
 126 C                             WORK(LIMIT*2+1) and WORK(LIMIT*3+1)
 127 C                             are set to zero.
 128 C                             WORK(1) is set to A and WORK(LIMIT+1)
 129 C                             to B (where LIMIT = (LENIW-NPTS2)/2).
 130 C
 131 C         DIMENSIONING PARAMETERS
 132 C            LENIW - Integer
 133 C                    Dimensioning parameter for IWORK
 134 C                    LENIW determines LIMIT = (LENIW-NPTS2)/2,
 135 C                    which is the maximum number of subintervals in the
 136 C                    partition of the given integration interval (A,B),
 137 C                    LENIW.GE.(3*NPTS2-2).
 138 C                    If LENIW.LT.(3*NPTS2-2), the routine will end with
 139 C                    IER = 6.
 140 C
 141 C            LENW  - Integer
 142 C                    Dimensioning parameter for WORK
 143 C                    LENW must be at least LENIW*2-NPTS2.
 144 C                    If LENW.LT.LENIW*2-NPTS2, the routine will end
 145 C                    with IER = 6.
 146 C
 147 C            LAST  - Integer
 148 C                    On return, LAST equals the number of subintervals
 149 C                    produced in the subdivision process, which
 150 C                    determines the number of significant elements
 151 C                    actually in the WORK ARRAYS.
 152 C
 153 C         WORK ARRAYS
 154 C            IWORK - Integer
 155 C                    Vector of dimension at least LENIW. on return,
 156 C                    the first K elements of which contain
 157 C                    pointers to the error estimates over the
 158 C                    subintervals, such that WORK(LIMIT*3+IWORK(1)),...,
 159 C                    WORK(LIMIT*3+IWORK(K)) form a decreasing
 160 C                    sequence, with K = LAST if LAST.LE.(LIMIT/2+2), and
 161 C                    K = LIMIT+1-LAST otherwise
 162 C                    IWORK(LIMIT+1), ...,IWORK(LIMIT+LAST) Contain the
 163 C                     subdivision levels of the subintervals, i.e.
 164 C                     if (AA,BB) is a subinterval of (P1,P2)
 165 C                     where P1 as well as P2 is a user-provided
 166 C                     break point or integration LIMIT, then (AA,BB) has
 167 C                     level L if ABS(BB-AA) = ABS(P2-P1)*2**(-L),
 168 C                    IWORK(LIMIT*2+1), ..., IWORK(LIMIT*2+NPTS2) have
 169 C                     no significance for the user,
 170 C                    note that LIMIT = (LENIW-NPTS2)/2.
 171 C
 172 C            WORK  - Double precision
 173 C                    Vector of dimension at least LENW
 174 C                    on return
 175 C                    WORK(1), ..., WORK(LAST) contain the left
 176 C                     end points of the subintervals in the
 177 C                     partition of (A,B),
 178 C                    WORK(LIMIT+1), ..., WORK(LIMIT+LAST) contain
 179 C                     the right end points,
 180 C                    WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) contain
 181 C                     the integral approximations over the subintervals,
 182 C                    WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST)
 183 C                     contain the corresponding error estimates,
 184 C                    WORK(LIMIT*4+1), ..., WORK(LIMIT*4+NPTS2)
 185 C                     contain the integration limits and the
 186 C                     break points sorted in an ascending sequence.
 187 C                    note that LIMIT = (LENIW-NPTS2)/2.
 188 C
 189 C***REFERENCES  (NONE)
 190 C***ROUTINES CALLED  DQAGPE, XERMSG
 191 C***REVISION HISTORY  (YYMMDD)
 192 C   800101  DATE WRITTEN
 193 C   890831  Modified array declarations.  (WRB)
 194 C   890831  REVISION DATE from Version 3.2
 195 C   891214  Prologue converted to Version 4.0 format.  (BAB)
 196 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
 197 C***END PROLOGUE  DQAGP
 198 C
 199       DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,F,POINTS,RESULT,WORK
 200       INTEGER IER,IWORK,LAST,LENIW,LENW,LIMIT,LVL,L1,L2,L3,L4,NEVAL,
 201      1  NPTS2
 202 C
 203       DIMENSION IWORK(*),POINTS(*),WORK(*)
 204 C
 205       EXTERNAL F
 206 C
 207 C         CHECK VALIDITY OF LIMIT AND LENW.
 208 C
 209 C***FIRST EXECUTABLE STATEMENT  DQAGP
 210       IER = 6
 211       NEVAL = 0
 212       LAST = 0
 213       RESULT = 0.0D+00
 214       ABSERR = 0.0D+00
 215       IF(LENIW.LT.(3*NPTS2-2).OR.LENW.LT.(LENIW*2-NPTS2).OR.NPTS2.LT.2)
 216      1  GO TO 10
 217 C
 218 C         PREPARE CALL FOR DQAGPE.
 219 C
 220       LIMIT = (LENIW-NPTS2)/2
 221       L1 = LIMIT+1
 222       L2 = LIMIT+L1
 223       L3 = LIMIT+L2
 224       L4 = LIMIT+L3
 225 C
 226       CALL DQAGPE(F,A,B,NPTS2,POINTS,EPSABS,EPSREL,LIMIT,RESULT,ABSERR,
 227      1  NEVAL,IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),WORK(L4),
 228      2  IWORK(1),IWORK(L1),IWORK(L2),LAST)
 229 C
 230 C         CALL ERROR HANDLER IF NECESSARY.
 231 C
 232       LVL = 0
 233 10    IF(IER.EQ.6) LVL = 1
 234       IF (IER .NE. 0) CALL XERMSG ('SLATEC', 'DQAGP',
 235      +   'ABNORMAL RETURN', IER, LVL)
 236       RETURN
 237       END