Mercurial > octave-nkf
view libcruft/quadpack/dqagi.f @ 12134:6c54ad0fde04 release-3-4-x ss-3-3-90
update copyright and version info for release branch
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Sat, 22 Jan 2011 13:50:32 -0500 |
parents | 5b781670e9ee |
children |
line wrap: on
line source
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