← Back to branch summary

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

« back to all changes in this revision

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