dqagi.f

Go to the documentation of this file.

      SUBROUTINE dqagi(F,BOUND,INF,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,
     *   ier,limit,lenw,last,iwork,work)
C***BEGIN PROLOGUE  DQAGI
C***DATE WRITTEN   800101   (YYMMDD)
C***REVISION DATE  830518   (YYMMDD)
C***CATEGORY NO.  H2A3A1,H2A4A1
C***KEYWORDS  AUTOMATIC INTEGRATOR, INFINITE INTERVALS,
C             GENERAL-PURPOSE, TRANSFORMATION, EXTRAPOLATION,
C             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            INTEGRAL   I = INTEGRAL OF F OVER (BOUND,+INFINITY)
C            OR I = INTEGRAL OF F OVER (-INFINITY,BOUND)
C            OR I = INTEGRAL OF F OVER (-INFINITY,+INFINITY)
C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
C***DESCRIPTION
C
C        INTEGRATION OVER INFINITE INTERVALS
C        STANDARD FORTRAN SUBROUTINE
C
C        PARAMETERS
C         ON ENTRY
C            F      - SUBROUTINE F(X,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            BOUND  - DOUBLE PRECISION
C                     FINITE BOUND OF INTEGRATION RANGE
C                     (HAS NO MEANING IF INTERVAL IS DOUBLY-INFINITE)
C
C            INF    - INTEGER
C                     INDICATING THE KIND OF INTEGRATION RANGE INVOLVED
C                     INF = 1 CORRESPONDS TO  (BOUND,+INFINITY),
C                     INF = -1            TO  (-INFINITY,BOUND),
C                     INF = 2             TO (-INFINITY,+INFINITY).
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
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. THE
C                             ESTIMATES FOR RESULT AND ERROR ARE LESS
C                             RELIABLE. IT IS ASSUMED THAT THE REQUESTED
C                             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 (E.G. SINGULARITY,
C                             DISCONTINUITY WITHIN THE INTERVAL) ONE
C                             WILL PROBABLY GAIN FROM SPLITTING UP THE
C                             INTERVAL AT THIS POINT AND CALLING THE
C                             INTEGRATOR ON THE SUBRANGES. IF POSSIBLE,
C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
C                             SHOULD 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 ASSUMED THAT THE REQUESTED TOLERANCE
C                             CANNOT BE ACHIEVED, AND THAT THE RETURNED
C                             RESULT IS THE BEST WHICH 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.
C                         = 6 THE INPUT IS INVALID, BECAUSE
C                             (EPSABS.LE.0 AND
C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))
C                              OR LIMIT.LT.1 OR LENIW.LT.LIMIT*4.
C                             RESULT, ABSERR, NEVAL, LAST ARE SET TO
C                             ZERO. EXEPT WHEN LIMIT OR LENIW IS
C                             INVALID, IWORK(1), WORK(LIMIT*2+1) AND
C                             WORK(LIMIT*3+1) ARE SET TO ZERO, WORK(1)
C                             IS SET TO A AND WORK(LIMIT+1) TO B.
C
C         DIMENSIONING PARAMETERS
C            LIMIT - INTEGER
C                    DIMENSIONING PARAMETER FOR IWORK
C                    LIMIT DETERMINES THE MAXIMUM NUMBER OF SUBINTERVALS
C                    IN THE PARTITION OF THE GIVEN INTEGRATION INTERVAL
C                    (A,B), LIMIT.GE.1.
C                    IF LIMIT.LT.1, THE ROUTINE WILL END WITH IER = 6.
C
C            LENW  - INTEGER
C                    DIMENSIONING PARAMETER FOR WORK
C                    LENW MUST BE AT LEAST LIMIT*4.
C                    IF LENW.LT.LIMIT*4, 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 LIMIT, THE FIRST
C                    K ELEMENTS OF WHICH CONTAIN POINTERS
C                    TO THE ERROR ESTIMATES OVER THE SUBINTERVALS,
C                    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
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 THE
C                     INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS,
C                    WORK(LIMIT*3+1), ..., WORK(LIMIT*3)
C                     CONTAIN THE ERROR ESTIMATES.
C***REFERENCES  (NONE)
C***ROUTINES CALLED  DQAGIE,XERROR
C***END PROLOGUE  DQAGI
C
      DOUBLE PRECISION abserr,bound,epsabs,epsrel,result,work
      INTEGER ier,inf,iwork,last,lenw,limit,lvl,l1,l2,l3,neval
C
      dimension iwork(limit),work(lenw)
C
      EXTERNAL f
C
C         CHECK VALIDITY OF LIMIT AND LENW.
C
C***FIRST EXECUTABLE STATEMENT  DQAGI
      ier = 6
      neval = 0
      last = 0
      result = 0.0d+00
      abserr = 0.0d+00
      IF(limit.LT.1.OR.lenw.LT.limit*4) go to 10
C
C         PREPARE CALL FOR DQAGIE.
C
      l1 = limit+1
      l2 = limit+l1
      l3 = limit+l2
C
      CALL dqagie(f,bound,inf,epsabs,epsrel,limit,result,abserr,
     *  neval,ier,work(1),work(l1),work(l2),work(l3),iwork,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 DQAGI',26,ier,lvl)
      RETURN
      END