1      subroutine qk15i(f,boun,inf,a,b,result,abserr,resabs,resasc,ierr)
2c***begin prologue  qk15i
3c***date written   800101   (yymmdd)
4c***revision date  830518   (yymmdd)
5c***category no.  h2a3a2,h2a4a2
6c***keywords  15-point transformed gauss-kronrod rules
7c***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
8c           de doncker,elise,appl. math. & progr. div. - k.u.leuven
9c***purpose  the original (infinite integration range is mapped
10c            onto the interval (0,1) and (a,b) is a part of (0,1).
11c            it is the purpose to compute
12c            i = integral of transformed integrand over (a,b),
13c            j = integral of abs(transformed integrand) over (a,b).
14c***description
15c
16c           integration rule
17c           standard fortran subroutine
18c           real version
19c
20c           parameters
21c            on entry
22c              f      - subroutine f(x,ierr,result) defining the integrand
23c                       function f(x). the actual name for f needs to be
24c                       declared e x t e r n a l in the calling program.
25c
26c              boun   - real
27c                       finite bound of original integration
28c                       range (set to zero if inf = +2)
29c
30c              inf    - integer
31c                       if inf = -1, the original interval is
32c                                   (-infinity,bound),
33c                       if inf = +1, the original interval is
34c                                   (bound,+infinity),
35c                       if inf = +2, the original interval is
36c                                   (-infinity,+infinity) and
37c                       the integral is computed as the sum of two
38c                       integrals, one over (-infinity,0) and one over
39c                       (0,+infinity).
40c
41c              a      - real
42c                       lower limit for integration over subrange
43c                       of (0,1)
44c
45c              b      - real
46c                       upper limit for integration over subrange
47c                       of (0,1)
48c
49c            on return
50c              result - real
51c                       approximation to the integral i
52c                       result is computed by applying the 15-point
53c                       kronrod rule(resk) obtained by optimal addition
54c                       of abscissae to the 7-point gauss rule(resg).
55c
56c              abserr - real
57c                       estimate of the modulus of the absolute error,
58c                       which should equal or exceed abs(i-result)
59c
60c              resabs - real
61c                       approximation to the integral j
62c
63c              resasc - real
64c                       approximation to the integral of
65c                       abs((transformed integrand)-i/(b-a)) over (a,b)
66c
67c***references  (none)
68c***routines called  r1mach
69c***end prologue  qk15i
70c
71      real a,absc,absc1,absc2,abserr,b,boun,centr,
72     *  dinf,r1mach,epmach,fc,fsum,fval1,fval2,fvalt,fv1,
73     *  fv2,hlgth,resabs,resasc,resg,resk,reskh,result,tabsc1,tabsc2,
74     *  uflow,wg,wgk,xgk
75      integer inf,j,min0
76      external f
77c
78      dimension fv1(7),fv2(7),xgk(8),wgk(8),wg(8)
79c
80c           the abscissae and weights are supplied for the interval
81c           (-1,1).  because of symmetry only the positive abscissae and
82c           their corresponding weights are given.
83c
84c           xgk    - abscissae of the 15-point kronrod rule
85c                    xgk(2), xgk(4), ... abscissae of the 7-point
86c                    gauss rule
87c                    xgk(1), xgk(3), ...  abscissae which are optimally
88c                    added to the 7-point gauss rule
89c
90c           wgk    - weights of the 15-point kronrod rule
91c
92c           wg     - weights of the 7-point gauss rule, corresponding
93c                    to the abscissae xgk(2), xgk(4), ...
94c                    wg(1), wg(3), ... are set to zero.
95c
96      data xgk(1),xgk(2),xgk(3),xgk(4),xgk(5),xgk(6),xgk(7),
97     *  xgk(8)/
98     *     0.9914553711208126e+00,     0.9491079123427585e+00,
99     *     0.8648644233597691e+00,     0.7415311855993944e+00,
100     *     0.5860872354676911e+00,     0.4058451513773972e+00,
101     *     0.2077849550078985e+00,     0.0000000000000000e+00/
102c
103      data wgk(1),wgk(2),wgk(3),wgk(4),wgk(5),wgk(6),wgk(7),
104     *  wgk(8)/
105     *     0.2293532201052922e-01,     0.6309209262997855e-01,
106     *     0.1047900103222502e+00,     0.1406532597155259e+00,
107     *     0.1690047266392679e+00,     0.1903505780647854e+00,
108     *     0.2044329400752989e+00,     0.2094821410847278e+00/
109c
110      data wg(1),wg(2),wg(3),wg(4),wg(5),wg(6),wg(7),wg(8)/
111     *     0.0000000000000000e+00,     0.1294849661688697e+00,
112     *     0.0000000000000000e+00,     0.2797053914892767e+00,
113     *     0.0000000000000000e+00,     0.3818300505051189e+00,
114     *     0.0000000000000000e+00,     0.4179591836734694e+00/
115c
116c
117c           list of major variables
118c           -----------------------
119c
120c           centr  - mid point of the interval
121c           hlgth  - half-length of the interval
122c           absc*  - abscissa
123c           tabsc* - transformed abscissa
124c           fval*  - function value
125c           resg   - result of the 7-point gauss formula
126c           resk   - result of the 15-point kronrod formula
127c           reskh  - approximation to the mean value of the transformed
128c                    integrand over (a,b), i.e. to i/(b-a)
129c
130c           machine dependent constants
131c           ---------------------------
132c
133c           epmach is the largest relative spacing.
134c           uflow is the smallest positive magnitude.
135c
136c***first executable statement  qk15i
137      epmach = r1mach(4)
138      uflow = r1mach(1)
139      dinf = min0(1,inf)
140c
141      centr = 0.5e+00*(a+b)
142      hlgth = 0.5e+00*(b-a)
143      tabsc1 = boun+dinf*(0.1e+01-centr)/centr
144      call f(tabsc1, ierr, fval1)
145      if (ierr.lt.0) return
146      if(inf.eq.2) then
147         call f(-tabsc1, ierr, fval1)
148         if (ierr.lt.0) return
149         fval1 = fval1 + fvalt
150      endif
151      fc = (fval1/centr)/centr
152c
153c           compute the 15-point kronrod approximation to
154c           the integral, and estimate the error.
155c
156      resg = wg(8)*fc
157      resk = wgk(8)*fc
158      resabs = abs(resk)
159      do 10 j=1,7
160        absc = hlgth*xgk(j)
161        absc1 = centr-absc
162        absc2 = centr+absc
163        tabsc1 = boun+dinf*(0.1e+01-absc1)/absc1
164        tabsc2 = boun+dinf*(0.1e+01-absc2)/absc2
165        call f(tabsc1, ierr, fval1)
166        if (ierr.lt.0) return
167        call f(tabsc2, ierr, fval2)
168        if (ierr.lt.0) return
169        if(inf.eq.2) then
170           call f(-tabsc1,ierr,fvalt)
171           if (ierr.lt.0) return
172           fval1 = fval1 + fvalt
173        endif
174        if(inf.eq.2) then
175           call f(-tabsc2,ierr,fvalt)
176           if (ierr.lt.0) return
177           fval2 = fval2 + fvalt
178        endif
179        fval1 = (fval1/absc1)/absc1
180        fval2 = (fval2/absc2)/absc2
181        fv1(j) = fval1
182        fv2(j) = fval2
183        fsum = fval1+fval2
184        resg = resg+wg(j)*fsum
185        resk = resk+wgk(j)*fsum
186        resabs = resabs+wgk(j)*(abs(fval1)+abs(fval2))
187   10 continue
188      reskh = resk*0.5e+00
189      resasc = wgk(8)*abs(fc-reskh)
190      do 20 j=1,7
191        resasc = resasc+wgk(j)*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh))
192   20 continue
193      result = resk*hlgth
194      resasc = resasc*hlgth
195      resabs = resabs*hlgth
196      abserr = abs((resk-resg)*hlgth)
197      if(resasc.ne.0.0e+00.and.abserr.ne.0.e0) abserr = resasc*
198     * amin1(0.1e+01,(0.2e+03*abserr/resasc)**1.5e+00)
199      if(resabs.gt.uflow/(0.5e+02*epmach)) abserr = amax1
200     * ((epmach*0.5e+02)*resabs,abserr)
201      return
202      end
203