1*DECK QK41
2      SUBROUTINE QK41 (F, A, B, RESULT, ABSERR, RESABS, RESASC)
3C***BEGIN PROLOGUE  QK41
4C***PURPOSE  To compute I = Integral of F over (A,B), with error
5C                           estimate
6C                       J = Integral of ABS(F) over (A,B)
7C***LIBRARY   SLATEC (QUADPACK)
8C***CATEGORY  H2A1A2
9C***TYPE      SINGLE PRECISION (QK41-S, DQK41-D)
10C***KEYWORDS  41-POINT GAUSS-KRONROD RULES, QUADPACK, QUADRATURE
11C***AUTHOR  Piessens, Robert
12C             Applied Mathematics and Programming Division
13C             K. U. Leuven
14C           de Doncker, Elise
15C             Applied Mathematics and Programming Division
16C             K. U. Leuven
17C***DESCRIPTION
18C
19C           Integration rules
20C           Standard fortran subroutine
21C           Real version
22C
23C           PARAMETERS
24C            ON ENTRY
25C              F      - Real
26C                       Function subprogram defining the integrand
27C                       FUNCTION F(X). The actual name for F needs to be
28C                       declared E X T E R N A L in the calling program.
29C
30C              A      - Real
31C                       Lower limit of integration
32C
33C              B      - Real
34C                       Upper limit of integration
35C
36C            ON RETURN
37C              RESULT - Real
38C                       Approximation to the integral I
39C                       RESULT is computed by applying the 41-POINT
40C                       GAUSS-KRONROD RULE (RESK) obtained by optimal
41C                       addition of abscissae to the 20-POINT GAUSS
42C                       RULE (RESG).
43C
44C              ABSERR - Real
45C                       Estimate of the modulus of the absolute error,
46C                       which should not exceed ABS(I-RESULT)
47C
48C              RESABS - Real
49C                       Approximation to the integral J
50C
51C              RESASC - Real
52C                       Approximation to the integral of ABS(F-I/(B-A))
53C                       over (A,B)
54C
55C***REFERENCES  (NONE)
56C***ROUTINES CALLED  R1MACH
57C***REVISION HISTORY  (YYMMDD)
58C   800101  DATE WRITTEN
59C   890531  Changed all specific intrinsics to generic.  (WRB)
60C   890531  REVISION DATE from Version 3.2
61C   891214  Prologue converted to Version 4.0 format.  (BAB)
62C***END PROLOGUE  QK41
63C
64      REAL A,ABSC,ABSERR,B,CENTR,DHLGTH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,
65     1  FV1,FV2,HLGTH,RESABS,
66     2  RESASC,RESG,RESK,RESKH,RESULT,R1MACH,UFLOW,
67     3  WG,WGK,XGK
68      INTEGER J,JTW,JTWM1
69      EXTERNAL F
70C
71      DIMENSION FV1(20),FV2(20),XGK(21),WGK(21),WG(10)
72C
73C           THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1).
74C           BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR
75C           CORRESPONDING WEIGHTS ARE GIVEN.
76C
77C           XGK    - ABSCISSAE OF THE 41-POINT GAUSS-KRONROD RULE
78C                    XGK(2), XGK(4), ...  ABSCISSAE OF THE 20-POINT
79C                    GAUSS RULE
80C                    XGK(1), XGK(3), ...  ABSCISSAE WHICH ARE OPTIMALLY
81C                    ADDED TO THE 20-POINT GAUSS RULE
82C
83C           WGK    - WEIGHTS OF THE 41-POINT GAUSS-KRONROD RULE
84C
85C           WG     - WEIGHTS OF THE 20-POINT GAUSS RULE
86C
87      SAVE XGK, WGK, WG
88      DATA XGK(1),XGK(2),XGK(3),XGK(4),XGK(5),XGK(6),XGK(7),XGK(8),
89     1  XGK(9),XGK(10),XGK(11),XGK(12),XGK(13),XGK(14),XGK(15),
90     2  XGK(16),XGK(17),XGK(18),XGK(19),XGK(20),XGK(21)/
91     3     0.9988590315882777E+00,   0.9931285991850949E+00,
92     4     0.9815078774502503E+00,   0.9639719272779138E+00,
93     5     0.9408226338317548E+00,   0.9122344282513259E+00,
94     6     0.8782768112522820E+00,   0.8391169718222188E+00,
95     7     0.7950414288375512E+00,   0.7463319064601508E+00,
96     8     0.6932376563347514E+00,   0.6360536807265150E+00,
97     9     0.5751404468197103E+00,   0.5108670019508271E+00,
98     1     0.4435931752387251E+00,   0.3737060887154196E+00,
99     2     0.3016278681149130E+00,   0.2277858511416451E+00,
100     3     0.1526054652409227E+00,   0.7652652113349733E-01,
101     4     0.0E+00               /
102      DATA WGK(1),WGK(2),WGK(3),WGK(4),WGK(5),WGK(6),WGK(7),WGK(8),
103     1  WGK(9),WGK(10),WGK(11),WGK(12),WGK(13),WGK(14),WGK(15),WGK(16),
104     2  WGK(17),WGK(18),WGK(19),WGK(20),WGK(21)/
105     3     0.3073583718520532E-02,   0.8600269855642942E-02,
106     4     0.1462616925697125E-01,   0.2038837346126652E-01,
107     5     0.2588213360495116E-01,   0.3128730677703280E-01,
108     6     0.3660016975820080E-01,   0.4166887332797369E-01,
109     7     0.4643482186749767E-01,   0.5094457392372869E-01,
110     8     0.5519510534828599E-01,   0.5911140088063957E-01,
111     9     0.6265323755478117E-01,   0.6583459713361842E-01,
112     1     0.6864867292852162E-01,   0.7105442355344407E-01,
113     2     0.7303069033278667E-01,   0.7458287540049919E-01,
114     3     0.7570449768455667E-01,   0.7637786767208074E-01,
115     4     0.7660071191799966E-01/
116      DATA WG(1),WG(2),WG(3),WG(4),WG(5),WG(6),WG(7),WG(8),WG(9),WG(10)/
117     1     0.1761400713915212E-01,    0.4060142980038694E-01,
118     2     0.6267204833410906E-01,    0.8327674157670475E-01,
119     3     0.1019301198172404E+00,    0.1181945319615184E+00,
120     4     0.1316886384491766E+00,    0.1420961093183821E+00,
121     5     0.1491729864726037E+00,    0.1527533871307259E+00/
122C
123C
124C           LIST OF MAJOR VARIABLES
125C           -----------------------
126C
127C           CENTR  - MID POINT OF THE INTERVAL
128C           HLGTH  - HALF-LENGTH OF THE INTERVAL
129C           ABSC   - ABSCISSA
130C           FVAL*  - FUNCTION VALUE
131C           RESG   - RESULT OF THE 20-POINT GAUSS FORMULA
132C           RESK   - RESULT OF THE 41-POINT KRONROD FORMULA
133C           RESKH  - APPROXIMATION TO MEAN VALUE OF F OVER (A,B), I.E.
134C                    TO I/(B-A)
135C
136C           MACHINE DEPENDENT CONSTANTS
137C           ---------------------------
138C
139C           EPMACH IS THE LARGEST RELATIVE SPACING.
140C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
141C
142C***FIRST EXECUTABLE STATEMENT  QK41
143      EPMACH = R1MACH(4)
144      UFLOW = R1MACH(1)
145C
146      CENTR = 0.5E+00*(A+B)
147      HLGTH = 0.5E+00*(B-A)
148      DHLGTH = ABS(HLGTH)
149C
150C           COMPUTE THE 41-POINT GAUSS-KRONROD APPROXIMATION TO
151C           THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR.
152C
153      RESG = 0.0E+00
154      FC = F(CENTR)
155      RESK = WGK(21)*FC
156      RESABS = ABS(RESK)
157      DO 10 J=1,10
158        JTW = J*2
159        ABSC = HLGTH*XGK(JTW)
160        FVAL1 = F(CENTR-ABSC)
161        FVAL2 = F(CENTR+ABSC)
162        FV1(JTW) = FVAL1
163        FV2(JTW) = FVAL2
164        FSUM = FVAL1+FVAL2
165        RESG = RESG+WG(J)*FSUM
166        RESK = RESK+WGK(JTW)*FSUM
167        RESABS = RESABS+WGK(JTW)*(ABS(FVAL1)+ABS(FVAL2))
168   10 CONTINUE
169      DO 15 J = 1,10
170        JTWM1 = J*2-1
171        ABSC = HLGTH*XGK(JTWM1)
172        FVAL1 = F(CENTR-ABSC)
173        FVAL2 = F(CENTR+ABSC)
174        FV1(JTWM1) = FVAL1
175        FV2(JTWM1) = FVAL2
176        FSUM = FVAL1+FVAL2
177        RESK = RESK+WGK(JTWM1)*FSUM
178        RESABS = RESABS+WGK(JTWM1)*(ABS(FVAL1)+ABS(FVAL2))
179   15 CONTINUE
180      RESKH = RESK*0.5E+00
181      RESASC = WGK(21)*ABS(FC-RESKH)
182      DO 20 J=1,20
183        RESASC = RESASC+WGK(J)*(ABS(FV1(J)-RESKH)+ABS(FV2(J)-RESKH))
184   20 CONTINUE
185      RESULT = RESK*HLGTH
186      RESABS = RESABS*DHLGTH
187      RESASC = RESASC*DHLGTH
188      ABSERR = ABS((RESK-RESG)*HLGTH)
189      IF(RESASC.NE.0.0E+00.AND.ABSERR.NE.0.E+00)
190     1  ABSERR = RESASC*MIN(0.1E+01,
191     2  (0.2E+03*ABSERR/RESASC)**1.5E+00)
192      IF(RESABS.GT.UFLOW/(0.5E+02*EPMACH)) ABSERR = MAX
193     1  ((EPMACH*0.5E+02)*RESABS,ABSERR)
194      RETURN
195      END
196