GNU Octave: libcruft/quadpack/dqk15i.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DQK15I(F,BOUN,INF,A,B,RESULT,ABSERR,RESABS,RESASC,
00002      1   IERR)
00003 C***BEGIN PROLOGUE  DQK15I
00004 C***DATE WRITTEN   800101   (YYMMDD)
00005 C***REVISION DATE  830518   (YYMMDD)
00006 C***CATEGORY NO.  H2A3A2,H2A4A2
00007 C***KEYWORDS  15-POINT TRANSFORMED GAUSS-KRONROD RULES
00008 C***AUTHOR  PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
00009 C           DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
00010 C***PURPOSE  THE ORIGINAL (INFINITE INTEGRATION RANGE IS MAPPED
00011 C            ONTO THE INTERVAL (0,1) AND (A,B) IS A PART OF (0,1).
00012 C            IT IS THE PURPOSE TO COMPUTE
00013 C            I = INTEGRAL OF TRANSFORMED INTEGRAND OVER (A,B),
00014 C            J = INTEGRAL OF ABS(TRANSFORMED INTEGRAND) OVER (A,B).
00015 C***DESCRIPTION
00016 C
00017 C           INTEGRATION RULE
00018 C           STANDARD FORTRAN SUBROUTINE
00019 C           DOUBLE PRECISION VERSION
00020 C
00021 C           PARAMETERS
00022 C            ON ENTRY
00023 C              F      - SUBROUTINE F(X,IERR,RESULT) DEFINING THE INTEGRAND
00024 C                       FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
00025 C                       DECLARED E X T E R N A L IN THE CALLING PROGRAM.
00026 C
00027 C              BOUN   - DOUBLE PRECISION
00028 C                       FINITE BOUND OF ORIGINAL INTEGRATION
00029 C                       RANGE (SET TO ZERO IF INF = +2)
00030 C
00031 C              INF    - INTEGER
00032 C                       IF INF = -1, THE ORIGINAL INTERVAL IS
00033 C                                   (-INFINITY,BOUND),
00034 C                       IF INF = +1, THE ORIGINAL INTERVAL IS
00035 C                                   (BOUND,+INFINITY),
00036 C                       IF INF = +2, THE ORIGINAL INTERVAL IS
00037 C                                   (-INFINITY,+INFINITY) AND
00038 C                       THE INTEGRAL IS COMPUTED AS THE SUM OF TWO
00039 C                       INTEGRALS, ONE OVER (-INFINITY,0) AND ONE OVER
00040 C                       (0,+INFINITY).
00041 C
00042 C              A      - DOUBLE PRECISION
00043 C                       LOWER LIMIT FOR INTEGRATION OVER SUBRANGE
00044 C                       OF (0,1)
00045 C
00046 C              B      - DOUBLE PRECISION
00047 C                       UPPER LIMIT FOR INTEGRATION OVER SUBRANGE
00048 C                       OF (0,1)
00049 C
00050 C            ON RETURN
00051 C              RESULT - DOUBLE PRECISION
00052 C                       APPROXIMATION TO THE INTEGRAL I
00053 C                       RESULT IS COMPUTED BY APPLYING THE 15-POINT
00054 C                       KRONROD RULE(RESK) OBTAINED BY OPTIMAL ADDITION
00055 C                       OF ABSCISSAE TO THE 7-POINT GAUSS RULE(RESG).
00056 C
00057 C              ABSERR - DOUBLE PRECISION
00058 C                       ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
00059 C                       WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
00060 C
00061 C              RESABS - DOUBLE PRECISION
00062 C                       APPROXIMATION TO THE INTEGRAL J
00063 C
00064 C              RESASC - DOUBLE PRECISION
00065 C                       APPROXIMATION TO THE INTEGRAL OF
00066 C                       ABS((TRANSFORMED INTEGRAND)-I/(B-A)) OVER (A,B)
00067 C
00068 C***REFERENCES  (NONE)
00069 C***ROUTINES CALLED  D1MACH
00070 C***END PROLOGUE  DQK15I
00071 C
00072       DOUBLE PRECISION A,ABSC,ABSC1,ABSC2,ABSERR,B,BOUN,CENTR,DABS,DINF,
00073      *  DMAX1,DMIN1,D1MACH,EPMACH,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,
00074      *  RESABS,RESASC,RESG,RESK,RESKH,RESULT,TABSC1,TABSC2,UFLOW,WG,WGK,
00075      *  XGK,FVALT
00076       INTEGER INF,J
00077       EXTERNAL F
00078 C
00079       DIMENSION FV1(7),FV2(7),XGK(8),WGK(8),WG(8)
00080 C
00081 C           THE ABSCISSAE AND WEIGHTS ARE SUPPLIED FOR THE INTERVAL
00082 C           (-1,1).  BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND
00083 C           THEIR CORRESPONDING WEIGHTS ARE GIVEN.
00084 C
00085 C           XGK    - ABSCISSAE OF THE 15-POINT KRONROD RULE
00086 C                    XGK(2), XGK(4), ... ABSCISSAE OF THE 7-POINT
00087 C                    GAUSS RULE
00088 C                    XGK(1), XGK(3), ...  ABSCISSAE WHICH ARE OPTIMALLY
00089 C                    ADDED TO THE 7-POINT GAUSS RULE
00090 C
00091 C           WGK    - WEIGHTS OF THE 15-POINT KRONROD RULE
00092 C
00093 C           WG     - WEIGHTS OF THE 7-POINT GAUSS RULE, CORRESPONDING
00094 C                    TO THE ABSCISSAE XGK(2), XGK(4), ...
00095 C                    WG(1), WG(3), ... ARE SET TO ZERO.
00096 C
00097       DATA WG(1) / 0.0D0 /
00098       DATA WG(2) / 0.1294849661 6886969327 0611432679 082D0 /
00099       DATA WG(3) / 0.0D0 /
00100       DATA WG(4) / 0.2797053914 8927666790 1467771423 780D0 /
00101       DATA WG(5) / 0.0D0 /
00102       DATA WG(6) / 0.3818300505 0511894495 0369775488 975D0 /
00103       DATA WG(7) / 0.0D0 /
00104       DATA WG(8) / 0.4179591836 7346938775 5102040816 327D0 /
00105 C
00106       DATA XGK(1) / 0.9914553711 2081263920 6854697526 329D0 /
00107       DATA XGK(2) / 0.9491079123 4275852452 6189684047 851D0 /
00108       DATA XGK(3) / 0.8648644233 5976907278 9712788640 926D0 /
00109       DATA XGK(4) / 0.7415311855 9939443986 3864773280 788D0 /
00110       DATA XGK(5) / 0.5860872354 6769113029 4144838258 730D0 /
00111       DATA XGK(6) / 0.4058451513 7739716690 6606412076 961D0 /
00112       DATA XGK(7) / 0.2077849550 0789846760 0689403773 245D0 /
00113       DATA XGK(8) / 0.0000000000 0000000000 0000000000 000D0 /
00114 C
00115       DATA WGK(1) / 0.0229353220 1052922496 3732008058 970D0 /
00116       DATA WGK(2) / 0.0630920926 2997855329 0700663189 204D0 /
00117       DATA WGK(3) / 0.1047900103 2225018383 9876322541 518D0 /
00118       DATA WGK(4) / 0.1406532597 1552591874 5189590510 238D0 /
00119       DATA WGK(5) / 0.1690047266 3926790282 6583426598 550D0 /
00120       DATA WGK(6) / 0.1903505780 6478540991 3256402421 014D0 /
00121       DATA WGK(7) / 0.2044329400 7529889241 4161999234 649D0 /
00122       DATA WGK(8) / 0.2094821410 8472782801 2999174891 714D0 /
00123 C
00124 C
00125 C           LIST OF MAJOR VARIABLES
00126 C           -----------------------
00127 C
00128 C           CENTR  - MID POINT OF THE INTERVAL
00129 C           HLGTH  - HALF-LENGTH OF THE INTERVAL
00130 C           ABSC*  - ABSCISSA
00131 C           TABSC* - TRANSFORMED ABSCISSA
00132 C           FVAL*  - FUNCTION VALUE
00133 C           RESG   - RESULT OF THE 7-POINT GAUSS FORMULA
00134 C           RESK   - RESULT OF THE 15-POINT KRONROD FORMULA
00135 C           RESKH  - APPROXIMATION TO THE MEAN VALUE OF THE TRANSFORMED
00136 C                    INTEGRAND OVER (A,B), I.E. TO I/(B-A)
00137 C
00138 C           MACHINE DEPENDENT CONSTANTS
00139 C           ---------------------------
00140 C
00141 C           EPMACH IS THE LARGEST RELATIVE SPACING.
00142 C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
00143 C
00144 C***FIRST EXECUTABLE STATEMENT  DQK15I
00145       EPMACH = D1MACH(4)
00146       UFLOW = D1MACH(1)
00147       DINF = MIN0(1,INF)
00148 C
00149       CENTR = 0.5D+00*(A+B)
00150       HLGTH = 0.5D+00*(B-A)
00151       TABSC1 = BOUN+DINF*(0.1D+01-CENTR)/CENTR
00152       IERR = 0
00153       CALL F(TABSC1,IERR,FVAL1)
00154       IF (IERR .LT. 0) RETURN
00155       IF(INF.EQ.2) THEN
00156         CALL F(-TABSC1,IERR,FVALT)
00157         IF (IERR .LT. 0) RETURN
00158         FVAL1 = FVAL1+FVALT
00159       ENDIF
00160       FC = (FVAL1/CENTR)/CENTR
00161 C
00162 C           COMPUTE THE 15-POINT KRONROD APPROXIMATION TO
00163 C           THE INTEGRAL, AND ESTIMATE THE ERROR.
00164 C
00165       RESG = WG(8)*FC
00166       RESK = WGK(8)*FC
00167       RESABS = DABS(RESK)
00168       DO 10 J=1,7
00169         ABSC = HLGTH*XGK(J)
00170         ABSC1 = CENTR-ABSC
00171         ABSC2 = CENTR+ABSC
00172         TABSC1 = BOUN+DINF*(0.1D+01-ABSC1)/ABSC1
00173         TABSC2 = BOUN+DINF*(0.1D+01-ABSC2)/ABSC2
00174         CALL F(TABSC1,IERR,FVAL1)
00175         IF (IERR .LT. 0) RETURN
00176         CALL F(TABSC2,IERR,FVAL2)
00177         IF (IERR .LT. 0) RETURN
00178         IF(INF.EQ.2) THEN
00179           CALL F(-TABSC1,IERR,FVALT)
00180           IF (IERR .LT. 0) RETURN
00181           FVAL1 = FVAL1+FVALT
00182         ENDIF
00183         IF(INF.EQ.2) THEN
00184           CALL F(-TABSC2,IERR,FVALT)
00185           IF (IERR .LT. 0) RETURN
00186           FVAL2 = FVAL2+FVALT
00187         ENDIF
00188         FVAL1 = (FVAL1/ABSC1)/ABSC1
00189         FVAL2 = (FVAL2/ABSC2)/ABSC2
00190         FV1(J) = FVAL1
00191         FV2(J) = FVAL2
00192         FSUM = FVAL1+FVAL2
00193         RESG = RESG+WG(J)*FSUM
00194         RESK = RESK+WGK(J)*FSUM
00195         RESABS = RESABS+WGK(J)*(DABS(FVAL1)+DABS(FVAL2))
00196    10 CONTINUE
00197       RESKH = RESK*0.5D+00
00198       RESASC = WGK(8)*DABS(FC-RESKH)
00199       DO 20 J=1,7
00200         RESASC = RESASC+WGK(J)*(DABS(FV1(J)-RESKH)+DABS(FV2(J)-RESKH))
00201    20 CONTINUE
00202       RESULT = RESK*HLGTH
00203       RESASC = RESASC*HLGTH
00204       RESABS = RESABS*HLGTH
00205       ABSERR = DABS((RESK-RESG)*HLGTH)
00206       IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.D0) ABSERR = RESASC*
00207      * DMIN1(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00)
00208       IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = DMAX1
00209      * ((EPMACH*0.5D+02)*RESABS,ABSERR)
00210       RETURN
00211       END