← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/interpolate/fitpack/cocosp.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/cocosp.f
1
 
      subroutine cocosp(m,x,y,w,n,t,e,maxtr,maxbin,c,sq,sx,bind,wrk,
2
 
     * 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 cocosp determines the weighted
5
 
c  least-squares cubic spline s(x) with given knots t(j),j=1,2,...,n
6
 
c  which satisfies the following concavity/convexity conditions
7
 
c      s''(t(j+3))*e(j) <= 0, j=1,2,...n-6
8
 
c  the fit is given in the b-spline representation( b-spline coef-
9
 
c  ficients c(j),j=1,2,...n-4) and can be evaluated by means of
10
 
c  subroutine splev.
11
 
c
12
 
c  calling sequence:
13
 
c     call cocosp(m,x,y,w,n,t,e,maxtr,maxbin,c,sq,sx,bind,wrk,
14
 
c    * lwrk,iwrk,kwrk,ier)
15
 
c
16
 
c  parameters:
17
 
c    m   : integer. on entry m must specify the number of data points.
18
 
c          m > 3. unchanged on exit.
19
 
c    x   : real array of dimension at least (m). before entry, x(i)
20
 
c          must be set to the i-th value of the independent variable x,
21
 
c          for i=1,2,...,m. these values must be supplied in strictly
22
 
c          ascending order. unchanged on exit.
23
 
c    y   : real array of dimension at least (m). before entry, y(i)
24
 
c          must be set to the i-th value of the dependent variable y,
25
 
c          for i=1,2,...,m. unchanged on exit.
26
 
c    w   : real array of dimension at least (m). before entry, w(i)
27
 
c          must be set to the i-th value in the set of weights. the
28
 
c          w(i) must be strictly positive. unchanged on exit.
29
 
c    n   : integer. on entry n must contain the total number of knots
30
 
c          of the cubic spline. m+4>=n>=8. unchanged on exit.
31
 
c    t   : real array of dimension at least (n). before entry, this
32
 
c          array must contain the knots of the spline, i.e. the position
33
 
c          of the interior knots t(5),t(6),...,t(n-4) as well as the
34
 
c          position of the boundary knots t(1),t(2),t(3),t(4) and t(n-3)
35
 
c          t(n-2),t(n-1),t(n) needed for the b-spline representation.
36
 
c          unchanged on exit. see also the restrictions (ier=10).
37
 
c    e   : real array of dimension at least (n). before entry, e(j)
38
 
c          must be set to 1 if s(x) must be locally concave at t(j+3),
39
 
c          to (-1) if s(x) must be locally convex at t(j+3) and to 0
40
 
c          if no convexity constraint is imposed at t(j+3),j=1,2,..,n-6.
41
 
c          e(n-5),...,e(n) are not used. unchanged on exit.
42
 
c  maxtr : integer. on entry maxtr must contain an over-estimate of the
43
 
c          total number of records in the used tree structure, to indic-
44
 
c          ate the storage space available to the routine. maxtr >=1
45
 
c          in most practical situation maxtr=100 will be sufficient.
46
 
c          always large enough is
47
 
c                         n-5       n-6
48
 
c              maxtr =  (     ) + (     )  with l the greatest
49
 
c                          l        l+1
50
 
c          integer <= (n-6)/2 . unchanged on exit.
51
 
c  maxbin: integer. on entry maxbin must contain an over-estimate of the
52
 
c          number of knots where s(x) will have a zero second derivative
53
 
c          maxbin >=1. in most practical situation maxbin = 10 will be
54
 
c          sufficient. always large enough is maxbin=n-6.
55
 
c          unchanged on exit.
56
 
c    c   : real array of dimension at least (n).
57
 
c          on succesful exit, this array will contain the coefficients
58
 
c          c(1),c(2),..,c(n-4) in the b-spline representation of s(x)
59
 
c    sq  : real. on succesful exit, sq contains the weighted sum of
60
 
c          squared residuals of the spline approximation returned.
61
 
c    sx  : real array of dimension at least m. on succesful exit
62
 
c          this array will contain the spline values s(x(i)),i=1,...,m
63
 
c   bind : logical array of dimension at least (n). on succesful exit
64
 
c          this array will indicate the knots where s''(x)=0, i.e.
65
 
c                s''(t(j+3)) .eq. 0 if  bind(j) = .true.
66
 
c                s''(t(j+3)) .ne. 0 if  bind(j) = .false., j=1,2,...,n-6
67
 
c   wrk  : real array of dimension at least  m*4+n*7+maxbin*(maxbin+n+1)
68
 
c          used as working space.
69
 
c   lwrk : integer. on entry,lwrk must specify the actual dimension of
70
 
c          the array wrk as declared in the calling (sub)program.lwrk
71
 
c          must not be too small (see wrk). unchanged on exit.
72
 
c   iwrk : integer array of dimension at least (maxtr*4+2*(maxbin+1))
73
 
c          used as working space.
74
 
c   kwrk : integer. on entry,kwrk must specify the actual dimension of
75
 
c          the array iwrk as declared in the calling (sub)program. kwrk
76
 
c          must not be too small (see iwrk). unchanged on exit.
77
 
c   ier   : integer. error flag
78
 
c      ier=0 : succesful exit.
79
 
c      ier>0 : abnormal termination: no approximation is returned
80
 
c        ier=1  : the number of knots where s''(x)=0 exceeds maxbin.
81
 
c                 probably causes : maxbin too small.
82
 
c        ier=2  : the number of records in the tree structure exceeds
83
 
c                 maxtr.
84
 
c                 probably causes : maxtr too small.
85
 
c        ier=3  : the algoritm finds no solution to the posed quadratic
86
 
c                 programming problem.
87
 
c                 probably causes : rounding errors.
88
 
c        ier=10 : on entry, the input data are controlled on validity.
89
 
c                 the following restrictions must be satisfied
90
 
c                   m>3, maxtr>=1, maxbin>=1, 8<=n<=m+4,w(i) > 0,
91
 
c                   x(1)<x(2)<...<x(m), t(1)<=t(2)<=t(3)<=t(4)<=x(1),
92
 
c                   x(1)<t(5)<t(6)<...<t(n-4)<x(m)<=t(n-3)<=...<=t(n),
93
 
c                   kwrk>=maxtr*4+2*(maxbin+1),
94
 
c                   lwrk>=m*4+n*7+maxbin*(maxbin+n+1),
95
 
c                   the schoenberg-whitney conditions, i.e. there must
96
 
c                   be a subset of data points xx(j) such that
97
 
c                     t(j) < xx(j) < t(j+4), j=1,2,...,n-4
98
 
c                 if one of these restrictions is found to be violated
99
 
c                 control is immediately repassed to the calling program
100
 
c
101
 
c
102
 
c  other subroutines required:
103
 
c    fpcosp,fpbspl,fpadno,fpdeno,fpseno,fpfrno,fpchec
104
 
c
105
 
c  references:
106
 
c   dierckx p. : an algorithm for cubic spline fitting with convexity
107
 
c                constraints, computing 24 (1980) 349-371.
108
 
c   dierckx p. : an algorithm for least-squares cubic spline fitting
109
 
c                with convexity and concavity constraints, report tw39,
110
 
c                dept. computer science, k.u.leuven, 1978.
111
 
c   dierckx p. : curve and surface fitting with splines, monographs on
112
 
c                numerical analysis, oxford university press, 1993.
113
 
c
114
 
c  author:
115
 
c   p. dierckx
116
 
c   dept. computer science, k.u.leuven
117
 
c   celestijnenlaan 200a, b-3001 heverlee, belgium.
118
 
c   e-mail : Paul.Dierckx@cs.kuleuven.ac.be
119
 
c
120
 
c  creation date : march 1978
121
 
c  latest update : march 1987.
122
 
c
123
 
c  ..
124
 
c  ..scalar arguments..
125
 
      real*8 sq
126
 
      integer m,n,maxtr,maxbin,lwrk,kwrk,ier
127
 
c  ..array arguments..
128
 
      real*8 x(m),y(m),w(m),t(n),e(n),c(n),sx(m),wrk(lwrk)
129
 
      integer iwrk(kwrk)
130
 
      logical bind(n)
131
 
c  ..local scalars..
132
 
      integer i,ia,ib,ic,iq,iu,iz,izz,ji,jib,jjb,jl,jr,ju,kwest,
133
 
     * lwest,mb,nm,n6
134
 
      real*8 one
135
 
c  ..
136
 
c  set constant
137
 
      one = 0.1e+01
138
 
c  before starting computations a data check is made. if the input data
139
 
c  are invalid, control is immediately repassed to the calling program.
140
 
      ier = 10
141
 
      if(m.lt.4 .or. n.lt.8) go to 40
142
 
      if(maxtr.lt.1 .or. maxbin.lt.1) go to 40
143
 
      lwest = 7*n+m*4+maxbin*(1+n+maxbin)
144
 
      kwest = 4*maxtr+2*(maxbin+1)
145
 
      if(lwrk.lt.lwest .or. kwrk.lt.kwest) go to 40
146
 
      if(w(1).le.0.) go to 40
147
 
      do 10 i=2,m
148
 
         if(x(i-1).ge.x(i) .or. w(i).le.0.) go to 40
149
 
  10  continue
150
 
      call fpchec(x,m,t,n,3,ier)
151
 
      if (ier.eq.0) go to 20
152
 
      go to 40
153
 
c  set numbers e(i)
154
 
  20  n6 = n-6
155
 
      do 30 i=1,n6
156
 
        if(e(i).gt.0.) e(i) = one
157
 
        if(e(i).lt.0.) e(i) = -one
158
 
  30  continue
159
 
c  we partition the working space and determine the spline approximation
160
 
      nm = n+maxbin
161
 
      mb = maxbin+1
162
 
      ia = 1
163
 
      ib = ia+4*n
164
 
      ic = ib+nm*maxbin
165
 
      iz = ic+n
166
 
      izz = iz+n
167
 
      iu = izz+n
168
 
      iq = iu+maxbin
169
 
      ji = 1
170
 
      ju = ji+maxtr
171
 
      jl = ju+maxtr
172
 
      jr = jl+maxtr
173
 
      jjb = jr+maxtr
174
 
      jib = jjb+mb
175
 
      call fpcosp(m,x,y,w,n,t,e,maxtr,maxbin,c,sq,sx,bind,nm,mb,wrk(ia),
176
 
     *
177
 
     * wrk(ib),wrk(ic),wrk(iz),wrk(izz),wrk(iu),wrk(iq),iwrk(ji),
178
 
     * iwrk(ju),iwrk(jl),iwrk(jr),iwrk(jjb),iwrk(jib),ier)
179
 
  40  return
180
 
      end