← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/integrate/quadpack/dqawc.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/dqawc.f
1
 
      subroutine dqawc(f,a,b,c,epsabs,epsrel,result,abserr,neval,ier,
2
 
     *   limit,lenw,last,iwork,work)
3
 
c***begin prologue  dqawc
4
 
c***date written   800101   (yymmdd)
5
 
c***revision date  830518   (yymmdd)
6
 
c***category no.  h2a2a1,j4
7
 
c***keywords  automatic integrator, special-purpose,
8
 
c             cauchy principal value,
9
 
c             clenshaw-curtis, globally adaptive
10
 
c***author  piessens,robert ,appl. math. & progr. div. - k.u.leuven
11
 
c           de doncker,elise,appl. math. & progr. div. - k.u.leuven
12
 
c***purpose  the routine calculates an approximation result to a
13
 
c            cauchy principal value i = integral of f*w over (a,b)
14
 
c            (w(x) = 1/((x-c), c.ne.a, c.ne.b), hopefully satisfying
15
 
c            following claim for accuracy
16
 
c            abs(i-result).le.max(epsabe,epsrel*abs(i)).
17
 
c***description
18
 
c
19
 
c        computation of a cauchy principal value
20
 
c        standard fortran subroutine
21
 
c        double precision version
22
 
c
23
 
c
24
 
c        parameters
25
 
c         on entry
26
 
c            f      - double precision
27
 
c                     function subprogram defining the integrand
28
 
c                     function f(x). the actual name for f needs to be
29
 
c                     declared e x t e r n a l in the driver program.
30
 
c
31
 
c            a      - double precision
32
 
c                     under limit of integration
33
 
c
34
 
c            b      - double precision
35
 
c                     upper limit of integration
36
 
c
37
 
c            c      - parameter in the weight function, c.ne.a, c.ne.b.
38
 
c                     if c = a or c = b, the routine will end with
39
 
c                     ier = 6 .
40
 
c
41
 
c            epsabs - double precision
42
 
c                     absolute accuracy requested
43
 
c            epsrel - double precision
44
 
c                     relative accuracy requested
45
 
c                     if  epsabs.le.0
46
 
c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
47
 
c                     the routine will end with ier = 6.
48
 
c
49
 
c         on return
50
 
c            result - double precision
51
 
c                     approximation to the integral
52
 
c
53
 
c            abserr - double precision
54
 
c                     estimate or the modulus of the absolute error,
55
 
c                     which should equal or exceed abs(i-result)
56
 
c
57
 
c            neval  - integer
58
 
c                     number of integrand evaluations
59
 
c
60
 
c            ier    - integer
61
 
c                     ier = 0 normal and reliable termination of the
62
 
c                             routine. it is assumed that the requested
63
 
c                             accuracy has been achieved.
64
 
c                     ier.gt.0 abnormal termination of the routine
65
 
c                             the estimates for integral and error are
66
 
c                             less reliable. it is assumed that the
67
 
c                             requested accuracy has not been achieved.
68
 
c            error messages
69
 
c                     ier = 1 maximum number of subdivisions allowed
70
 
c                             has been achieved. one can allow more sub-
71
 
c                             divisions by increasing the value of limit
72
 
c                             (and taking the according dimension
73
 
c                             adjustments into account). however, if
74
 
c                             this yields no improvement it is advised
75
 
c                             to analyze the integrand in order to
76
 
c                             determine the integration difficulties.
77
 
c                             if the position of a local difficulty
78
 
c                             can be determined (e.g. singularity,
79
 
c                             discontinuity within the interval) one
80
 
c                             will probably gain from splitting up the
81
 
c                             interval at this point and calling
82
 
c                             appropriate integrators on the subranges.
83
 
c                         = 2 the occurrence of roundoff error is detec-
84
 
c                             ted, which prevents the requested
85
 
c                             tolerance from being achieved.
86
 
c                         = 3 extremely bad integrand behaviour occurs
87
 
c                             at some points of the integration
88
 
c                             interval.
89
 
c                         = 6 the input is invalid, because
90
 
c                             c = a or c = b or
91
 
c                             (epsabs.le.0 and
92
 
c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
93
 
c                             or limit.lt.1 or lenw.lt.limit*4.
94
 
c                             result, abserr, neval, last are set to
95
 
c                             zero. exept when lenw or limit is invalid,
96
 
c                             iwork(1), work(limit*2+1) and
97
 
c                             work(limit*3+1) are set to zero, work(1)
98
 
c                             is set to a and work(limit+1) to b.
99
 
c
100
 
c         dimensioning parameters
101
 
c            limit - integer
102
 
c                    dimensioning parameter for iwork
103
 
c                    limit determines the maximum number of subintervals
104
 
c                    in the partition of the given integration interval
105
 
c                    (a,b), limit.ge.1.
106
 
c                    if limit.lt.1, the routine will end with ier = 6.
107
 
c
108
 
c           lenw   - integer
109
 
c                    dimensioning parameter for work
110
 
c                    lenw must be at least limit*4.
111
 
c                    if lenw.lt.limit*4, the routine will end with
112
 
c                    ier = 6.
113
 
c
114
 
c            last  - integer
115
 
c                    on return, last equals the number of subintervals
116
 
c                    produced in the subdivision process, which
117
 
c                    determines the number of significant elements
118
 
c                    actually in the work arrays.
119
 
c
120
 
c         work arrays
121
 
c            iwork - integer
122
 
c                    vector of dimension at least limit, the first k
123
 
c                    elements of which contain pointers
124
 
c                    to the error estimates over the subintervals,
125
 
c                    such that work(limit*3+iwork(1)), ... ,
126
 
c                    work(limit*3+iwork(k)) form a decreasing
127
 
c                    sequence, with k = last if last.le.(limit/2+2),
128
 
c                    and k = limit+1-last otherwise
129
 
c
130
 
c            work  - double precision
131
 
c                    vector of dimension at least lenw
132
 
c                    on return
133
 
c                    work(1), ..., work(last) contain the left
134
 
c                     end points of the subintervals in the
135
 
c                     partition of (a,b),
136
 
c                    work(limit+1), ..., work(limit+last) contain
137
 
c                     the right end points,
138
 
c                    work(limit*2+1), ..., work(limit*2+last) contain
139
 
c                     the integral approximations over the subintervals,
140
 
c                    work(limit*3+1), ..., work(limit*3+last)
141
 
c                     contain the error estimates.
142
 
c
143
 
c***references  (none)
144
 
c***routines called  dqawce,xerror
145
 
c***end prologue  dqawc
146
 
c
147
 
      double precision a,abserr,b,c,epsabs,epsrel,f,result,work
148
 
      integer ier,iwork,last,lenw,limit,lvl,l1,l2,l3,neval
149
 
c
150
 
      dimension iwork(limit),work(lenw)
151
 
c
152
 
      external f
153
 
c
154
 
c         check validity of limit and lenw.
155
 
c
156
 
c***first executable statement  dqawc
157
 
      ier = 6
158
 
      neval = 0
159
 
      last = 0
160
 
      result = 0.0d+00
161
 
      abserr = 0.0d+00
162
 
      if(limit.lt.1.or.lenw.lt.limit*4) go to 10
163
 
c
164
 
c         prepare call for dqawce.
165
 
c
166
 
      l1 = limit+1
167
 
      l2 = limit+l1
168
 
      l3 = limit+l2
169
 
      call dqawce(f,a,b,c,epsabs,epsrel,limit,result,abserr,neval,ier,
170
 
     *  work(1),work(l1),work(l2),work(l3),iwork,last)
171
 
c
172
 
c         call error handler if necessary.
173
 
c
174
 
      lvl = 0
175
 
10    if(ier.eq.6) lvl = 1
176
 
      if(ier.ne.0) call xerror('abnormal return from dqawc',26,ier,lvl)
177
 
      return
178
 
      end