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