dqk15i.f

Go to the documentation of this file.

      SUBROUTINE dqk15i(F,BOUN,INF,A,B,RESULT,ABSERR,RESABS,RESASC,
     1   IERR)
C***BEGIN PROLOGUE  DQK15I
C***DATE WRITTEN   800101   (YYMMDD)
C***REVISION DATE  830518   (YYMMDD)
C***CATEGORY NO.  H2A3A2,H2A4A2
C***KEYWORDS  15-POINT TRANSFORMED GAUSS-KRONROD RULES
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 ORIGINAL (INFINITE INTEGRATION RANGE IS MAPPED
C            ONTO THE INTERVAL (0,1) AND (A,B) IS A PART OF (0,1).
C            IT IS THE PURPOSE TO COMPUTE
C            I = INTEGRAL OF TRANSFORMED INTEGRAND OVER (A,B),
C            J = INTEGRAL OF ABS(TRANSFORMED INTEGRAND) OVER (A,B).
C***DESCRIPTION
C
C           INTEGRATION RULE
C           STANDARD FORTRAN SUBROUTINE
C           DOUBLE PRECISION VERSION
C
C           PARAMETERS
C            ON ENTRY
C              F      - SUBROUTINE F(X,IERR,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 CALLING PROGRAM.
C
C              BOUN   - DOUBLE PRECISION
C                       FINITE BOUND OF ORIGINAL INTEGRATION
C                       RANGE (SET TO ZERO IF INF = +2)
C
C              INF    - INTEGER
C                       IF INF = -1, THE ORIGINAL INTERVAL IS
C                                   (-INFINITY,BOUND),
C                       IF INF = +1, THE ORIGINAL INTERVAL IS
C                                   (BOUND,+INFINITY),
C                       IF INF = +2, THE ORIGINAL INTERVAL IS
C                                   (-INFINITY,+INFINITY) AND
C                       THE INTEGRAL IS COMPUTED AS THE SUM OF TWO
C                       INTEGRALS, ONE OVER (-INFINITY,0) AND ONE OVER
C                       (0,+INFINITY).
C
C              A      - DOUBLE PRECISION
C                       LOWER LIMIT FOR INTEGRATION OVER SUBRANGE
C                       OF (0,1)
C
C              B      - DOUBLE PRECISION
C                       UPPER LIMIT FOR INTEGRATION OVER SUBRANGE
C                       OF (0,1)
C
C            ON RETURN
C              RESULT - DOUBLE PRECISION
C                       APPROXIMATION TO THE INTEGRAL I
C                       RESULT IS COMPUTED BY APPLYING THE 15-POINT
C                       KRONROD RULE(RESK) OBTAINED BY OPTIMAL ADDITION
C                       OF ABSCISSAE TO THE 7-POINT GAUSS RULE(RESG).
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              RESABS - DOUBLE PRECISION
C                       APPROXIMATION TO THE INTEGRAL J
C
C              RESASC - DOUBLE PRECISION
C                       APPROXIMATION TO THE INTEGRAL OF
C                       ABS((TRANSFORMED INTEGRAND)-I/(B-A)) OVER (A,B)
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH
C***END PROLOGUE  DQK15I
C
      DOUBLE PRECISION A,ABSC,ABSC1,ABSC2,ABSERR,B,BOUN,CENTR,DABS,DINF,
     *  dmax1,dmin1,d1mach,epmach,fc,fsum,fval1,fval2,fv1,fv2,hlgth,
     *  resabs,resasc,resg,resk,reskh,result,tabsc1,tabsc2,uflow,wg,wgk,
     *  xgk,fvalt
      INTEGER INF,J
      EXTERNAL f
C
      dimension fv1(7),fv2(7),xgk(8),wgk(8),wg(8)
C
C           THE ABSCISSAE AND WEIGHTS ARE SUPPLIED FOR THE INTERVAL
C           (-1,1).  BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND
C           THEIR CORRESPONDING WEIGHTS ARE GIVEN.
C
C           XGK    - ABSCISSAE OF THE 15-POINT KRONROD RULE
C                    XGK(2), XGK(4), ... ABSCISSAE OF THE 7-POINT
C                    GAUSS RULE
C                    XGK(1), XGK(3), ...  ABSCISSAE WHICH ARE OPTIMALLY
C                    ADDED TO THE 7-POINT GAUSS RULE
C
C           WGK    - WEIGHTS OF THE 15-POINT KRONROD RULE
C
C           WG     - WEIGHTS OF THE 7-POINT GAUSS RULE, CORRESPONDING
C                    TO THE ABSCISSAE XGK(2), XGK(4), ...
C                    WG(1), WG(3), ... ARE SET TO ZERO.
C
      DATA wg(1) / 0.0d0 /
      DATA wg(2) / 0.1294849661 6886969327 0611432679 082d0 /
      DATA wg(3) / 0.0d0 /
      DATA wg(4) / 0.2797053914 8927666790 1467771423 780d0 /
      DATA wg(5) / 0.0d0 /
      DATA wg(6) / 0.3818300505 0511894495 0369775488 975d0 /
      DATA wg(7) / 0.0d0 /
      DATA wg(8) / 0.4179591836 7346938775 5102040816 327d0 /
C
      DATA xgk(1) / 0.9914553711 2081263920 6854697526 329d0 /
      DATA xgk(2) / 0.9491079123 4275852452 6189684047 851d0 /
      DATA xgk(3) / 0.8648644233 5976907278 9712788640 926d0 /
      DATA xgk(4) / 0.7415311855 9939443986 3864773280 788d0 /
      DATA xgk(5) / 0.5860872354 6769113029 4144838258 730d0 /
      DATA xgk(6) / 0.4058451513 7739716690 6606412076 961d0 /
      DATA xgk(7) / 0.2077849550 0789846760 0689403773 245d0 /
      DATA xgk(8) / 0.0000000000 0000000000 0000000000 000d0 /
C
      DATA wgk(1) / 0.0229353220 1052922496 3732008058 970d0 /
      DATA wgk(2) / 0.0630920926 2997855329 0700663189 204d0 /
      DATA wgk(3) / 0.1047900103 2225018383 9876322541 518d0 /
      DATA wgk(4) / 0.1406532597 1552591874 5189590510 238d0 /
      DATA wgk(5) / 0.1690047266 3926790282 6583426598 550d0 /
      DATA wgk(6) / 0.1903505780 6478540991 3256402421 014d0 /
      DATA wgk(7) / 0.2044329400 7529889241 4161999234 649d0 /
      DATA wgk(8) / 0.2094821410 8472782801 2999174891 714d0 /
C
C
C           LIST OF MAJOR VARIABLES
C           -----------------------
C
C           CENTR  - MID POINT OF THE INTERVAL
C           HLGTH  - HALF-LENGTH OF THE INTERVAL
C           ABSC*  - ABSCISSA
C           TABSC* - TRANSFORMED ABSCISSA
C           FVAL*  - FUNCTION VALUE
C           RESG   - RESULT OF THE 7-POINT GAUSS FORMULA
C           RESK   - RESULT OF THE 15-POINT KRONROD FORMULA
C           RESKH  - APPROXIMATION TO THE MEAN VALUE OF THE TRANSFORMED
C                    INTEGRAND OVER (A,B), I.E. TO I/(B-A)
C
C           MACHINE DEPENDENT CONSTANTS
C           ---------------------------
C
C           EPMACH IS THE LARGEST RELATIVE SPACING.
C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
C
C***FIRST EXECUTABLE STATEMENT  DQK15I
      epmach = d1mach(4)
      uflow = d1mach(1)
      dinf = min0(1,inf)
C
      centr = 0.5d+00*(a+b)
      hlgth = 0.5d+00*(b-a)
      tabsc1 = boun+dinf*(0.1d+01-centr)/centr
      ierr = 0
      CALL f(tabsc1,ierr,fval1)
      IF (ierr .LT. 0) RETURN
      IF(inf.EQ.2) THEN
        CALL f(-tabsc1,ierr,fvalt)
        IF (ierr .LT. 0) RETURN
        fval1 = fval1+fvalt
      ENDIF
      fc = (fval1/centr)/centr
C
C           COMPUTE THE 15-POINT KRONROD APPROXIMATION TO
C           THE INTEGRAL, AND ESTIMATE THE ERROR.
C
      resg = wg(8)*fc
      resk = wgk(8)*fc
      resabs = dabs(resk)
      DO 10 j=1,7
        absc = hlgth*xgk(j)
        absc1 = centr-absc
        absc2 = centr+absc
        tabsc1 = boun+dinf*(0.1d+01-absc1)/absc1
        tabsc2 = boun+dinf*(0.1d+01-absc2)/absc2
        CALL f(tabsc1,ierr,fval1)
        IF (ierr .LT. 0) RETURN
        CALL f(tabsc2,ierr,fval2)
        IF (ierr .LT. 0) RETURN
        IF(inf.EQ.2) THEN
          CALL f(-tabsc1,ierr,fvalt)
          IF (ierr .LT. 0) RETURN
          fval1 = fval1+fvalt
        ENDIF
        IF(inf.EQ.2) THEN
          CALL f(-tabsc2,ierr,fvalt)
          IF (ierr .LT. 0) RETURN
          fval2 = fval2+fvalt
        ENDIF
        fval1 = (fval1/absc1)/absc1
        fval2 = (fval2/absc2)/absc2
        fv1(j) = fval1
        fv2(j) = fval2
        fsum = fval1+fval2
        resg = resg+wg(j)*fsum
        resk = resk+wgk(j)*fsum
        resabs = resabs+wgk(j)*(dabs(fval1)+dabs(fval2))
   10 CONTINUE
      reskh = resk*0.5d+00
      resasc = wgk(8)*dabs(fc-reskh)
      DO 20 j=1,7
        resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
   20 CONTINUE
      result = resk*hlgth
      resasc = resasc*hlgth
      resabs = resabs*hlgth
      abserr = dabs((resk-resg)*hlgth)
      IF(resasc.NE.0.0d+00.AND.abserr.NE.0.d0) abserr = resasc*
     * dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
      IF(resabs.GT.uflow/(0.5d+02*epmach)) abserr = dmax1
     * ((epmach*0.5d+02)*resabs,abserr)
      RETURN
      END