1      SUBROUTINE DQK15I(F,BOUN,INF,A,B,RESULT,ABSERR,RESABS,RESASC)
2C***BEGIN PROLOGUE  DQK15I
3C***DATE WRITTEN   800101   (YYMMDD)
4C***REVISION DATE  830518   (YYMMDD)
5C***REVISION HISTORY (YYMMDD)
6C   000601   Changed DMAX1/DMIN1/DABS to generic MAX/MIN/ABS
7C***CATEGORY NO.  H2A3A2,H2A4A2
8C***KEYWORDS  15-POINT TRANSFORMED GAUSS-KRONROD RULES
9C***AUTHOR  PIESSENS, ROBERT, APPLIED MATH. AND PROGR. DIV. -
10C             K. U. LEUVEN
11C           DE DONCKER, ELISE, APPLIED MATH. AND PROGR. DIV. -
12C             K. U. LEUVEN
13C***PURPOSE  The original (infinite integration range is mapped
14C            onto the interval (0,1) and (A,B) is a part of (0,1).
15C            it is the purpose to compute
16C            I = Integral of transformed integrand over (A,B),
17C            J = Integral of ABS(Transformed Integrand) over (A,B).
18C***DESCRIPTION
19C
20C           Integration Rule
21C           Standard Fortran subroutine
22C           Double precision version
23C
24C           PARAMETERS
25C            ON ENTRY
26C              F      - Double precision
27C                       Fuction subprogram defining the integrand
28C                       FUNCTION F(X). The actual name for F needs to be
29C                       Declared E X T E R N A L in the calling program.
30C
31C              BOUN   - Double precision
32C                       Finite bound of original integration
33C                       Range (SET TO ZERO IF INF = +2)
34C
35C              INF    - Integer
36C                       If INF = -1, the original interval is
37C                                   (-INFINITY,BOUND),
38C                       If INF = +1, the original interval is
39C                                   (BOUND,+INFINITY),
40C                       If INF = +2, the original interval is
41C                                   (-INFINITY,+INFINITY) AND
42C                       The integral is computed as the sum of two
43C                       integrals, one over (-INFINITY,0) and one over
44C                       (0,+INFINITY).
45C
46C              A      - Double precision
47C                       Lower limit for integration over subrange
48C                       of (0,1)
49C
50C              B      - Double precision
51C                       Upper limit for integration over subrange
52C                       of (0,1)
53C
54C            ON RETURN
55C              RESULT - Double precision
56C                       Approximation to the integral I
57C                       Result is computed by applying the 15-POINT
58C                       KRONROD RULE(RESK) obtained by optimal addition
59C                       of abscissae to the 7-POINT GAUSS RULE(RESG).
60C
61C              ABSERR - Double precision
62C                       Estimate of the modulus of the absolute error,
63C                       WHICH SHOULD EQUAL or EXCEED ABS(I-RESULT)
64C
65C              RESABS - Double precision
66C                       Approximation to the integral J
67C
68C              RESASC - Double precision
69C                       Approximation to the integral of
70C                       ABS((TRANSFORMED INTEGRAND)-I/(B-A)) over (A,B)
71C***REFERENCES  (NONE)
72C***ROUTINES CALLED  D1MACH
73C***END PROLOGUE  DQK15I
74C
75      DOUBLE PRECISION A,ABSC,ABSC1,ABSC2,ABSERR,B,BOUN,CENTR,DINF,
76     1  D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,
77     2  RESABS,RESASC,RESG,RESK,RESKH,RESULT,TABSC1,TABSC2,UFLOW,WG,WGK,
78     3  XGK
79      INTEGER INF,J
80      EXTERNAL F
81C
82      DIMENSION FV1(7),FV2(7),XGK(8),WGK(8),WG(8)
83C
84C           THE ABSCISSAE AND WEIGHTS ARE SUPPLIED FOR THE INTERVAL
85C           (-1,1).  BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND
86C           THEIR CORRESPONDING WEIGHTS ARE GIVEN.
87C
88C           XGK    - ABSCISSAE OF THE 15-POINT KRONROD RULE
89C                    XGK(2), XGK(4), ... ABSCISSAE OF THE 7-POINT
90C                    GAUSS RULE
91C                    XGK(1), XGK(3), ...  ABSCISSAE WHICH ARE OPTIMALLY
92C                    ADDED TO THE 7-POINT GAUSS RULE
93C
94C           WGK    - WEIGHTS OF THE 15-POINT KRONROD RULE
95C
96C           WG     - WEIGHTS OF THE 7-POINT GAUSS RULE, CORRESPONDING
97C                    TO THE ABSCISSAE XGK(2), XGK(4), ...
98C                    WG(1), WG(3), ... ARE SET TO ZERO.
99C
100      DATA XGK(1),XGK(2),XGK(3),XGK(4),XGK(5),XGK(6),XGK(7),XGK(8)/
101     1     0.9914553711208126D+00,     0.9491079123427585D+00,
102     2     0.8648644233597691D+00,     0.7415311855993944D+00,
103     3     0.5860872354676911D+00,     0.4058451513773972D+00,
104     4     0.2077849550078985D+00,     0.0000000000000000D+00/
105C
106      DATA WGK(1),WGK(2),WGK(3),WGK(4),WGK(5),WGK(6),WGK(7),WGK(8)/
107     1     0.2293532201052922D-01,     0.6309209262997855D-01,
108     2     0.1047900103222502D+00,     0.1406532597155259D+00,
109     3     0.1690047266392679D+00,     0.1903505780647854D+00,
110     4     0.2044329400752989D+00,     0.2094821410847278D+00/
111C
112      DATA WG(1),WG(2),WG(3),WG(4),WG(5),WG(6),WG(7),WG(8)/
113     1     0.0000000000000000D+00,     0.1294849661688697D+00,
114     2     0.0000000000000000D+00,     0.2797053914892767D+00,
115     3     0.0000000000000000D+00,     0.3818300505051189D+00,
116     4     0.0000000000000000D+00,     0.4179591836734694D+00/
117C
118C
119C           LIST OF MAJOR VARIABLES
120C           -----------------------
121C
122C           CENTR  - MID POINT OF THE INTERVAL
123C           HLGTH  - HALF-LENGTH OF THE INTERVAL
124C           ABSC*  - ABSCISSA
125C           TABSC* - TRANSFORMED ABSCISSA
126C           FVAL*  - FUNCTION VALUE
127C           RESG   - RESULT OF THE 7-POINT GAUSS FORMULA
128C           RESK   - RESULT OF THE 15-POINT KRONROD FORMULA
129C           RESKH  - APPROXIMATION TO THE MEAN VALUE OF THE TRANSFORMED
130C                    INTEGRAND OVER (A,B), I.E. TO I/(B-A)
131C
132C           MACHINE DEPENDENT CONSTANTS
133C           ---------------------------
134C
135C           EPMACH IS THE LARGEST RELATIVE SPACING.
136C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
137C
138C***FIRST EXECUTABLE STATEMENT  DQK15I
139      EPMACH = D1MACH(4)
140      UFLOW = D1MACH(1)
141      DINF = MIN0(1,INF)
142C
143      CENTR = 0.5D+00*(A+B)
144      HLGTH = 0.5D+00*(B-A)
145      TABSC1 = BOUN+DINF*(0.1D+01-CENTR)/CENTR
146      FVAL1 = F(TABSC1)
147      IF(INF.EQ.2) FVAL1 = FVAL1+F(-TABSC1)
148      FC = (FVAL1/CENTR)/CENTR
149C
150C           COMPUTE THE 15-POINT KRONROD APPROXIMATION TO
151C           THE INTEGRAL, AND ESTIMATE THE ERROR.
152C
153      RESG = WG(8)*FC
154      RESK = WGK(8)*FC
155      RESABS = ABS(RESK)
156      DO 10 J=1,7
157        ABSC = HLGTH*XGK(J)
158        ABSC1 = CENTR-ABSC
159        ABSC2 = CENTR+ABSC
160        TABSC1 = BOUN+DINF*(0.1D+01-ABSC1)/ABSC1
161        TABSC2 = BOUN+DINF*(0.1D+01-ABSC2)/ABSC2
162        FVAL1 = F(TABSC1)
163        FVAL2 = F(TABSC2)
164        IF(INF.EQ.2) FVAL1 = FVAL1+F(-TABSC1)
165        IF(INF.EQ.2) FVAL2 = FVAL2+F(-TABSC2)
166        FVAL1 = (FVAL1/ABSC1)/ABSC1
167        FVAL2 = (FVAL2/ABSC2)/ABSC2
168        FV1(J) = FVAL1
169        FV2(J) = FVAL2
170        FSUM = FVAL1+FVAL2
171        RESG = RESG+WG(J)*FSUM
172        RESK = RESK+WGK(J)*FSUM
173        RESABS = RESABS+WGK(J)*(ABS(FVAL1)+ABS(FVAL2))
174   10 CONTINUE
175      RESKH = RESK*0.5D+00
176      RESASC = WGK(8)*ABS(FC-RESKH)
177      DO 20 J=1,7
178        RESASC = RESASC+WGK(J)*(ABS(FV1(J)-RESKH)+ABS(FV2(J)-RESKH))
179   20 CONTINUE
180      RESULT = RESK*HLGTH
181      RESASC = RESASC*HLGTH
182      RESABS = RESABS*HLGTH
183      ABSERR = ABS((RESK-RESG)*HLGTH)
184      IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.D0) ABSERR = RESASC*
185     1 MIN(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00)
186      IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = MAX
187     1 ((EPMACH*0.5D+02)*RESABS,ABSERR)
188      RETURN
189      END
190