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
1.7.1
