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