← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/integrate/quadpack/dqk15w.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/dqk15w.f
1
 
      subroutine dqk15w(f,w,p1,p2,p3,p4,kp,a,b,result,abserr,
2
 
     *   resabs,resasc)
3
 
c***begin prologue  dqk15w
4
 
c***date written   810101   (yymmdd)
5
 
c***revision date  830518   (mmddyy)
6
 
c***category no.  h2a2a2
7
 
c***keywords  15-point gauss-kronrod rules
8
 
c***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
9
 
c           de doncker,elise,appl. math. & progr. div. - k.u.leuven
10
 
c***purpose  to compute i = integral of f*w over (a,b), with error
11
 
c                           estimate
12
 
c                       j = integral of abs(f*w) over (a,b)
13
 
c***description
14
 
c
15
 
c           integration rules
16
 
c           standard fortran subroutine
17
 
c           double precision version
18
 
c
19
 
c           parameters
20
 
c             on entry
21
 
c              f      - double precision
22
 
c                       function subprogram defining the integrand
23
 
c                       function f(x). the actual name for f needs to be
24
 
c                       declared e x t e r n a l in the driver program.
25
 
c
26
 
c              w      - double precision
27
 
c                       function subprogram defining the integrand
28
 
c                       weight function w(x). the actual name for w
29
 
c                       needs to be declared e x t e r n a l in the
30
 
c                       calling program.
31
 
c
32
 
c              p1, p2, p3, p4 - double precision
33
 
c                       parameters in the weight function
34
 
c
35
 
c              kp     - integer
36
 
c                       key for indicating the type of weight function
37
 
c
38
 
c              a      - double precision
39
 
c                       lower limit of integration
40
 
c
41
 
c              b      - double precision
42
 
c                       upper limit of integration
43
 
c
44
 
c            on return
45
 
c              result - double precision
46
 
c                       approximation to the integral i
47
 
c                       result is computed by applying the 15-point
48
 
c                       kronrod rule (resk) obtained by optimal addition
49
 
c                       of abscissae to the 7-point gauss rule (resg).
50
 
c
51
 
c              abserr - double precision
52
 
c                       estimate of the modulus of the absolute error,
53
 
c                       which should equal or exceed abs(i-result)
54
 
c
55
 
c              resabs - double precision
56
 
c                       approximation to the integral of abs(f)
57
 
c
58
 
c              resasc - double precision
59
 
c                       approximation to the integral of abs(f-i/(b-a))
60
 
c
61
 
c
62
 
c***references  (none)
63
 
c***routines called  d1mach
64
 
c***end prologue  dqk15w
65
 
c
66
 
      double precision a,absc,absc1,absc2,abserr,b,centr,dabs,dhlgth,
67
 
     *  dmax1,dmin1,d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,
68
 
     *  p1,p2,p3,p4,resabs,resasc,resg,resk,reskh,result,uflow,w,wg,wgk,
69
 
     *  xgk
70
 
      integer j,jtw,jtwm1,kp
71
 
      external f,w
72
 
c
73
 
      dimension fv1(7),fv2(7),xgk(8),wgk(8),wg(4)
74
 
c
75
 
c           the abscissae and weights are given for the interval (-1,1).
76
 
c           because of symmetry only the positive abscissae and their
77
 
c           corresponding weights are given.
78
 
c
79
 
c           xgk    - abscissae of the 15-point gauss-kronrod rule
80
 
c                    xgk(2), xgk(4), ... abscissae of the 7-point
81
 
c                    gauss rule
82
 
c                    xgk(1), xgk(3), ... abscissae which are optimally
83
 
c                    added to the 7-point gauss rule
84
 
c
85
 
c           wgk    - weights of the 15-point gauss-kronrod rule
86
 
c
87
 
c           wg     - weights of the 7-point gauss rule
88
 
c
89
 
      data xgk(1),xgk(2),xgk(3),xgk(4),xgk(5),xgk(6),xgk(7),xgk(8)/
90
 
     *     0.9914553711208126d+00,     0.9491079123427585d+00,
91
 
     *     0.8648644233597691d+00,     0.7415311855993944d+00,
92
 
     *     0.5860872354676911d+00,     0.4058451513773972d+00,
93
 
     *     0.2077849550078985d+00,     0.0000000000000000d+00/
94
 
c
95
 
      data wgk(1),wgk(2),wgk(3),wgk(4),wgk(5),wgk(6),wgk(7),wgk(8)/
96
 
     *     0.2293532201052922d-01,     0.6309209262997855d-01,
97
 
     *     0.1047900103222502d+00,     0.1406532597155259d+00,
98
 
     *     0.1690047266392679d+00,     0.1903505780647854d+00,
99
 
     *     0.2044329400752989d+00,     0.2094821410847278d+00/
100
 
c
101
 
      data wg(1),wg(2),wg(3),wg(4)/
102
 
     *     0.1294849661688697d+00,    0.2797053914892767d+00,
103
 
     *     0.3818300505051889d+00,    0.4179591836734694d+00/
104
 
c
105
 
c
106
 
c           list of major variables
107
 
c           -----------------------
108
 
c
109
 
c           centr  - mid point of the interval
110
 
c           hlgth  - half-length of the interval
111
 
c           absc*  - abscissa
112
 
c           fval*  - function value
113
 
c           resg   - result of the 7-point gauss formula
114
 
c           resk   - result of the 15-point kronrod formula
115
 
c           reskh  - approximation to the mean value of f*w over (a,b),
116
 
c                    i.e. to i/(b-a)
117
 
c
118
 
c           machine dependent constants
119
 
c           ---------------------------
120
 
c
121
 
c           epmach is the largest relative spacing.
122
 
c           uflow is the smallest positive magnitude.
123
 
c
124
 
c***first executable statement  dqk15w
125
 
      epmach = d1mach(4)
126
 
      uflow = d1mach(1)
127
 
c
128
 
      centr = 0.5d+00*(a+b)
129
 
      hlgth = 0.5d+00*(b-a)
130
 
      dhlgth = dabs(hlgth)
131
 
c
132
 
c           compute the 15-point kronrod approximation to the
133
 
c           integral, and estimate the error.
134
 
c
135
 
      fc = f(centr)*w(centr,p1,p2,p3,p4,kp)
136
 
      resg = wg(4)*fc
137
 
      resk = wgk(8)*fc
138
 
      resabs = dabs(resk)
139
 
      do 10 j=1,3
140
 
        jtw = j*2
141
 
        absc = hlgth*xgk(jtw)
142
 
        absc1 = centr-absc
143
 
        absc2 = centr+absc
144
 
        fval1 = f(absc1)*w(absc1,p1,p2,p3,p4,kp)
145
 
        fval2 = f(absc2)*w(absc2,p1,p2,p3,p4,kp)
146
 
        fv1(jtw) = fval1
147
 
        fv2(jtw) = fval2
148
 
        fsum = fval1+fval2
149
 
        resg = resg+wg(j)*fsum
150
 
        resk = resk+wgk(jtw)*fsum
151
 
        resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
152
 
   10 continue
153
 
      do 15 j=1,4
154
 
        jtwm1 = j*2-1
155
 
        absc = hlgth*xgk(jtwm1)
156
 
        absc1 = centr-absc
157
 
        absc2 = centr+absc
158
 
        fval1 = f(absc1)*w(absc1,p1,p2,p3,p4,kp)
159
 
        fval2 = f(absc2)*w(absc2,p1,p2,p3,p4,kp)
160
 
        fv1(jtwm1) = fval1
161
 
        fv2(jtwm1) = fval2
162
 
        fsum = fval1+fval2
163
 
        resk = resk+wgk(jtwm1)*fsum
164
 
        resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
165
 
   15 continue
166
 
      reskh = resk*0.5d+00
167
 
      resasc = wgk(8)*dabs(fc-reskh)
168
 
      do 20 j=1,7
169
 
        resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
170
 
   20 continue
171
 
      result = resk*hlgth
172
 
      resabs = resabs*dhlgth
173
 
      resasc = resasc*dhlgth
174
 
      abserr = dabs((resk-resg)*hlgth)
175
 
      if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
176
 
     *  abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
177
 
      if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1((epmach*
178
 
     *  0.5d+02)*resabs,abserr)
179
 
      return
180
 
      end