← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/integrate/quadpack/dqk21.f

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/integrate/quadpack/dqk21.f
1
 
      subroutine dqk21(f,a,b,result,abserr,resabs,resasc)
2
 
c***begin prologue  dqk21
3
 
c***date written   800101   (yymmdd)
4
 
c***revision date  830518   (yymmdd)
5
 
c***category no.  h2a1a2
6
 
c***keywords  21-point gauss-kronrod rules
7
 
c***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
8
 
c           de doncker,elise,appl. math. & progr. div. - k.u.leuven
9
 
c***purpose  to compute i = integral of f over (a,b), with error
10
 
c                           estimate
11
 
c                       j = integral of abs(f) over (a,b)
12
 
c***description
13
 
c
14
 
c           integration rules
15
 
c           standard fortran subroutine
16
 
c           double precision version
17
 
c
18
 
c           parameters
19
 
c            on entry
20
 
c              f      - double precision
21
 
c                       function subprogram defining the integrand
22
 
c                       function f(x). the actual name for f needs to be
23
 
c                       declared e x t e r n a l in the driver program.
24
 
c
25
 
c              a      - double precision
26
 
c                       lower limit of integration
27
 
c
28
 
c              b      - double precision
29
 
c                       upper limit of integration
30
 
c
31
 
c            on return
32
 
c              result - double precision
33
 
c                       approximation to the integral i
34
 
c                       result is computed by applying the 21-point
35
 
c                       kronrod rule (resk) obtained by optimal addition
36
 
c                       of abscissae to the 10-point gauss rule (resg).
37
 
c
38
 
c              abserr - double precision
39
 
c                       estimate of the modulus of the absolute error,
40
 
c                       which should not exceed abs(i-result)
41
 
c
42
 
c              resabs - double precision
43
 
c                       approximation to the integral j
44
 
c
45
 
c              resasc - double precision
46
 
c                       approximation to the integral of abs(f-i/(b-a))
47
 
c                       over (a,b)
48
 
c
49
 
c***references  (none)
50
 
c***routines called  d1mach
51
 
c***end prologue  dqk21
52
 
c
53
 
      double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
54
 
     *  d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
55
 
     *  resg,resk,reskh,result,uflow,wg,wgk,xgk
56
 
      integer j,jtw,jtwm1
57
 
      external f
58
 
c
59
 
      dimension fv1(10),fv2(10),wg(5),wgk(11),xgk(11)
60
 
c
61
 
c           the abscissae and weights are given for the interval (-1,1).
62
 
c           because of symmetry only the positive abscissae and their
63
 
c           corresponding weights are given.
64
 
c
65
 
c           xgk    - abscissae of the 21-point kronrod rule
66
 
c                    xgk(2), xgk(4), ...  abscissae of the 10-point
67
 
c                    gauss rule
68
 
c                    xgk(1), xgk(3), ...  abscissae which are optimally
69
 
c                    added to the 10-point gauss rule
70
 
c
71
 
c           wgk    - weights of the 21-point kronrod rule
72
 
c
73
 
c           wg     - weights of the 10-point gauss rule
74
 
c
75
 
c
76
 
c gauss quadrature weights and kronron quadrature abscissae and weights
77
 
c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
78
 
c bell labs, nov. 1981.
79
 
c
80
 
      data wg  (  1) / 0.0666713443 0868813759 3568809893 332 d0 /
81
 
      data wg  (  2) / 0.1494513491 5058059314 5776339657 697 d0 /
82
 
      data wg  (  3) / 0.2190863625 1598204399 5534934228 163 d0 /
83
 
      data wg  (  4) / 0.2692667193 0999635509 1226921569 469 d0 /
84
 
      data wg  (  5) / 0.2955242247 1475287017 3892994651 338 d0 /
85
 
c
86
 
      data xgk (  1) / 0.9956571630 2580808073 5527280689 003 d0 /
87
 
      data xgk (  2) / 0.9739065285 1717172007 7964012084 452 d0 /
88
 
      data xgk (  3) / 0.9301574913 5570822600 1207180059 508 d0 /
89
 
      data xgk (  4) / 0.8650633666 8898451073 2096688423 493 d0 /
90
 
      data xgk (  5) / 0.7808177265 8641689706 3717578345 042 d0 /
91
 
      data xgk (  6) / 0.6794095682 9902440623 4327365114 874 d0 /
92
 
      data xgk (  7) / 0.5627571346 6860468333 9000099272 694 d0 /
93
 
      data xgk (  8) / 0.4333953941 2924719079 9265943165 784 d0 /
94
 
      data xgk (  9) / 0.2943928627 0146019813 1126603103 866 d0 /
95
 
      data xgk ( 10) / 0.1488743389 8163121088 4826001129 720 d0 /
96
 
      data xgk ( 11) / 0.0000000000 0000000000 0000000000 000 d0 /
97
 
c
98
 
      data wgk (  1) / 0.0116946388 6737187427 8064396062 192 d0 /
99
 
      data wgk (  2) / 0.0325581623 0796472747 8818972459 390 d0 /
100
 
      data wgk (  3) / 0.0547558965 7435199603 1381300244 580 d0 /
101
 
      data wgk (  4) / 0.0750396748 1091995276 7043140916 190 d0 /
102
 
      data wgk (  5) / 0.0931254545 8369760553 5065465083 366 d0 /
103
 
      data wgk (  6) / 0.1093871588 0229764189 9210590325 805 d0 /
104
 
      data wgk (  7) / 0.1234919762 6206585107 7958109831 074 d0 /
105
 
      data wgk (  8) / 0.1347092173 1147332592 8054001771 707 d0 /
106
 
      data wgk (  9) / 0.1427759385 7706008079 7094273138 717 d0 /
107
 
      data wgk ( 10) / 0.1477391049 0133849137 4841515972 068 d0 /
108
 
      data wgk ( 11) / 0.1494455540 0291690566 4936468389 821 d0 /
109
 
c
110
 
c
111
 
c           list of major variables
112
 
c           -----------------------
113
 
c
114
 
c           centr  - mid point of the interval
115
 
c           hlgth  - half-length of the interval
116
 
c           absc   - abscissa
117
 
c           fval*  - function value
118
 
c           resg   - result of the 10-point gauss formula
119
 
c           resk   - result of the 21-point kronrod formula
120
 
c           reskh  - approximation to the mean value of f over (a,b),
121
 
c                    i.e. to i/(b-a)
122
 
c
123
 
c
124
 
c           machine dependent constants
125
 
c           ---------------------------
126
 
c
127
 
c           epmach is the largest relative spacing.
128
 
c           uflow is the smallest positive magnitude.
129
 
c
130
 
c***first executable statement  dqk21
131
 
      epmach = d1mach(4)
132
 
      uflow = d1mach(1)
133
 
c
134
 
      centr = 0.5d+00*(a+b)
135
 
      hlgth = 0.5d+00*(b-a)
136
 
      dhlgth = dabs(hlgth)
137
 
c
138
 
c           compute the 21-point kronrod approximation to
139
 
c           the integral, and estimate the absolute error.
140
 
c
141
 
      resg = 0.0d+00
142
 
      fc = f(centr)
143
 
      resk = wgk(11)*fc
144
 
      resabs = dabs(resk)
145
 
      do 10 j=1,5
146
 
        jtw = 2*j
147
 
        absc = hlgth*xgk(jtw)
148
 
        fval1 = f(centr-absc)
149
 
        fval2 = f(centr+absc)
150
 
        fv1(jtw) = fval1
151
 
        fv2(jtw) = fval2
152
 
        fsum = fval1+fval2
153
 
        resg = resg+wg(j)*fsum
154
 
        resk = resk+wgk(jtw)*fsum
155
 
        resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
156
 
   10 continue
157
 
      do 15 j = 1,5
158
 
        jtwm1 = 2*j-1
159
 
        absc = hlgth*xgk(jtwm1)
160
 
        fval1 = f(centr-absc)
161
 
        fval2 = f(centr+absc)
162
 
        fv1(jtwm1) = fval1
163
 
        fv2(jtwm1) = fval2
164
 
        fsum = fval1+fval2
165
 
        resk = resk+wgk(jtwm1)*fsum
166
 
        resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
167
 
   15 continue
168
 
      reskh = resk*0.5d+00
169
 
      resasc = wgk(11)*dabs(fc-reskh)
170
 
      do 20 j=1,10
171
 
        resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
172
 
   20 continue
173
 
      result = resk*hlgth
174
 
      resabs = resabs*dhlgth
175
 
      resasc = resasc*dhlgth
176
 
      abserr = dabs((resk-resg)*hlgth)
177
 
      if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
178
 
     *  abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
179
 
      if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
180
 
     *  ((epmach*0.5d+02)*resabs,abserr)
181
 
      return
182
 
      end