d0/df8/dqagp_8f_source.html

      SUBROUTINE dqagp(F,A,B,NPTS2,POINTS,EPSABS,EPSREL,RESULT,ABSERR,

     *   NEVAL,IER,LENIW,LENW,LAST,IWORK,WORK)

C***BEGIN PROLOGUE  DQAGP

C***DATE WRITTEN   800101   (YYMMDD)

C***REVISION DATE  830518   (YYMMDD)

C***CATEGORY NO.  H2A2A1

C***KEYWORDS  AUTOMATIC INTEGRATOR, GENERAL-PURPOSE,

C             SINGULARITIES AT USER SPECIFIED POINTS,

C             EXTRAPOLATION, GLOBALLY ADAPTIVE

C***AUTHOR  PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV - K.U.LEUVEN

C           DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN

C***PURPOSE  THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN

C            DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B),

C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY

C            BREAK POINTS OF THE INTEGRATION INTERVAL, WHERE LOCAL

C            DIFFICULTIES OF THE INTEGRAND MAY OCCUR (E.G.

C            SINGULARITIES, DISCONTINUITIES), ARE PROVIDED BY THE USER.

C***DESCRIPTION

C

C        COMPUTATION OF A DEFINITE INTEGRAL

C        STANDARD FORTRAN SUBROUTINE

C        DOUBLE PRECISION VERSION

C

C        PARAMETERS

C         ON ENTRY

C            F      - SUBROUTINE F(X,IERR,RESULT) DEFINING THE INTEGRAND

C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE

C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.

C

C            A      - DOUBLE PRECISION

C                     LOWER LIMIT OF INTEGRATION

C

C            B      - DOUBLE PRECISION

C                     UPPER LIMIT OF INTEGRATION

C

C            NPTS2  - INTEGER

C                     NUMBER EQUAL TO TWO MORE THAN THE NUMBER OF

C                     USER-SUPPLIED BREAK POINTS WITHIN THE INTEGRATION

C                     RANGE, NPTS.GE.2.

C                     IF NPTS2.LT.2, THE ROUTINE WILL END WITH IER = 6.

C

C            POINTS - DOUBLE PRECISION

C                     VECTOR OF DIMENSION NPTS2, THE FIRST (NPTS2-2)

C                     ELEMENTS OF WHICH ARE THE USER PROVIDED BREAK

C                     POINTS. IF THESE POINTS DO NOT CONSTITUTE AN

C                     ASCENDING SEQUENCE THERE WILL BE AN AUTOMATIC

C                     SORTING.

C

C            EPSABS - DOUBLE PRECISION

C                     ABSOLUTE ACCURACY REQUESTED

C            EPSREL - DOUBLE PRECISION

C                     RELATIVE ACCURACY REQUESTED

C                     IF  EPSABS.LE.0

C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),

C                     THE ROUTINE WILL END WITH IER = 6.

C

C         ON RETURN

C            RESULT - DOUBLE PRECISION

C                     APPROXIMATION TO THE INTEGRAL

C

C            ABSERR - DOUBLE PRECISION

C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,

C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)

C

C            NEVAL  - INTEGER

C                     NUMBER OF INTEGRAND EVALUATIONS

C

C            IER    - INTEGER

C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE

C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED

C                             ACCURACY HAS BEEN ACHIEVED.

C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE.

C                             THE ESTIMATES FOR INTEGRAL AND ERROR ARE

C                             LESS RELIABLE. IT IS ASSUMED THAT THE

C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.

C            ERROR MESSAGES

C                     IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED

C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE

C                             SUBDIVISIONS BY INCREASING THE VALUE OF

C                             LIMIT (AND TAKING THE ACCORDING DIMENSION

C                             ADJUSTMENTS INTO ACCOUNT). HOWEVER, IF

C                             THIS YIELDS NO IMPROVEMENT IT IS ADVISED

C                             TO ANALYZE THE INTEGRAND IN ORDER TO

C                             DETERMINE THE INTEGRATION DIFFICULTIES. IF

C                             THE POSITION OF A LOCAL DIFFICULTY CAN BE

C                             DETERMINED (I.E. SINGULARITY,

C                             DISCONTINUITY WITHIN THE INTERVAL), IT

C                             SHOULD BE SUPPLIED TO THE ROUTINE AS AN

C                             ELEMENT OF THE VECTOR POINTS. IF NECESSARY

C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR

C                             MUST BE USED, WHICH IS DESIGNED FOR

C                             HANDLING THE TYPE OF DIFFICULTY INVOLVED.

C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS

C                             DETECTED, WHICH PREVENTS THE REQUESTED

C                             TOLERANCE FROM BEING ACHIEVED.

C                             THE ERROR MAY BE UNDER-ESTIMATED.

C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS

C                             AT SOME POINTS OF THE INTEGRATION

C                             INTERVAL.

C                         = 4 THE ALGORITHM DOES NOT CONVERGE.

C                             ROUNDOFF ERROR IS DETECTED IN THE

C                             EXTRAPOLATION TABLE.

C                             IT IS PRESUMED THAT THE REQUESTED

C                             TOLERANCE CANNOT BE ACHIEVED, AND THAT

C                             THE RETURNED RESULT IS THE BEST WHICH

C                             CAN BE OBTAINED.

C                         = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR

C                             SLOWLY CONVERGENT. IT MUST BE NOTED THAT

C                             DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE

C                             OF IER.GT.0.

C                         = 6 THE INPUT IS INVALID BECAUSE

C                             NPTS2.LT.2 OR

C                             BREAK POINTS ARE SPECIFIED OUTSIDE

C                             THE INTEGRATION RANGE OR

C                             (EPSABS.LE.0 AND

C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))

C                             RESULT, ABSERR, NEVAL, LAST ARE SET TO

C                             ZERO. EXEPT WHEN LENIW OR LENW OR NPTS2 IS

C                             INVALID, IWORK(1), IWORK(LIMIT+1),

C                             WORK(LIMIT*2+1) AND WORK(LIMIT*3+1)

C                             ARE SET TO ZERO.

C                             WORK(1) IS SET TO A AND WORK(LIMIT+1)

C                             TO B (WHERE LIMIT = (LENIW-NPTS2)/2).

C

C         DIMENSIONING PARAMETERS

C            LENIW - INTEGER

C                    DIMENSIONING PARAMETER FOR IWORK

C                    LENIW DETERMINES LIMIT = (LENIW-NPTS2)/2,

C                    WHICH IS THE MAXIMUM NUMBER OF SUBINTERVALS IN THE

C                    PARTITION OF THE GIVEN INTEGRATION INTERVAL (A,B),

C                    LENIW.GE.(3*NPTS2-2).

C                    IF LENIW.LT.(3*NPTS2-2), THE ROUTINE WILL END WITH

C                    IER = 6.

C

C            LENW  - INTEGER

C                    DIMENSIONING PARAMETER FOR WORK

C                    LENW MUST BE AT LEAST LENIW*2-NPTS2.

C                    IF LENW.LT.LENIW*2-NPTS2, THE ROUTINE WILL END

C                    WITH IER = 6.

C

C            LAST  - INTEGER

C                    ON RETURN, LAST EQUALS THE NUMBER OF SUBINTERVALS

C                    PRODUCED IN THE SUBDIVISION PROCESS, WHICH

C                    DETERMINES THE NUMBER OF SIGNIFICANT ELEMENTS

C                    ACTUALLY IN THE WORK ARRAYS.

C

C         WORK ARRAYS

C            IWORK - INTEGER

C                    VECTOR OF DIMENSION AT LEAST LENIW. ON RETURN,

C                    THE FIRST K ELEMENTS OF WHICH CONTAIN

C                    POINTERS TO THE ERROR ESTIMATES OVER THE

C                    SUBINTERVALS, SUCH THAT WORK(LIMIT*3+IWORK(1)),...,

C                    WORK(LIMIT*3+IWORK(K)) FORM A DECREASING

C                    SEQUENCE, WITH K = LAST IF LAST.LE.(LIMIT/2+2), AND

C                    K = LIMIT+1-LAST OTHERWISE

C                    IWORK(LIMIT+1), ...,IWORK(LIMIT+LAST) CONTAIN THE

C                     SUBDIVISION LEVELS OF THE SUBINTERVALS, I.E.

C                     IF (AA,BB) IS A SUBINTERVAL OF (P1,P2)

C                     WHERE P1 AS WELL AS P2 IS A USER-PROVIDED

C                     BREAK POINT OR INTEGRATION LIMIT, THEN (AA,BB) HAS

C                     LEVEL L IF ABS(BB-AA) = ABS(P2-P1)*2**(-L),

C                    IWORK(LIMIT*2+1), ..., IWORK(LIMIT*2+NPTS2) HAVE

C                     NO SIGNIFICANCE FOR THE USER,

C                    NOTE THAT LIMIT = (LENIW-NPTS2)/2.

C

C            WORK  - DOUBLE PRECISION

C                    VECTOR OF DIMENSION AT LEAST LENW

C                    ON RETURN

C                    WORK(1), ..., WORK(LAST) CONTAIN THE LEFT

C                     END POINTS OF THE SUBINTERVALS IN THE

C                     PARTITION OF (A,B),

C                    WORK(LIMIT+1), ..., WORK(LIMIT+LAST) CONTAIN

C                     THE RIGHT END POINTS,

C                    WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) CONTAIN

C                     THE INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS,

C                    WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST)

C                     CONTAIN THE CORRESPONDING ERROR ESTIMATES,

C                    WORK(LIMIT*4+1), ..., WORK(LIMIT*4+NPTS2)

C                     CONTAIN THE INTEGRATION LIMITS AND THE

C                     BREAK POINTS SORTED IN AN ASCENDING SEQUENCE.

C                    NOTE THAT LIMIT = (LENIW-NPTS2)/2.

C

C***REFERENCES  (NONE)

C***ROUTINES CALLED  DQAGPE,XERROR

C***END PROLOGUE  DQAGP

C

      DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,POINTS,RESULT,WORK

      INTEGER IER,IWORK,LAST,LENIW,LENW,LIMIT,LVL,L1,L2,L3,L4,NEVAL,

     *  npts2

C

      dimension iwork(leniw),points(npts2),work(lenw)

C

      EXTERNAL f

C

C         CHECK VALIDITY OF LIMIT AND LENW.

C

C***FIRST EXECUTABLE STATEMENT  DQAGP

      ier = 6

      neval = 0

      last = 0

      result = 0.0d+00

      abserr = 0.0d+00

      IF(leniw.LT.(3*npts2-2).OR.lenw.LT.(leniw*2-npts2).OR.npts2.LT.2)

     *  GO TO 10

C

C         PREPARE CALL FOR DQAGPE.

C

      limit = (leniw-npts2)/2

      l1 = limit+1

      l2 = limit+l1

      l3 = limit+l2

      l4 = limit+l3

C

      CALL dqagpe(f,a,b,npts2,points,epsabs,epsrel,limit,result,abserr,

     *  neval,ier,work(1),work(l1),work(l2),work(l3),work(l4),

     *  iwork(1),iwork(l1),iwork(l2),last)

C

C         CALL ERROR HANDLER IF NECESSARY.

C

      lvl = 0

10    IF(ier.EQ.6) lvl = 1

      IF(ier.GT.0) CALL xerror('ABNORMAL RETURN FROM DQAGP',26,ier,lvl)

      RETURN


      END

dqagp
subroutine dqagp(f, a, b, npts2, points, epsabs, epsrel, result, abserr, neval, ier, leniw, lenw, last, iwork, work)
Definition dqagp.f:3

dqagpe
subroutine dqagpe(f, a, b, npts2, points, epsabs, epsrel, limit, result, abserr, neval, ier, alist, blist, rlist, elist, pts, iord, level, ndin, last)
Definition dqagpe.f:4

xerror
subroutine xerror(messg, nmessg, nerr, level)
Definition xerror.f:2