00001 SUBROUTINE DQAGP(F,A,B,NPTS2,POINTS,EPSABS,EPSREL,RESULT,ABSERR, 00002 * NEVAL,IER,LENIW,LENW,LAST,IWORK,WORK) 00003 C***BEGIN PROLOGUE DQAGP 00004 C***DATE WRITTEN 800101 (YYMMDD) 00005 C***REVISION DATE 830518 (YYMMDD) 00006 C***CATEGORY NO. H2A2A1 00007 C***KEYWORDS AUTOMATIC INTEGRATOR, GENERAL-PURPOSE, 00008 C SINGULARITIES AT USER SPECIFIED POINTS, 00009 C EXTRAPOLATION, 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 DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B), 00014 C HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY 00015 C BREAK POINTS OF THE INTEGRATION INTERVAL, WHERE LOCAL 00016 C DIFFICULTIES OF THE INTEGRAND MAY OCCUR (E.G. 00017 C SINGULARITIES, DISCONTINUITIES), ARE PROVIDED BY THE USER. 00018 C***DESCRIPTION 00019 C 00020 C COMPUTATION OF A DEFINITE INTEGRAL 00021 C STANDARD FORTRAN SUBROUTINE 00022 C DOUBLE PRECISION VERSION 00023 C 00024 C PARAMETERS 00025 C ON ENTRY 00026 C F - SUBROUTINE F(X,IERR,RESULT) DEFINING THE INTEGRAND 00027 C FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE 00028 C DECLARED E X T E R N A L IN THE DRIVER PROGRAM. 00029 C 00030 C A - DOUBLE PRECISION 00031 C LOWER LIMIT OF INTEGRATION 00032 C 00033 C B - DOUBLE PRECISION 00034 C UPPER LIMIT OF INTEGRATION 00035 C 00036 C NPTS2 - INTEGER 00037 C NUMBER EQUAL TO TWO MORE THAN THE NUMBER OF 00038 C USER-SUPPLIED BREAK POINTS WITHIN THE INTEGRATION 00039 C RANGE, NPTS.GE.2. 00040 C IF NPTS2.LT.2, THE ROUTINE WILL END WITH IER = 6. 00041 C 00042 C POINTS - DOUBLE PRECISION 00043 C VECTOR OF DIMENSION NPTS2, THE FIRST (NPTS2-2) 00044 C ELEMENTS OF WHICH ARE THE USER PROVIDED BREAK 00045 C POINTS. IF THESE POINTS DO NOT CONSTITUTE AN 00046 C ASCENDING SEQUENCE THERE WILL BE AN AUTOMATIC 00047 C SORTING. 00048 C 00049 C EPSABS - DOUBLE PRECISION 00050 C ABSOLUTE ACCURACY REQUESTED 00051 C EPSREL - DOUBLE PRECISION 00052 C RELATIVE ACCURACY REQUESTED 00053 C IF EPSABS.LE.0 00054 C AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28), 00055 C THE ROUTINE WILL END WITH IER = 6. 00056 C 00057 C ON RETURN 00058 C RESULT - DOUBLE PRECISION 00059 C APPROXIMATION TO THE INTEGRAL 00060 C 00061 C ABSERR - DOUBLE PRECISION 00062 C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR, 00063 C WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT) 00064 C 00065 C NEVAL - INTEGER 00066 C NUMBER OF INTEGRAND EVALUATIONS 00067 C 00068 C IER - INTEGER 00069 C IER = 0 NORMAL AND RELIABLE TERMINATION OF THE 00070 C ROUTINE. IT IS ASSUMED THAT THE REQUESTED 00071 C ACCURACY HAS BEEN ACHIEVED. 00072 C IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE. 00073 C THE ESTIMATES FOR INTEGRAL AND ERROR ARE 00074 C LESS RELIABLE. IT IS ASSUMED THAT THE 00075 C REQUESTED ACCURACY HAS NOT BEEN ACHIEVED. 00076 C ERROR MESSAGES 00077 C IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED 00078 C HAS BEEN ACHIEVED. ONE CAN ALLOW MORE 00079 C SUBDIVISIONS BY INCREASING THE VALUE OF 00080 C LIMIT (AND TAKING THE ACCORDING DIMENSION 00081 C ADJUSTMENTS INTO ACCOUNT). HOWEVER, IF 00082 C THIS YIELDS NO IMPROVEMENT IT IS ADVISED 00083 C TO ANALYZE THE INTEGRAND IN ORDER TO 00084 C DETERMINE THE INTEGRATION DIFFICULTIES. IF 00085 C THE POSITION OF A LOCAL DIFFICULTY CAN BE 00086 C DETERMINED (I.E. SINGULARITY, 00087 C DISCONTINUITY WITHIN THE INTERVAL), IT 00088 C SHOULD BE SUPPLIED TO THE ROUTINE AS AN 00089 C ELEMENT OF THE VECTOR POINTS. IF NECESSARY 00090 C AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR 00091 C MUST BE USED, WHICH IS DESIGNED FOR 00092 C HANDLING THE TYPE OF DIFFICULTY INVOLVED. 00093 C = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS 00094 C DETECTED, WHICH PREVENTS THE REQUESTED 00095 C TOLERANCE FROM BEING ACHIEVED. 00096 C THE ERROR MAY BE UNDER-ESTIMATED. 00097 C = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS 00098 C AT SOME POINTS OF THE INTEGRATION 00099 C INTERVAL. 00100 C = 4 THE ALGORITHM DOES NOT CONVERGE. 00101 C ROUNDOFF ERROR IS DETECTED IN THE 00102 C EXTRAPOLATION TABLE. 00103 C IT IS PRESUMED THAT THE REQUESTED 00104 C TOLERANCE CANNOT BE ACHIEVED, AND THAT 00105 C THE RETURNED RESULT IS THE BEST WHICH 00106 C CAN BE OBTAINED. 00107 C = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR 00108 C SLOWLY CONVERGENT. IT MUST BE NOTED THAT 00109 C DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE 00110 C OF IER.GT.0. 00111 C = 6 THE INPUT IS INVALID BECAUSE 00112 C NPTS2.LT.2 OR 00113 C BREAK POINTS ARE SPECIFIED OUTSIDE 00114 C THE INTEGRATION RANGE OR 00115 C (EPSABS.LE.0 AND 00116 C EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28)) 00117 C RESULT, ABSERR, NEVAL, LAST ARE SET TO 00118 C ZERO. EXEPT WHEN LENIW OR LENW OR NPTS2 IS 00119 C INVALID, IWORK(1), IWORK(LIMIT+1), 00120 C WORK(LIMIT*2+1) AND WORK(LIMIT*3+1) 00121 C ARE SET TO ZERO. 00122 C WORK(1) IS SET TO A AND WORK(LIMIT+1) 00123 C TO B (WHERE LIMIT = (LENIW-NPTS2)/2). 00124 C 00125 C DIMENSIONING PARAMETERS 00126 C LENIW - INTEGER 00127 C DIMENSIONING PARAMETER FOR IWORK 00128 C LENIW DETERMINES LIMIT = (LENIW-NPTS2)/2, 00129 C WHICH IS THE MAXIMUM NUMBER OF SUBINTERVALS IN THE 00130 C PARTITION OF THE GIVEN INTEGRATION INTERVAL (A,B), 00131 C LENIW.GE.(3*NPTS2-2). 00132 C IF LENIW.LT.(3*NPTS2-2), THE ROUTINE WILL END WITH 00133 C IER = 6. 00134 C 00135 C LENW - INTEGER 00136 C DIMENSIONING PARAMETER FOR WORK 00137 C LENW MUST BE AT LEAST LENIW*2-NPTS2. 00138 C IF LENW.LT.LENIW*2-NPTS2, THE ROUTINE WILL END 00139 C WITH IER = 6. 00140 C 00141 C LAST - INTEGER 00142 C ON RETURN, LAST EQUALS THE NUMBER OF SUBINTERVALS 00143 C PRODUCED IN THE SUBDIVISION PROCESS, WHICH 00144 C DETERMINES THE NUMBER OF SIGNIFICANT ELEMENTS 00145 C ACTUALLY IN THE WORK ARRAYS. 00146 C 00147 C WORK ARRAYS 00148 C IWORK - INTEGER 00149 C VECTOR OF DIMENSION AT LEAST LENIW. ON RETURN, 00150 C THE FIRST K ELEMENTS OF WHICH CONTAIN 00151 C POINTERS TO THE ERROR ESTIMATES OVER THE 00152 C SUBINTERVALS, SUCH THAT WORK(LIMIT*3+IWORK(1)),..., 00153 C WORK(LIMIT*3+IWORK(K)) FORM A DECREASING 00154 C SEQUENCE, WITH K = LAST IF LAST.LE.(LIMIT/2+2), AND 00155 C K = LIMIT+1-LAST OTHERWISE 00156 C IWORK(LIMIT+1), ...,IWORK(LIMIT+LAST) CONTAIN THE 00157 C SUBDIVISION LEVELS OF THE SUBINTERVALS, I.E. 00158 C IF (AA,BB) IS A SUBINTERVAL OF (P1,P2) 00159 C WHERE P1 AS WELL AS P2 IS A USER-PROVIDED 00160 C BREAK POINT OR INTEGRATION LIMIT, THEN (AA,BB) HAS 00161 C LEVEL L IF ABS(BB-AA) = ABS(P2-P1)*2**(-L), 00162 C IWORK(LIMIT*2+1), ..., IWORK(LIMIT*2+NPTS2) HAVE 00163 C NO SIGNIFICANCE FOR THE USER, 00164 C NOTE THAT LIMIT = (LENIW-NPTS2)/2. 00165 C 00166 C WORK - DOUBLE PRECISION 00167 C VECTOR OF DIMENSION AT LEAST LENW 00168 C ON RETURN 00169 C WORK(1), ..., WORK(LAST) CONTAIN THE LEFT 00170 C END POINTS OF THE SUBINTERVALS IN THE 00171 C PARTITION OF (A,B), 00172 C WORK(LIMIT+1), ..., WORK(LIMIT+LAST) CONTAIN 00173 C THE RIGHT END POINTS, 00174 C WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) CONTAIN 00175 C THE INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS, 00176 C WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST) 00177 C CONTAIN THE CORRESPONDING ERROR ESTIMATES, 00178 C WORK(LIMIT*4+1), ..., WORK(LIMIT*4+NPTS2) 00179 C CONTAIN THE INTEGRATION LIMITS AND THE 00180 C BREAK POINTS SORTED IN AN ASCENDING SEQUENCE. 00181 C NOTE THAT LIMIT = (LENIW-NPTS2)/2. 00182 C 00183 C***REFERENCES (NONE) 00184 C***ROUTINES CALLED DQAGPE,XERROR 00185 C***END PROLOGUE DQAGP 00186 C 00187 DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,POINTS,RESULT,WORK 00188 INTEGER IER,IWORK,LAST,LENIW,LENW,LIMIT,LVL,L1,L2,L3,L4,NEVAL, 00189 * NPTS2 00190 C 00191 DIMENSION IWORK(LENIW),POINTS(NPTS2),WORK(LENW) 00192 C 00193 EXTERNAL F 00194 C 00195 C CHECK VALIDITY OF LIMIT AND LENW. 00196 C 00197 C***FIRST EXECUTABLE STATEMENT DQAGP 00198 IER = 6 00199 NEVAL = 0 00200 LAST = 0 00201 RESULT = 0.0D+00 00202 ABSERR = 0.0D+00 00203 IF(LENIW.LT.(3*NPTS2-2).OR.LENW.LT.(LENIW*2-NPTS2).OR.NPTS2.LT.2) 00204 * GO TO 10 00205 C 00206 C PREPARE CALL FOR DQAGPE. 00207 C 00208 LIMIT = (LENIW-NPTS2)/2 00209 L1 = LIMIT+1 00210 L2 = LIMIT+L1 00211 L3 = LIMIT+L2 00212 L4 = LIMIT+L3 00213 C 00214 CALL DQAGPE(F,A,B,NPTS2,POINTS,EPSABS,EPSREL,LIMIT,RESULT,ABSERR, 00215 * NEVAL,IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),WORK(L4), 00216 * IWORK(1),IWORK(L1),IWORK(L2),LAST) 00217 C 00218 C CALL ERROR HANDLER IF NECESSARY. 00219 C 00220 LVL = 0 00221 10 IF(IER.EQ.6) LVL = 1 00222 IF(IER.GT.0) CALL XERROR('ABNORMAL RETURN FROM DQAGP',26,IER,LVL) 00223 RETURN 00224 END
1.7.1
