← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/interpolate/fitpack/concon.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/interpolate/fitpack/concon.f
1
 
      subroutine concon(iopt,m,x,y,w,v,s,nest,maxtr,maxbin,n,t,c,sq,
2
 
     * sx,bind,wrk,lwrk,iwrk,kwrk,ier)
3
 
c  given the set of data points (x(i),y(i)) and the set of positive
4
 
c  numbers w(i), i=1,2,...,m,subroutine concon determines a cubic spline
5
 
c  approximation s(x) which satisfies the following local convexity
6
 
c  constraints  s''(x(i))*v(i) <= 0, i=1,2,...,m.
7
 
c  the number of knots n and the position t(j),j=1,2,...n is chosen
8
 
c  automatically by the routine in a way that
9
 
c       sq = sum((w(i)*(y(i)-s(x(i))))**2) be <= s.
10
 
c  the fit is given in the b-spline representation (b-spline coef-
11
 
c  ficients c(j),j=1,2,...n-4) and can be evaluated by means of
12
 
c  subroutine splev.
13
 
c
14
 
c  calling sequence:
15
 
c
16
 
c     call concon(iopt,m,x,y,w,v,s,nest,maxtr,maxbin,n,t,c,sq,
17
 
c    * sx,bind,wrk,lwrk,iwrk,kwrk,ier)
18
 
c
19
 
c  parameters:
20
 
c    iopt: integer flag.
21
 
c          if iopt=0, the routine will start with the minimal number of
22
 
c          knots to guarantee that the convexity conditions will be
23
 
c          satisfied. if iopt=1, the routine will continue with the set
24
 
c          of knots found at the last call of the routine.
25
 
c          attention: a call with iopt=1 must always be immediately
26
 
c          preceded by another call with iopt=1 or iopt=0.
27
 
c          unchanged on exit.
28
 
c    m   : integer. on entry m must specify the number of data points.
29
 
c          m > 3. unchanged on exit.
30
 
c    x   : real array of dimension at least (m). before entry, x(i)
31
 
c          must be set to the i-th value of the independent variable x,
32
 
c          for i=1,2,...,m. these values must be supplied in strictly
33
 
c          ascending order. unchanged on exit.
34
 
c    y   : real array of dimension at least (m). before entry, y(i)
35
 
c          must be set to the i-th value of the dependent variable y,
36
 
c          for i=1,2,...,m. unchanged on exit.
37
 
c    w   : real array of dimension at least (m). before entry, w(i)
38
 
c          must be set to the i-th value in the set of weights. the
39
 
c          w(i) must be strictly positive. unchanged on exit.
40
 
c    v   : real array of dimension at least (m). before entry, v(i)
41
 
c          must be set to 1 if s(x) must be locally concave at x(i),
42
 
c          to (-1) if s(x) must be locally convex at x(i) and to 0
43
 
c          if no convexity constraint is imposed at x(i).
44
 
c    s   : real. on entry s must specify an over-estimate for the
45
 
c          the weighted sum of squared residuals sq of the requested
46
 
c          spline. s >=0. unchanged on exit.
47
 
c   nest : integer. on entry nest must contain an over-estimate of the
48
 
c          total number of knots of the spline returned, to indicate
49
 
c          the storage space available to the routine. nest >=8.
50
 
c          in most practical situation nest=m/2 will be sufficient.
51
 
c          always large enough is  nest=m+4. unchanged on exit.
52
 
c  maxtr : integer. on entry maxtr must contain an over-estimate of the
53
 
c          total number of records in the used tree structure, to indic-
54
 
c          ate the storage space available to the routine. maxtr >=1
55
 
c          in most practical situation maxtr=100 will be sufficient.
56
 
c          always large enough is
57
 
c                         nest-5      nest-6
58
 
c              maxtr =  (       ) + (        )  with l the greatest
59
 
c                           l          l+1
60
 
c          integer <= (nest-6)/2 . unchanged on exit.
61
 
c  maxbin: integer. on entry maxbin must contain an over-estimate of the
62
 
c          number of knots where s(x) will have a zero second derivative
63
 
c          maxbin >=1. in most practical situation maxbin = 10 will be
64
 
c          sufficient. always large enough is maxbin=nest-6.
65
 
c          unchanged on exit.
66
 
c    n   : integer.
67
 
c          on exit with ier <=0, n will contain the total number of
68
 
c          knots of the spline approximation returned. if the comput-
69
 
c          ation mode iopt=1 is used this value of n should be left
70
 
c          unchanged between subsequent calls.
71
 
c    t   : real array of dimension at least (nest).
72
 
c          on exit with ier<=0, this array will contain the knots of the
73
 
c          spline,i.e. the position of the interior knots t(5),t(6),...,
74
 
c          t(n-4) as well as the position of the additional knots
75
 
c          t(1)=t(2)=t(3)=t(4)=x(1) and t(n-3)=t(n-2)=t(n-1)=t(n)=x(m)
76
 
c          needed for the the b-spline representation.
77
 
c          if the computation mode iopt=1 is used, the values of t(1),
78
 
c          t(2),...,t(n) should be left unchanged between subsequent
79
 
c          calls.
80
 
c    c   : real array of dimension at least (nest).
81
 
c          on succesful exit, this array will contain the coefficients
82
 
c          c(1),c(2),..,c(n-4) in the b-spline representation of s(x)
83
 
c    sq  : real. unless ier>0 , sq contains the weighted sum of
84
 
c          squared residuals of the spline approximation returned.
85
 
c    sx  : real array of dimension at least m. on exit with ier<=0
86
 
c          this array will contain the spline values s(x(i)),i=1,...,m
87
 
c          if the computation mode iopt=1 is used, the values of sx(1),
88
 
c          sx(2),...,sx(m) should be left unchanged between subsequent
89
 
c          calls.
90
 
c    bind: logical array of dimension at least nest. on exit with ier<=0
91
 
c          this array will indicate the knots where s''(x)=0, i.e.
92
 
c                s''(t(j+3)) .eq. 0 if  bind(j) = .true.
93
 
c                s''(t(j+3)) .ne. 0 if  bind(j) = .false., j=1,2,...,n-6
94
 
c          if the computation mode iopt=1 is used, the values of bind(1)
95
 
c          ,...,bind(n-6) should be left unchanged between subsequent
96
 
c          calls.
97
 
c   wrk  : real array of dimension at least (m*4+nest*8+maxbin*(maxbin+
98
 
c          nest+1)). used as working space.
99
 
c   lwrk : integer. on entry,lwrk must specify the actual dimension of
100
 
c          the array wrk as declared in the calling (sub)program.lwrk
101
 
c          must not be too small (see wrk). unchanged on exit.
102
 
c   iwrk : integer array of dimension at least (maxtr*4+2*(maxbin+1))
103
 
c          used as working space.
104
 
c   kwrk : integer. on entry,kwrk must specify the actual dimension of
105
 
c          the array iwrk as declared in the calling (sub)program. kwrk
106
 
c          must not be too small (see iwrk). unchanged on exit.
107
 
c   ier   : integer. error flag
108
 
c      ier=0 : normal return, s(x) satisfies the concavity/convexity
109
 
c              constraints and sq <= s.
110
 
c      ier<0 : abnormal termination: s(x) satisfies the concavity/
111
 
c              convexity constraints but sq > s.
112
 
c        ier=-3 : the requested storage space exceeds the available
113
 
c                 storage space as specified by the parameter nest.
114
 
c                 probably causes: nest too small. if nest is already
115
 
c                 large (say nest > m/2), it may also indicate that s
116
 
c                 is too small.
117
 
c                 the approximation returned is the least-squares cubic
118
 
c                 spline according to the knots t(1),...,t(n) (n=nest)
119
 
c                 which satisfies the convexity constraints.
120
 
c        ier=-2 : the maximal number of knots n=m+4 has been reached.
121
 
c                 probably causes: s too small.
122
 
c        ier=-1 : the number of knots n is less than the maximal number
123
 
c                 m+4 but concon finds that adding one or more knots
124
 
c                 will not further reduce the value of sq.
125
 
c                 probably causes : s too small.
126
 
c      ier>0 : abnormal termination: no approximation is returned
127
 
c        ier=1  : the number of knots where s''(x)=0 exceeds maxbin.
128
 
c                 probably causes : maxbin too small.
129
 
c        ier=2  : the number of records in the tree structure exceeds
130
 
c                 maxtr.
131
 
c                 probably causes : maxtr too small.
132
 
c        ier=3  : the algoritm finds no solution to the posed quadratic
133
 
c                 programming problem.
134
 
c                 probably causes : rounding errors.
135
 
c        ier=4  : the minimum number of knots (given by n) to guarantee
136
 
c                 that the concavity/convexity conditions will be
137
 
c                 satisfied is greater than nest.
138
 
c                 probably causes: nest too small.
139
 
c        ier=5  : the minimum number of knots (given by n) to guarantee
140
 
c                 that the concavity/convexity conditions will be
141
 
c                 satisfied is greater than m+4.
142
 
c                 probably causes: strongly alternating convexity and
143
 
c                 concavity conditions. normally the situation can be
144
 
c                 coped with by adding n-m-4 extra data points (found
145
 
c                 by linear interpolation e.g.) with a small weight w(i)
146
 
c                 and a v(i) number equal to zero.
147
 
c        ier=10 : on entry, the input data are controlled on validity.
148
 
c                 the following restrictions must be satisfied
149
 
c                   0<=iopt<=1, m>3, nest>=8, s>=0, maxtr>=1, maxbin>=1,
150
 
c                   kwrk>=maxtr*4+2*(maxbin+1), w(i)>0, x(i) < x(i+1),
151
 
c                   lwrk>=m*4+nest*8+maxbin*(maxbin+nest+1)
152
 
c                 if one of these restrictions is found to be violated
153
 
c                 control is immediately repassed to the calling program
154
 
c
155
 
c  further comments:
156
 
c    as an example of the use of the computation mode iopt=1, the
157
 
c    following program segment will cause concon to return control
158
 
c    each time a spline with a new set of knots has been computed.
159
 
c     .............
160
 
c     iopt = 0
161
 
c     s = 0.1e+60  (s very large)
162
 
c     do 10 i=1,m
163
 
c       call concon(iopt,m,x,y,w,v,s,nest,maxtr,maxbin,n,t,c,sq,sx,
164
 
c    *  bind,wrk,lwrk,iwrk,kwrk,ier)
165
 
c       ......
166
 
c       s = sq
167
 
c       iopt=1
168
 
c 10  continue
169
 
c     .............
170
 
c
171
 
c  other subroutines required:
172
 
c    fpcoco,fpcosp,fpbspl,fpadno,fpdeno,fpseno,fpfrno
173
 
c
174
 
c  references:
175
 
c   dierckx p. : an algorithm for cubic spline fitting with convexity
176
 
c                constraints, computing 24 (1980) 349-371.
177
 
c   dierckx p. : an algorithm for least-squares cubic spline fitting
178
 
c                with convexity and concavity constraints, report tw39,
179
 
c                dept. computer science, k.u.leuven, 1978.
180
 
c   dierckx p. : curve and surface fitting with splines, monographs on
181
 
c                numerical analysis, oxford university press, 1993.
182
 
c
183
 
c  author:
184
 
c   p. dierckx
185
 
c   dept. computer science, k.u.leuven
186
 
c   celestijnenlaan 200a, b-3001 heverlee, belgium.
187
 
c   e-mail : Paul.Dierckx@cs.kuleuven.ac.be
188
 
c
189
 
c  creation date : march 1978
190
 
c  latest update : march 1987.
191
 
c
192
 
c  ..
193
 
c  ..scalar arguments..
194
 
      real*8 s,sq
195
 
      integer iopt,m,nest,maxtr,maxbin,n,lwrk,kwrk,ier
196
 
c  ..array arguments..
197
 
      real*8 x(m),y(m),w(m),v(m),t(nest),c(nest),sx(m),wrk(lwrk)
198
 
      integer iwrk(kwrk)
199
 
      logical bind(nest)
200
 
c  ..local scalars..
201
 
      integer i,lwest,kwest,ie,iw,lww
202
 
      real*8 one
203
 
c  ..
204
 
c  set constant
205
 
      one = 0.1e+01
206
 
c  before starting computations a data check is made. if the input data
207
 
c  are invalid, control is immediately repassed to the calling program.
208
 
      ier = 10
209
 
      if(iopt.lt.0 .or. iopt.gt.1) go to 30
210
 
      if(m.lt.4 .or. nest.lt.8) go to 30
211
 
      if(s.lt.0.) go to 30
212
 
      if(maxtr.lt.1 .or. maxbin.lt.1) go to 30
213
 
      lwest = 8*nest+m*4+maxbin*(1+nest+maxbin)
214
 
      kwest = 4*maxtr+2*(maxbin+1)
215
 
      if(lwrk.lt.lwest .or. kwrk.lt.kwest) go to 30
216
 
      if(iopt.gt.0) go to 20
217
 
      if(w(1).le.0.) go to 30
218
 
      if(v(1).gt.0.) v(1) = one
219
 
      if(v(1).lt.0.) v(1) = -one
220
 
      do 10 i=2,m
221
 
         if(x(i-1).ge.x(i) .or. w(i).le.0.) go to 30
222
 
         if(v(i).gt.0.) v(i) = one
223
 
         if(v(i).lt.0.) v(i) = -one
224
 
  10  continue
225
 
  20  ier = 0
226
 
c  we partition the working space and determine the spline approximation
227
 
      ie = 1
228
 
      iw = ie+nest
229
 
      lww = lwrk-nest
230
 
      call fpcoco(iopt,m,x,y,w,v,s,nest,maxtr,maxbin,n,t,c,sq,sx,
231
 
     * bind,wrk(ie),wrk(iw),lww,iwrk,kwrk,ier)
232
 
  30  return
233
 
      end