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

« back to all changes in this revision

Viewing changes to Lib/interpolate/fitpack/fpknot.f

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
ImportĀ upstreamĀ versionĀ 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
      subroutine fpknot(x,m,t,n,fpint,nrdata,nrint,nest,istart)
 
2
c  subroutine fpknot locates an additional knot for a spline of degree
 
3
c  k and adjusts the corresponding parameters,i.e.
 
4
c    t     : the position of the knots.
 
5
c    n     : the number of knots.
 
6
c    nrint : the number of knotintervals.
 
7
c    fpint : the sum of squares of residual right hand sides
 
8
c            for each knot interval.
 
9
c    nrdata: the number of data points inside each knot interval.
 
10
c  istart indicates that the smallest data point at which the new knot
 
11
c  may be added is x(istart+1)
 
12
c  ..
 
13
c  ..scalar arguments..
 
14
      integer m,n,nrint,nest,istart
 
15
c  ..array arguments..
 
16
      real*8 x(m),t(nest),fpint(nest)
 
17
      integer nrdata(nest)
 
18
c  ..local scalars..
 
19
      real*8 an,am,fpmax
 
20
      integer ihalf,j,jbegin,jj,jk,jpoint,k,maxbeg,maxpt,
 
21
     * next,nrx,number
 
22
c  ..
 
23
      k = (n-nrint-1)/2
 
24
c  search for knot interval t(number+k) <= x <= t(number+k+1) where
 
25
c  fpint(number) is maximal on the condition that nrdata(number)
 
26
c  not equals zero.
 
27
      fpmax = 0.
 
28
      jbegin = istart
 
29
      do 20 j=1,nrint
 
30
        jpoint = nrdata(j)
 
31
        if(fpmax.ge.fpint(j) .or. jpoint.eq.0) go to 10
 
32
        fpmax = fpint(j)
 
33
        number = j
 
34
        maxpt = jpoint
 
35
        maxbeg = jbegin
 
36
  10    jbegin = jbegin+jpoint+1
 
37
  20  continue
 
38
c  let coincide the new knot t(number+k+1) with a data point x(nrx)
 
39
c  inside the old knot interval t(number+k) <= x <= t(number+k+1).
 
40
      ihalf = maxpt/2+1
 
41
      nrx = maxbeg+ihalf
 
42
      next = number+1
 
43
      if(next.gt.nrint) go to 40
 
44
c  adjust the different parameters.
 
45
      do 30 j=next,nrint
 
46
        jj = next+nrint-j
 
47
        fpint(jj+1) = fpint(jj)
 
48
        nrdata(jj+1) = nrdata(jj)
 
49
        jk = jj+k
 
50
        t(jk+1) = t(jk)
 
51
  30  continue
 
52
  40  nrdata(number) = ihalf-1
 
53
      nrdata(next) = maxpt-ihalf
 
54
      am = maxpt
 
55
      an = nrdata(number)
 
56
      fpint(number) = fpmax*an/am
 
57
      an = nrdata(next)
 
58
      fpint(next) = fpmax*an/am
 
59
      jk = next+k
 
60
      t(jk) = x(nrx)
 
61
      n = n+1
 
62
      nrint = nrint+1
 
63
      return
 
64
      end