dqagi.f

Go to the documentation of this file.
00001       SUBROUTINE DQAGI(F,BOUND,INF,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,
00002      *   IER,LIMIT,LENW,LAST,IWORK,WORK)
00003 C***BEGIN PROLOGUE  DQAGI
00004 C***DATE WRITTEN   800101   (YYMMDD)
00005 C***REVISION DATE  830518   (YYMMDD)
00006 C***CATEGORY NO.  H2A3A1,H2A4A1
00007 C***KEYWORDS  AUTOMATIC INTEGRATOR, INFINITE INTERVALS,
00008 C             GENERAL-PURPOSE, TRANSFORMATION, EXTRAPOLATION,
00009 C             GLOBALLY ADAPTIVE
00010 C***AUTHOR  PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
00011 C           DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. -K.U.LEUVEN
00012 C***PURPOSE  THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
00013 C            INTEGRAL   I = INTEGRAL OF F OVER (BOUND,+INFINITY)
00014 C            OR I = INTEGRAL OF F OVER (-INFINITY,BOUND)
00015 C            OR I = INTEGRAL OF F OVER (-INFINITY,+INFINITY)
00016 C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
00017 C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
00018 C***DESCRIPTION
00019 C
00020 C        INTEGRATION OVER INFINITE INTERVALS
00021 C        STANDARD FORTRAN SUBROUTINE
00022 C
00023 C        PARAMETERS
00024 C         ON ENTRY
00025 C            F      - SUBROUTINE F(X,RESULT) DEFINING THE INTEGRAND
00026 C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
00027 C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
00028 C
00029 C            BOUND  - DOUBLE PRECISION
00030 C                     FINITE BOUND OF INTEGRATION RANGE
00031 C                     (HAS NO MEANING IF INTERVAL IS DOUBLY-INFINITE)
00032 C
00033 C            INF    - INTEGER
00034 C                     INDICATING THE KIND OF INTEGRATION RANGE INVOLVED
00035 C                     INF = 1 CORRESPONDS TO  (BOUND,+INFINITY),
00036 C                     INF = -1            TO  (-INFINITY,BOUND),
00037 C                     INF = 2             TO (-INFINITY,+INFINITY).
00038 C
00039 C            EPSABS - DOUBLE PRECISION
00040 C                     ABSOLUTE ACCURACY REQUESTED
00041 C            EPSREL - DOUBLE PRECISION
00042 C                     RELATIVE ACCURACY REQUESTED
00043 C                     IF  EPSABS.LE.0
00044 C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
00045 C                     THE ROUTINE WILL END WITH IER = 6.
00046 C
00047 C
00048 C         ON RETURN
00049 C            RESULT - DOUBLE PRECISION
00050 C                     APPROXIMATION TO THE INTEGRAL
00051 C
00052 C            ABSERR - DOUBLE PRECISION
00053 C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
00054 C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
00055 C
00056 C            NEVAL  - INTEGER
00057 C                     NUMBER OF INTEGRAND EVALUATIONS
00058 C
00059 C            IER    - INTEGER
00060 C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
00061 C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
00062 C                             ACCURACY HAS BEEN ACHIEVED.
00063 C                   - IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE. THE
00064 C                             ESTIMATES FOR RESULT AND ERROR ARE LESS
00065 C                             RELIABLE. IT IS ASSUMED THAT THE REQUESTED
00066 C                             ACCURACY HAS NOT BEEN ACHIEVED.
00067 C            ERROR MESSAGES
00068 C                     IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
00069 C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE
00070 C                             SUBDIVISIONS BY INCREASING THE VALUE OF
00071 C                             LIMIT (AND TAKING THE ACCORDING DIMENSION
00072 C                             ADJUSTMENTS INTO ACCOUNT). HOWEVER, IF
00073 C                             THIS YIELDS NO IMPROVEMENT IT IS ADVISED
00074 C                             TO ANALYZE THE INTEGRAND IN ORDER TO
00075 C                             DETERMINE THE INTEGRATION DIFFICULTIES. IF
00076 C                             THE POSITION OF A LOCAL DIFFICULTY CAN BE
00077 C                             DETERMINED (E.G. SINGULARITY,
00078 C                             DISCONTINUITY WITHIN THE INTERVAL) ONE
00079 C                             WILL PROBABLY GAIN FROM SPLITTING UP THE
00080 C                             INTERVAL AT THIS POINT AND CALLING THE
00081 C                             INTEGRATOR ON THE SUBRANGES. IF POSSIBLE,
00082 C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
00083 C                             SHOULD BE USED, WHICH IS DESIGNED FOR
00084 C                             HANDLING THE TYPE OF DIFFICULTY INVOLVED.
00085 C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS
00086 C                             DETECTED, WHICH PREVENTS THE REQUESTED
00087 C                             TOLERANCE FROM BEING ACHIEVED.
00088 C                             THE ERROR MAY BE UNDER-ESTIMATED.
00089 C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS
00090 C                             AT SOME POINTS OF THE INTEGRATION
00091 C                             INTERVAL.
00092 C                         = 4 THE ALGORITHM DOES NOT CONVERGE.
00093 C                             ROUNDOFF ERROR IS DETECTED IN THE
00094 C                             EXTRAPOLATION TABLE.
00095 C                             IT IS ASSUMED THAT THE REQUESTED TOLERANCE
00096 C                             CANNOT BE ACHIEVED, AND THAT THE RETURNED
00097 C                             RESULT IS THE BEST WHICH CAN BE OBTAINED.
00098 C                         = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
00099 C                             SLOWLY CONVERGENT. IT MUST BE NOTED THAT
00100 C                             DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE
00101 C                             OF IER.
00102 C                         = 6 THE INPUT IS INVALID, BECAUSE
00103 C                             (EPSABS.LE.0 AND
00104 C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))
00105 C                              OR LIMIT.LT.1 OR LENIW.LT.LIMIT*4.
00106 C                             RESULT, ABSERR, NEVAL, LAST ARE SET TO
00107 C                             ZERO. EXEPT WHEN LIMIT OR LENIW IS
00108 C                             INVALID, IWORK(1), WORK(LIMIT*2+1) AND
00109 C                             WORK(LIMIT*3+1) ARE SET TO ZERO, WORK(1)
00110 C                             IS SET TO A AND WORK(LIMIT+1) TO B.
00111 C
00112 C         DIMENSIONING PARAMETERS
00113 C            LIMIT - INTEGER
00114 C                    DIMENSIONING PARAMETER FOR IWORK
00115 C                    LIMIT DETERMINES THE MAXIMUM NUMBER OF SUBINTERVALS
00116 C                    IN THE PARTITION OF THE GIVEN INTEGRATION INTERVAL
00117 C                    (A,B), LIMIT.GE.1.
00118 C                    IF LIMIT.LT.1, THE ROUTINE WILL END WITH IER = 6.
00119 C
00120 C            LENW  - INTEGER
00121 C                    DIMENSIONING PARAMETER FOR WORK
00122 C                    LENW MUST BE AT LEAST LIMIT*4.
00123 C                    IF LENW.LT.LIMIT*4, THE ROUTINE WILL END
00124 C                    WITH IER = 6.
00125 C
00126 C            LAST  - INTEGER
00127 C                    ON RETURN, LAST EQUALS THE NUMBER OF SUBINTERVALS
00128 C                    PRODUCED IN THE SUBDIVISION PROCESS, WHICH
00129 C                    DETERMINES THE NUMBER OF SIGNIFICANT ELEMENTS
00130 C                    ACTUALLY IN THE WORK ARRAYS.
00131 C
00132 C         WORK ARRAYS
00133 C            IWORK - INTEGER
00134 C                    VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
00135 C                    K ELEMENTS OF WHICH CONTAIN POINTERS
00136 C                    TO THE ERROR ESTIMATES OVER THE SUBINTERVALS,
00137 C                    SUCH THAT WORK(LIMIT*3+IWORK(1)),... ,
00138 C                    WORK(LIMIT*3+IWORK(K)) FORM A DECREASING
00139 C                    SEQUENCE, WITH K = LAST IF LAST.LE.(LIMIT/2+2), AND
00140 C                    K = LIMIT+1-LAST OTHERWISE
00141 C
00142 C            WORK  - DOUBLE PRECISION
00143 C                    VECTOR OF DIMENSION AT LEAST LENW
00144 C                    ON RETURN
00145 C                    WORK(1), ..., WORK(LAST) CONTAIN THE LEFT
00146 C                     END POINTS OF THE SUBINTERVALS IN THE
00147 C                     PARTITION OF (A,B),
00148 C                    WORK(LIMIT+1), ..., WORK(LIMIT+LAST) CONTAIN
00149 C                     THE RIGHT END POINTS,
00150 C                    WORK(LIMIT*2+1), ...,WORK(LIMIT*2+LAST) CONTAIN THE
00151 C                     INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS,
00152 C                    WORK(LIMIT*3+1), ..., WORK(LIMIT*3)
00153 C                     CONTAIN THE ERROR ESTIMATES.
00154 C***REFERENCES  (NONE)
00155 C***ROUTINES CALLED  DQAGIE,XERROR
00156 C***END PROLOGUE  DQAGI
00157 C
00158       DOUBLE PRECISION ABSERR,BOUND,EPSABS,EPSREL,RESULT,WORK
00159       INTEGER IER,INF,IWORK,LAST,LENW,LIMIT,LVL,L1,L2,L3,NEVAL
00160 C
00161       DIMENSION IWORK(LIMIT),WORK(LENW)
00162 C
00163       EXTERNAL F
00164 C
00165 C         CHECK VALIDITY OF LIMIT AND LENW.
00166 C
00167 C***FIRST EXECUTABLE STATEMENT  DQAGI
00168       IER = 6
00169       NEVAL = 0
00170       LAST = 0
00171       RESULT = 0.0D+00
00172       ABSERR = 0.0D+00
00173       IF(LIMIT.LT.1.OR.LENW.LT.LIMIT*4) GO TO 10
00174 C
00175 C         PREPARE CALL FOR DQAGIE.
00176 C
00177       L1 = LIMIT+1
00178       L2 = LIMIT+L1
00179       L3 = LIMIT+L2
00180 C
00181       CALL DQAGIE(F,BOUND,INF,EPSABS,EPSREL,LIMIT,RESULT,ABSERR,
00182      *  NEVAL,IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST)
00183 C
00184 C         CALL ERROR HANDLER IF NECESSARY.
00185 C
00186        LVL = 0
00187 10    IF(IER.EQ.6) LVL = 1
00188       IF(IER.GT.0) CALL XERROR('ABNORMAL RETURN FROM DQAGI',26,IER,LVL)
00189       RETURN
00190       END
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines