[Scipy-svn] r4550 - in branches/Interpolate1D: . fitpack
scipy-svn@scip...
scipy-svn@scip...
Fri Jul 18 14:44:23 CDT 2008
Author: fcady
Date: 2008-07-18 14:44:12 -0500 (Fri, 18 Jul 2008)
New Revision: 4550
Added:
branches/Interpolate1D/fitpack/
branches/Interpolate1D/fitpack/Makefile
branches/Interpolate1D/fitpack/README
branches/Interpolate1D/fitpack/bispev.f
branches/Interpolate1D/fitpack/clocur.f
branches/Interpolate1D/fitpack/cocosp.f
branches/Interpolate1D/fitpack/concon.f
branches/Interpolate1D/fitpack/concur.f
branches/Interpolate1D/fitpack/cualde.f
branches/Interpolate1D/fitpack/curev.f
branches/Interpolate1D/fitpack/curfit.f
branches/Interpolate1D/fitpack/dblint.f
branches/Interpolate1D/fitpack/evapol.f
branches/Interpolate1D/fitpack/fourco.f
branches/Interpolate1D/fitpack/fpader.f
branches/Interpolate1D/fitpack/fpadno.f
branches/Interpolate1D/fitpack/fpadpo.f
branches/Interpolate1D/fitpack/fpback.f
branches/Interpolate1D/fitpack/fpbacp.f
branches/Interpolate1D/fitpack/fpbfout.f
branches/Interpolate1D/fitpack/fpbisp.f
branches/Interpolate1D/fitpack/fpbspl.f
branches/Interpolate1D/fitpack/fpchec.f
branches/Interpolate1D/fitpack/fpched.f
branches/Interpolate1D/fitpack/fpchep.f
branches/Interpolate1D/fitpack/fpclos.f
branches/Interpolate1D/fitpack/fpcoco.f
branches/Interpolate1D/fitpack/fpcons.f
branches/Interpolate1D/fitpack/fpcosp.f
branches/Interpolate1D/fitpack/fpcsin.f
branches/Interpolate1D/fitpack/fpcurf.f
branches/Interpolate1D/fitpack/fpcuro.f
branches/Interpolate1D/fitpack/fpcyt1.f
branches/Interpolate1D/fitpack/fpcyt2.f
branches/Interpolate1D/fitpack/fpdeno.f
branches/Interpolate1D/fitpack/fpdisc.f
branches/Interpolate1D/fitpack/fpfrno.f
branches/Interpolate1D/fitpack/fpgivs.f
branches/Interpolate1D/fitpack/fpgrdi.f
branches/Interpolate1D/fitpack/fpgrpa.f
branches/Interpolate1D/fitpack/fpgrre.f
branches/Interpolate1D/fitpack/fpgrsp.f
branches/Interpolate1D/fitpack/fpinst.f
branches/Interpolate1D/fitpack/fpintb.f
branches/Interpolate1D/fitpack/fpknot.f
branches/Interpolate1D/fitpack/fpopdi.f
branches/Interpolate1D/fitpack/fpopsp.f
branches/Interpolate1D/fitpack/fporde.f
branches/Interpolate1D/fitpack/fppara.f
branches/Interpolate1D/fitpack/fppasu.f
branches/Interpolate1D/fitpack/fpperi.f
branches/Interpolate1D/fitpack/fppocu.f
branches/Interpolate1D/fitpack/fppogr.f
branches/Interpolate1D/fitpack/fppola.f
branches/Interpolate1D/fitpack/fprank.f
branches/Interpolate1D/fitpack/fprati.f
branches/Interpolate1D/fitpack/fpregr.f
branches/Interpolate1D/fitpack/fprota.f
branches/Interpolate1D/fitpack/fprppo.f
branches/Interpolate1D/fitpack/fprpsp.f
branches/Interpolate1D/fitpack/fpseno.f
branches/Interpolate1D/fitpack/fpspgr.f
branches/Interpolate1D/fitpack/fpsphe.f
branches/Interpolate1D/fitpack/fpsuev.f
branches/Interpolate1D/fitpack/fpsurf.f
branches/Interpolate1D/fitpack/fpsysy.f
branches/Interpolate1D/fitpack/fptrnp.f
branches/Interpolate1D/fitpack/fptrpe.f
branches/Interpolate1D/fitpack/insert.f
branches/Interpolate1D/fitpack/parcur.f
branches/Interpolate1D/fitpack/parder.f
branches/Interpolate1D/fitpack/parsur.f
branches/Interpolate1D/fitpack/percur.f
branches/Interpolate1D/fitpack/pogrid.f
branches/Interpolate1D/fitpack/polar.f
branches/Interpolate1D/fitpack/profil.f
branches/Interpolate1D/fitpack/regrid.f
branches/Interpolate1D/fitpack/spalde.f
branches/Interpolate1D/fitpack/spgrid.f
branches/Interpolate1D/fitpack/sphere.f
branches/Interpolate1D/fitpack/splder.f
branches/Interpolate1D/fitpack/splev.f
branches/Interpolate1D/fitpack/splint.f
branches/Interpolate1D/fitpack/sproot.f
branches/Interpolate1D/fitpack/surev.f
branches/Interpolate1D/fitpack/surfit.f
branches/Interpolate1D/multipack.h
Removed:
branches/Interpolate1D/Interpolate1D.pyc
branches/Interpolate1D/__init__fit.py
branches/Interpolate1D/_interpolate.pyd
branches/Interpolate1D/build/
branches/Interpolate1D/dfitpack.py
branches/Interpolate1D/dfitpack.pyd
branches/Interpolate1D/fitpack_wrapper.pyc
branches/Interpolate1D/info_fit.py
branches/Interpolate1D/interpolate_wrapper.pyc
Modified:
branches/Interpolate1D/Interpolate1D.py
branches/Interpolate1D/__init__.py
branches/Interpolate1D/fitpack_wrapper.py
branches/Interpolate1D/interpolate_wrapper.py
branches/Interpolate1D/setup.py
Log:
Interpolate1D builds correctly, unnecessary files removed, other improvements.
Modified: branches/Interpolate1D/Interpolate1D.py
===================================================================
--- branches/Interpolate1D/Interpolate1D.py 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/Interpolate1D.py 2008-07-18 19:44:12 UTC (rev 4550)
@@ -11,6 +11,7 @@
# fixme: use this to ensure proper type of all inputs and outputs in Interpolate1D
def make_array_safe(ary, typecode=np.float64):
+ # fixme: could pick correct typecode
ary = np.atleast_1d(np.asarray(ary, typecode))
if not ary.flags['CONTIGUOUS']:
ary = ary.copy()
@@ -18,44 +19,98 @@
class Interpolate1D(object):
- # see enthought.interpolate
- # fixme: Handle other data types.
- def __init__(self, x, y, k=1, kind='linear', low=None, high=None):
-
+ def __init__(self, x, y, kind='linear', low=np.NaN, high=np.NaN, kindkw={}, lowkw={}, highkw={}, missing_data=[None, np.NaN]):
+ """
+ Object for interpolation of 1D data.
+
+ REQUIRED ARGUMENTS:
+
+ x -- list or NumPy array
+ x includes the x-values for the data set to
+ interpolate from. It must be sorted in
+ ascending order
+
+ y -- list or NumPy array
+ y includes the y-values for the data set to
+ interpolate from.
+
+ OPTIONAL ARGUMENTS:
+
+ kind -- Usu. function or string. But can be any type.
+ Specifies the type of extrapolation to use for values within
+ the range of x. If a string is passed, it will look for an object
+ or function with that name and call it when evaluating. If
+ a function or object is passed, it will be called when interpolating.
+ If nothing else, assumes the argument is intended as a value
+ to be returned for all arguments. Defaults to linear interpolation.
+
+ kindkw -- dictionary
+ If kind is a class, function or string, additional keyword arguments
+ may be needed (example: if you want a 2nd order spline, kind = 'spline'
+ and kindkw = {'k' : 2}.
+
+ low (high) -- same as for kind
+ Same options as for 'kind'. Defaults to returning numpy.NaN ('not
+ a number') for all values outside the range of x.
+
+
+
+ """
# fixme: Handle checking if they are the correct size.
self._x = make_array_safe(x).copy()
self._y = make_array_safe(y).copy()
- self._xdtype = type(self._x[0])
- self._ydtype = type(self._y[0])
-
- assert( len(x) == len(y) , "x and y must be of the same length" )
- assert( x.ndim == 1 , "x must be one-dimensional" )
- assert( y.ndim == 1 , "y must be one-dimensional" )
+
+ assert len(x) == len(y) , "x and y must be of the same length"
+ assert x.ndim == 1 , "x must be one-dimensional"
+ assert y.ndim == 1 , "y must be one-dimensional"
# fixme: let y be 2-dimensional. Involves reworking of Interpolate1D.__call__
# because Spline enumerates y along the last, rather then first, axis,
# while concatenate works along first axis
- self.kind = self._init_interp_method(self._x, self._y, k, kind)
- self.low = self._init_interp_method(self._x, self._y, k, low)
- self.high = self._init_interp_method(self._x, self._y, k, high)
+ self.kind = self._init_interp_method(self._x, self._y, kind, kindkw)
+ self.low = self._init_interp_method(self._x, self._y, low, lowkw)
+ self.high = self._init_interp_method(self._x, self._y, high, highkw)
- def _init_interp_method(self, x, y, k, interp_arg):
+ def _format_array(x, y, missing_data=[None, np.NaN]):
+ # fixme: don't allow copying multiple times.
+
+ assert len(x) > 0 and len(y) > 0 , "interpolation does not support\
+ array of length 0"
+ assert len(x) == len(y) , "x and y must be of the same length"
+ mask = [((xi not in missing_data) and (y[i] not in missing_data)) \
+ for i, xi in enumerate(x) ]
+ if isinstance(x, list):
+ x = [x[i] for (i, good_data) in enumerate(mask) if good_data]
+ else:
+ x = x[mask]
+ if isinstance(y, list):
+ y = [y[i] for (i, good_data) in enumerate(mask) if good_data]
+ else:
+ y = y[mask]
+ self._xdtype = type(x[0])
+ self._x = make_array_safe(x, _xdtype).copy()
+ self._ydtype = type(y[0])
+ self._y = make_array_safe(y, _ydtype).copy()
+
+ assert self._x.ndim == 1 , "x must be one-dimensional"
+ assert self._y.ndim == 1 , "y must be one-dimensional"
+
+ def _init_interp_method(self, x, y, interp_arg, kw):
from inspect import isclass, isfunction
if interp_arg in ['linear', 'logarithmic', 'block', 'block_average_above']:
func = {'linear':linear, 'logarithmic':logarithmic, 'block':block, \
'block_average_above':block_average_above}[interp_arg]
- result = lambda new_x : func(self._x, self._y, new_x)
+ result = lambda new_x : func(self._x, self._y, new_x, **kw)
elif interp_arg in ['Spline', Spline, 'spline']:
- result = Spline(self._x, self._y, k=k)
+ result = Spline(self._x, self._y, **kw)
elif isfunction(interp_arg):
- result = interp_arg
+ result = lambda new_x : interp_arg(new_x, **kw)
elif isclass(interp_arg):
- result = interp_arg(x, y)
+ result = interp_arg(x, y, **kw)
else:
- print "warning: defaulting on extrapolation"
result = np.vectorize(lambda new_x : interp_arg)
return result
@@ -65,10 +120,9 @@
low_mask = x<self._x[0]
high_mask = x>self._x[-1]
interp_mask = (~low_mask) & (~high_mask)
-
- # hack, since getting an error when self.low or self.high gets 0-length array
- # and they return None or NaN
- if len(x[low_mask]) == 0: new_low=np.array([])
+
+ if len(x[low_mask]) == 0: new_low=np.array([]) # hack, since vectorize is failing
+ # work on lists/arrays of length 0
else: new_low = self.low(x[low_mask])
if len(x[interp_mask])==0: new_interp=np.array([])
else: new_interp = self.kind(x[interp_mask])
@@ -86,45 +140,89 @@
def assertAllclose(self, x, y):
self.assert_(np.allclose(make_array_safe(x), make_array_safe(y)))
- # fixme: run the test contained in the wrapper modules
+ def test__interpolate_wrapper(self):
+ print "\n\nTESTING _interpolate_wrapper MODULE"
+ from interpolate_wrapper import Test
+ T = Test()
+ T.runTest()
- def test_Interp_linearSpl(self):
- #return
+ def test__fitpack_wrapper(self):
+ print "\n\nTESTING _fitpack_wrapper MODULE"
+ from fitpack_wrapper import Test
+ T = Test()
+ T.runTest()
+
+ def test_spline1_defaultExt(self):
+ # make sure : spline order 1 (linear) interpolation works correctly
+ # make sure : default extrapolation works
print "\n\nTESTING LINEAR (1st ORDER) SPLINE"
- N = 7
+ N = 7 # must be > 5
x = np.arange(N)
y = np.arange(N)
+ interp_func = Interpolate1D(x, y, kind='Spline', kindkw={'k':1}, low=None, high=599.73)
+ new_x = np.arange(N+1)-0.5
+ new_y = interp_func(new_x)
+
+ self.assertAllclose(new_y[1:5], [0.5, 1.5, 2.5, 3.5])
+ self.assert_(new_y[0] == None)
+ self.assert_(new_y[-1] == 599.73)
+
+ def test_spline2(self):
+ print "\n\nTESTING 2nd ORDER SPLINE"
+ # make sure : order-2 splines work on linear data
+ N = 7 #must be > 5
+ x = np.arange(N)
+ y = np.arange(N)
T1 = time.clock()
- interp_func = Interpolate1D(x, y, k=1, kind='Spline', low=None, high=None)
+ interp_func = Interpolate1D(x, y, kind='Spline', kindkw={'k':2}, low='spline', high='spline')
T2 = time.clock()
- print 'time to create linear interp function: ', T2 - T1
- new_x = np.arange(N)-0.5
+ print "time to create 2nd order spline interp function with N = %i: " % N, T2 - T1
+ new_x = np.arange(N+1)-0.5
t1 = time.clock()
new_y = interp_func(new_x)
t2 = time.clock()
- print '1d interp (sec):', t2 - t1
+ print "time to evaluate 2nd order spline interp function with N = %i: " % N, t2 - t1
+ self.assertAllclose(new_x, new_y)
- print "new_y: ", new_y
- self.assertAllclose(new_y[1:5], [0.5, 1.5, 2.5, 3.5])
- self.assert_(new_y[0] == None)
+ # make sure for non-linear data
+ N = 7
+ x = np.arange(N)
+ y = x**2
+ interp_func = Interpolate1D(x, y, kind='Spline', kindkw={'k':2}, low='spline', high='spline')
+ new_x = np.arange(N+1)-0.5
+ new_y = interp_func(new_x)
+ self.assertAllclose(new_x**2, new_y)
+
def test_linear(self):
+ # make sure : linear interpolation works
+ # make sure : linear extrapolation works
print "\n\nTESTING LINEAR INTERPOLATION"
N = 7
x = arange(N)
y = arange(N)
new_x = arange(N+1)-0.5
T1 = time.clock()
- interp_func = Interpolate1D(x, y, kind='linear', low=None, high=None)
+ interp_func = Interpolate1D(x, y, kind='linear', low='linear', high='linear')
T2 = time.clock()
- print 'time to create linear interp function: ', T2 - T1
+ print "time to create linear interp function with N = %i: " % N, T2 - T1
t1 = time.clock()
new_y = interp_func(new_x)
t2 = time.clock()
- print '1d interp (sec):', t2 - t1
+ print "time to create linear interp function with N = %i: " % N, t2 - t1
- self.assertAllclose(new_y[1:5], [0.5, 1.5, 2.5, 3.5])
- self.assert_(new_y[0] == None)
+ self.assertAllclose(new_x, new_y)
+ def test_noLow(self):
+ # make sure : having the out-of-range elements in new_x is fine
+ # there was a bug with this
+ N = 5
+ x = arange(N)
+ y = arange(N)
+ new_x = arange(1,N-1)+.2
+ interp_func = Interpolate1D(x, y, kind='linear', low='linear', high=np.NaN)
+ new_y = interp_func(new_x)
+ self.assertAllclose(new_x, new_y)
+
if __name__ == '__main__':
unittest.main()
\ No newline at end of file
Deleted: branches/Interpolate1D/Interpolate1D.pyc
===================================================================
(Binary files differ)
Modified: branches/Interpolate1D/__init__.py
===================================================================
--- branches/Interpolate1D/__init__.py 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/__init__.py 2008-07-18 19:44:12 UTC (rev 4550)
@@ -1,3 +1,5 @@
-from interpolate_wrapper import linear, logarithmicm, block_average_above
+
+from interpolate_wrapper import linear, logarithmic, block, block_average_above
+from fitpack_wrapper import Spline
from Interpolate1D import Interpolate1D
\ No newline at end of file
Deleted: branches/Interpolate1D/__init__fit.py
===================================================================
--- branches/Interpolate1D/__init__fit.py 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/__init__fit.py 2008-07-18 19:44:12 UTC (rev 4550)
@@ -1,15 +0,0 @@
-#
-# interpolate - Interpolation Tools
-#
-
-from info import __doc__
-
-#from interpolate import *
-#from fitpack import *
-
-# New interface to fitpack library:
-from fitpack2 import *
-
-__all__ = filter(lambda s:not s.startswith('_'),dir())
-from numpy.testing import NumpyTest
-test = NumpyTest().test
Deleted: branches/Interpolate1D/_interpolate.pyd
===================================================================
(Binary files differ)
Deleted: branches/Interpolate1D/dfitpack.py
===================================================================
--- branches/Interpolate1D/dfitpack.py 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/dfitpack.py 2008-07-18 19:44:12 UTC (rev 4550)
@@ -1,7 +0,0 @@
-def __bootstrap__():
- global __bootstrap__, __loader__, __file__
- import sys, pkg_resources, imp
- __file__ = pkg_resources.resource_filename(__name__,'dfitpack.pyd')
- del __bootstrap__, __loader__
- imp.load_dynamic(__name__,__file__)
-__bootstrap__()
Deleted: branches/Interpolate1D/dfitpack.pyd
===================================================================
(Binary files differ)
Added: branches/Interpolate1D/fitpack/Makefile
===================================================================
--- branches/Interpolate1D/fitpack/Makefile 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/fitpack/Makefile 2008-07-18 19:44:12 UTC (rev 4550)
@@ -0,0 +1,19 @@
+# Makefile that builts a library lib$(LIB).a from all
+# of the Fortran files found in the current directory.
+# Usage: make LIB=<libname>
+# Pearu
+
+OBJ=$(patsubst %.f,%.o,$(shell ls *.f))
+all: lib$(LIB).a
+$(OBJ):
+ $(FC) -c $(FFLAGS) $(FSHARED) $(patsubst %.o,%.f,$(@F)) -o $@
+lib$(LIB).a: $(OBJ)
+ $(AR) rus lib$(LIB).a $?
+clean:
+ rm *.o
+
+
+
+
+
+
Added: branches/Interpolate1D/fitpack/README
===================================================================
--- branches/Interpolate1D/fitpack/README 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/fitpack/README 2008-07-18 19:44:12 UTC (rev 4550)
@@ -0,0 +1,3 @@
+- ddierckx is a 'real*8' version of dierckx
+ generated by Pearu Peterson <pearu@ioc.ee>.
+- dierckx (in netlib) is fitpack by P. Dierckx
Added: branches/Interpolate1D/fitpack/bispev.f
===================================================================
--- branches/Interpolate1D/fitpack/bispev.f 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/fitpack/bispev.f 2008-07-18 19:44:12 UTC (rev 4550)
@@ -0,0 +1,103 @@
+ subroutine bispev(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk,lwrk,
+ * iwrk,kwrk,ier)
+c subroutine bispev evaluates on a grid (x(i),y(j)),i=1,...,mx; j=1,...
+c ,my a bivariate spline s(x,y) of degrees kx and ky, given in the
+c b-spline representation.
+c
+c calling sequence:
+c call bispev(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk,lwrk,
+c * iwrk,kwrk,ier)
+c
+c input parameters:
+c tx : real array, length nx, which contains the position of the
+c knots in the x-direction.
+c nx : integer, giving the total number of knots in the x-direction
+c ty : real array, length ny, which contains the position of the
+c knots in the y-direction.
+c ny : integer, giving the total number of knots in the y-direction
+c c : real array, length (nx-kx-1)*(ny-ky-1), which contains the
+c b-spline coefficients.
+c kx,ky : integer values, giving the degrees of the spline.
+c x : real array of dimension (mx).
+c before entry x(i) must be set to the x co-ordinate of the
+c i-th grid point along the x-axis.
+c tx(kx+1)<=x(i-1)<=x(i)<=tx(nx-kx), i=2,...,mx.
+c mx : on entry mx must specify the number of grid points along
+c the x-axis. mx >=1.
+c y : real array of dimension (my).
+c before entry y(j) must be set to the y co-ordinate of the
+c j-th grid point along the y-axis.
+c ty(ky+1)<=y(j-1)<=y(j)<=ty(ny-ky), j=2,...,my.
+c my : on entry my must specify the number of grid points along
+c the y-axis. my >=1.
+c wrk : real array of dimension lwrk. used as workspace.
+c lwrk : integer, specifying the dimension of wrk.
+c lwrk >= mx*(kx+1)+my*(ky+1)
+c iwrk : integer array of dimension kwrk. used as workspace.
+c kwrk : integer, specifying the dimension of iwrk. kwrk >= mx+my.
+c
+c output parameters:
+c z : real array of dimension (mx*my).
+c on succesful exit z(my*(i-1)+j) contains the value of s(x,y)
+c at the point (x(i),y(j)),i=1,...,mx;j=1,...,my.
+c ier : integer error flag
+c ier=0 : normal return
+c ier=10: invalid input data (see restrictions)
+c
+c restrictions:
+c mx >=1, my >=1, lwrk>=mx*(kx+1)+my*(ky+1), kwrk>=mx+my
+c tx(kx+1) <= x(i-1) <= x(i) <= tx(nx-kx), i=2,...,mx
+c ty(ky+1) <= y(j-1) <= y(j) <= ty(ny-ky), j=2,...,my
+c
+c other subroutines required:
+c fpbisp,fpbspl
+c
+c references :
+c de boor c : on calculating with b-splines, j. approximation theory
+c 6 (1972) 50-62.
+c cox m.g. : the numerical evaluation of b-splines, j. inst. maths
+c applics 10 (1972) 134-149.
+c dierckx p. : curve and surface fitting with splines, monographs on
+c numerical analysis, oxford university press, 1993.
+c
+c author :
+c p.dierckx
+c dept. computer science, k.u.leuven
+c celestijnenlaan 200a, b-3001 heverlee, belgium.
+c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
+c
+c latest update : march 1987
+c
+c ..scalar arguments..
+ integer nx,ny,kx,ky,mx,my,lwrk,kwrk,ier
+c ..array arguments..
+ integer iwrk(kwrk)
+ real*8 tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1)),x(mx),y(my),z(mx*my),
+ * wrk(lwrk)
+c ..local scalars..
+ integer i,iw,lwest
+c ..
+c before starting computations a data check is made. if the input data
+c are invalid control is immediately repassed to the calling program.
+ ier = 10
+ lwest = (kx+1)*mx+(ky+1)*my
+ if(lwrk.lt.lwest) go to 100
+ if(kwrk.lt.(mx+my)) go to 100
+ if (mx.lt.1) go to 100
+ if (mx.eq.1) go to 30
+ go to 10
+ 10 do 20 i=2,mx
+ if(x(i).lt.x(i-1)) go to 100
+ 20 continue
+ 30 if (my.lt.1) go to 100
+ if (my.eq.1) go to 60
+ go to 40
+ 40 do 50 i=2,my
+ if(y(i).lt.y(i-1)) go to 100
+ 50 continue
+ 60 ier = 0
+ iw = mx*(kx+1)+1
+ call fpbisp(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk(1),wrk(iw),
+ * iwrk(1),iwrk(mx+1))
+ 100 return
+ end
Added: branches/Interpolate1D/fitpack/clocur.f
===================================================================
--- branches/Interpolate1D/fitpack/clocur.f 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/fitpack/clocur.f 2008-07-18 19:44:12 UTC (rev 4550)
@@ -0,0 +1,352 @@
+ subroutine clocur(iopt,ipar,idim,m,u,mx,x,w,k,s,nest,n,t,nc,c,fp,
+ * wrk,lwrk,iwrk,ier)
+c given the ordered set of m points x(i) in the idim-dimensional space
+c with x(1)=x(m), and given also a corresponding set of strictly in-
+c creasing values u(i) and the set of positive numbers w(i),i=1,2,...,m
+c subroutine clocur determines a smooth approximating closed spline
+c curve s(u), i.e.
+c x1 = s1(u)
+c x2 = s2(u) u(1) <= u <= u(m)
+c .........
+c xidim = sidim(u)
+c with sj(u),j=1,2,...,idim periodic spline functions of degree k with
+c common knots t(j),j=1,2,...,n.
+c if ipar=1 the values u(i),i=1,2,...,m must be supplied by the user.
+c if ipar=0 these values are chosen automatically by clocur as
+c v(1) = 0
+c v(i) = v(i-1) + dist(x(i),x(i-1)) ,i=2,3,...,m
+c u(i) = v(i)/v(m) ,i=1,2,...,m
+c if iopt=-1 clocur calculates the weighted least-squares closed spline
+c curve according to a given set of knots.
+c if iopt>=0 the number of knots of the splines sj(u) and the position
+c t(j),j=1,2,...,n is chosen automatically by the routine. the smooth-
+c ness of s(u) is then achieved by minimalizing the discontinuity
+c jumps of the k-th derivative of s(u) at the knots t(j),j=k+2,k+3,...,
+c n-k-1. the amount of smoothness is determined by the condition that
+c f(p)=sum((w(i)*dist(x(i),s(u(i))))**2) be <= s, with s a given non-
+c negative constant, called the smoothing factor.
+c the fit s(u) is given in the b-spline representation and can be
+c evaluated by means of subroutine curev.
+c
+c calling sequence:
+c call clocur(iopt,ipar,idim,m,u,mx,x,w,k,s,nest,n,t,nc,c,
+c * fp,wrk,lwrk,iwrk,ier)
+c
+c parameters:
+c iopt : integer flag. on entry iopt must specify whether a weighted
+c least-squares closed spline curve (iopt=-1) or a smoothing
+c closed spline curve (iopt=0 or 1) must be determined. if
+c iopt=0 the routine will start with an initial set of knots
+c t(i)=u(1)+(u(m)-u(1))*(i-k-1),i=1,2,...,2*k+2. if iopt=1 the
+c routine will continue with the knots found at the last call.
+c attention: a call with iopt=1 must always be immediately
+c preceded by another call with iopt=1 or iopt=0.
+c unchanged on exit.
+c ipar : integer flag. on entry ipar must specify whether (ipar=1)
+c the user will supply the parameter values u(i),or whether
+c (ipar=0) these values are to be calculated by clocur.
+c unchanged on exit.
+c idim : integer. on entry idim must specify the dimension of the
+c curve. 0 < idim < 11.
+c unchanged on exit.
+c m : integer. on entry m must specify the number of data points.
+c m > 1. unchanged on exit.
+c u : real array of dimension at least (m). in case ipar=1,before
+c entry, u(i) must be set to the i-th value of the parameter
+c variable u for i=1,2,...,m. these values must then be
+c supplied in strictly ascending order and will be unchanged
+c on exit. in case ipar=0, on exit,the array will contain the
+c values u(i) as determined by clocur.
+c mx : integer. on entry mx must specify the actual dimension of
+c the array x as declared in the calling (sub)program. mx must
+c not be too small (see x). unchanged on exit.
+c x : real array of dimension at least idim*m.
+c before entry, x(idim*(i-1)+j) must contain the j-th coord-
+c inate of the i-th data point for i=1,2,...,m and j=1,2,...,
+c idim. since first and last data point must coincide it
+c means that x(j)=x(idim*(m-1)+j),j=1,2,...,idim.
+c unchanged on exit.
+c w : real array of dimension at least (m). before entry, w(i)
+c must be set to the i-th value in the set of weights. the
+c w(i) must be strictly positive. w(m) is not used.
+c unchanged on exit. see also further comments.
+c k : integer. on entry k must specify the degree of the splines.
+c 1<=k<=5. it is recommended to use cubic splines (k=3).
+c the user is strongly dissuaded from choosing k even,together
+c with a small s-value. unchanged on exit.
+c s : real.on entry (in case iopt>=0) s must specify the smoothing
+c factor. s >=0. unchanged on exit.
+c for advice on the choice of s see further comments.
+c nest : integer. on entry nest must contain an over-estimate of the
+c total number of knots of the splines returned, to indicate
+c the storage space available to the routine. nest >=2*k+2.
+c in most practical situation nest=m/2 will be sufficient.
+c always large enough is nest=m+2*k, the number of knots
+c needed for interpolation (s=0). unchanged on exit.
+c n : integer.
+c unless ier = 10 (in case iopt >=0), n will contain the
+c total number of knots of the smoothing spline curve returned
+c if the computation mode iopt=1 is used this value of n
+c should be left unchanged between subsequent calls.
+c in case iopt=-1, the value of n must be specified on entry.
+c t : real array of dimension at least (nest).
+c on succesful exit, this array will contain the knots of the
+c spline curve,i.e. the position of the interior knots t(k+2),
+c t(k+3),..,t(n-k-1) as well as the position of the additional
+c t(1),t(2),..,t(k+1)=u(1) and u(m)=t(n-k),...,t(n) needed for
+c the b-spline representation.
+c if the computation mode iopt=1 is used, the values of t(1),
+c t(2),...,t(n) should be left unchanged between subsequent
+c calls. if the computation mode iopt=-1 is used, the values
+c t(k+2),...,t(n-k-1) must be supplied by the user, before
+c entry. see also the restrictions (ier=10).
+c nc : integer. on entry nc must specify the actual dimension of
+c the array c as declared in the calling (sub)program. nc
+c must not be too small (see c). unchanged on exit.
+c c : real array of dimension at least (nest*idim).
+c on succesful exit, this array will contain the coefficients
+c in the b-spline representation of the spline curve s(u),i.e.
+c the b-spline coefficients of the spline sj(u) will be given
+c in c(n*(j-1)+i),i=1,2,...,n-k-1 for j=1,2,...,idim.
+c fp : real. unless ier = 10, fp contains the weighted sum of
+c squared residuals of the spline curve returned.
+c wrk : real array of dimension at least m*(k+1)+nest*(7+idim+5*k).
+c used as working space. if the computation mode iopt=1 is
+c used, the values wrk(1),...,wrk(n) should be left unchanged
+c between subsequent calls.
+c lwrk : integer. on entry,lwrk must specify the actual dimension of
+c the array wrk as declared in the calling (sub)program. lwrk
+c must not be too small (see wrk). unchanged on exit.
+c iwrk : integer array of dimension at least (nest).
+c used as working space. if the computation mode iopt=1 is
+c used,the values iwrk(1),...,iwrk(n) should be left unchanged
+c between subsequent calls.
+c ier : integer. unless the routine detects an error, ier contains a
+c non-positive value on exit, i.e.
+c ier=0 : normal return. the close curve returned has a residual
+c sum of squares fp such that abs(fp-s)/s <= tol with tol a
+c relative tolerance set to 0.001 by the program.
+c ier=-1 : normal return. the curve returned is an interpolating
+c spline curve (fp=0).
+c ier=-2 : normal return. the curve returned is the weighted least-
+c squares point,i.e. each spline sj(u) is a constant. in
+c this extreme case fp gives the upper bound fp0 for the
+c smoothing factor s.
+c ier=1 : error. the required storage space exceeds the available
+c storage space, as specified by the parameter nest.
+c probably causes : nest too small. if nest is already
+c large (say nest > m/2), it may also indicate that s is
+c too small
+c the approximation returned is the least-squares closed
+c curve according to the knots t(1),t(2),...,t(n). (n=nest)
+c the parameter fp gives the corresponding weighted sum of
+c squared residuals (fp>s).
+c ier=2 : error. a theoretically impossible result was found during
+c the iteration proces for finding a smoothing curve with
+c fp = s. probably causes : s too small.
+c there is an approximation returned but the corresponding
+c weighted sum of squared residuals does not satisfy the
+c condition abs(fp-s)/s < tol.
+c ier=3 : error. the maximal number of iterations maxit (set to 20
+c by the program) allowed for finding a smoothing curve
+c with fp=s has been reached. probably causes : s too small
+c there is an approximation returned but the corresponding
+c weighted sum of squared residuals does not satisfy the
+c condition abs(fp-s)/s < tol.
+c ier=10 : error. on entry, the input data are controlled on validity
+c the following restrictions must be satisfied.
+c -1<=iopt<=1, 1<=k<=5, m>1, nest>2*k+2, w(i)>0,i=1,2,...,m
+c 0<=ipar<=1, 0<idim<=10, lwrk>=(k+1)*m+nest*(7+idim+5*k),
+c nc>=nest*idim, x(j)=x(idim*(m-1)+j), j=1,2,...,idim
+c if ipar=0: sum j=1,idim (x(i*idim+j)-x((i-1)*idim+j))**2>0
+c i=1,2,...,m-1.
+c if ipar=1: u(1)<u(2)<...<u(m)
+c if iopt=-1: 2*k+2<=n<=min(nest,m+2*k)
+c u(1)<t(k+2)<t(k+3)<...<t(n-k-1)<u(m)
+c (u(1)=0 and u(m)=1 in case ipar=0)
+c the schoenberg-whitney conditions, i.e. there
+c must be a subset of data points uu(j) with
+c uu(j) = u(i) or u(i)+(u(m)-u(1)) such that
+c t(j) < uu(j) < t(j+k+1), j=k+1,...,n-k-1
+c if iopt>=0: s>=0
+c if s=0 : nest >= m+2*k
+c if one of these conditions is found to be violated,control
+c is immediately repassed to the calling program. in that
+c case there is no approximation returned.
+c
+c further comments:
+c by means of the parameter s, the user can control the tradeoff
+c between closeness of fit and smoothness of fit of the approximation.
+c if s is too large, the curve will be too smooth and signal will be
+c lost ; if s is too small the curve will pick up too much noise. in
+c the extreme cases the program will return an interpolating curve if
+c s=0 and the weighted least-squares point if s is very large.
+c between these extremes, a properly chosen s will result in a good
+c compromise between closeness of fit and smoothness of fit.
+c to decide whether an approximation, corresponding to a certain s is
+c satisfactory the user is highly recommended to inspect the fits
+c graphically.
+c recommended values for s depend on the weights w(i). if these are
+c taken as 1/d(i) with d(i) an estimate of the standard deviation of
+c x(i), a good s-value should be found in the range (m-sqrt(2*m),m+
+c sqrt(2*m)). if nothing is known about the statistical error in x(i)
+c each w(i) can be set equal to one and s determined by trial and
+c error, taking account of the comments above. the best is then to
+c start with a very large value of s ( to determine the weighted
+c least-squares point and the upper bound fp0 for s) and then to
+c progressively decrease the value of s ( say by a factor 10 in the
+c beginning, i.e. s=fp0/10, fp0/100,...and more carefully as the
+c approximating curve shows more detail) to obtain closer fits.
+c to economize the search for a good s-value the program provides with
+c different modes of computation. at the first call of the routine, or
+c whenever he wants to restart with the initial set of knots the user
+c must set iopt=0.
+c if iopt=1 the program will continue with the set of knots found at
+c the last call of the routine. this will save a lot of computation
+c time if clocur is called repeatedly for different values of s.
+c the number of knots of the spline returned and their location will
+c depend on the value of s and on the complexity of the shape of the
+c curve underlying the data. but, if the computation mode iopt=1 is
+c used, the knots returned may also depend on the s-values at previous
+c calls (if these were smaller). therefore, if after a number of
+c trials with different s-values and iopt=1, the user can finally
+c accept a fit as satisfactory, it may be worthwhile for him to call
+c clocur once more with the selected value for s but now with iopt=0.
+c indeed, clocur may then return an approximation of the same quality
+c of fit but with fewer knots and therefore better if data reduction
+c is also an important objective for the user.
+c
+c the form of the approximating curve can strongly be affected by
+c the choice of the parameter values u(i). if there is no physical
+c reason for choosing a particular parameter u, often good results
+c will be obtained with the choice of clocur(in case ipar=0), i.e.
+c v(1)=0, v(i)=v(i-1)+q(i), i=2,...,m, u(i)=v(i)/v(m), i=1,..,m
+c where
+c q(i)= sqrt(sum j=1,idim (xj(i)-xj(i-1))**2 )
+c other possibilities for q(i) are
+c q(i)= sum j=1,idim (xj(i)-xj(i-1))**2
+c q(i)= sum j=1,idim abs(xj(i)-xj(i-1))
+c q(i)= max j=1,idim abs(xj(i)-xj(i-1))
+c q(i)= 1
+c
+c
+c other subroutines required:
+c fpbacp,fpbspl,fpchep,fpclos,fpdisc,fpgivs,fpknot,fprati,fprota
+c
+c references:
+c dierckx p. : algorithms for smoothing data with periodic and
+c parametric splines, computer graphics and image
+c processing 20 (1982) 171-184.
+c dierckx p. : algorithms for smoothing data with periodic and param-
+c etric splines, report tw55, dept. computer science,
+c k.u.leuven, 1981.
+c dierckx p. : curve and surface fitting with splines, monographs on
+c numerical analysis, oxford university press, 1993.
+c
+c author:
+c p.dierckx
+c dept. computer science, k.u. leuven
+c celestijnenlaan 200a, b-3001 heverlee, belgium.
+c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
+c
+c creation date : may 1979
+c latest update : march 1987
+c
+c ..
+c ..scalar arguments..
+ real*8 s,fp
+ integer iopt,ipar,idim,m,mx,k,nest,n,nc,lwrk,ier
+c ..array arguments..
+ real*8 u(m),x(mx),w(m),t(nest),c(nc),wrk(lwrk)
+ integer iwrk(nest)
+c ..local scalars..
+ real*8 per,tol,dist
+ integer i,ia1,ia2,ib,ifp,ig1,ig2,iq,iz,i1,i2,j1,j2,k1,k2,lwest,
+ * maxit,m1,nmin,ncc,j
+c ..function references..
+ real*8 sqrt
+c we set up the parameters tol and maxit
+ maxit = 20
+ tol = 0.1e-02
+c before starting computations a data check is made. if the input data
+c are invalid, control is immediately repassed to the calling program.
+ ier = 10
+ if(iopt.lt.(-1) .or. iopt.gt.1) go to 90
+ if(ipar.lt.0 .or. ipar.gt.1) go to 90
+ if(idim.le.0 .or. idim.gt.10) go to 90
+ if(k.le.0 .or. k.gt.5) go to 90
+ k1 = k+1
+ k2 = k1+1
+ nmin = 2*k1
+ if(m.lt.2 .or. nest.lt.nmin) go to 90
+ ncc = nest*idim
+ if(mx.lt.m*idim .or. nc.lt.ncc) go to 90
+ lwest = m*k1+nest*(7+idim+5*k)
+ if(lwrk.lt.lwest) go to 90
+ i1 = idim
+ i2 = m*idim
+ do 5 j=1,idim
+ if(x(i1).ne.x(i2)) go to 90
+ i1 = i1-1
+ i2 = i2-1
+ 5 continue
+ if(ipar.ne.0 .or. iopt.gt.0) go to 40
+ i1 = 0
+ i2 = idim
+ u(1) = 0.
+ do 20 i=2,m
+ dist = 0.
+ do 10 j1=1,idim
+ i1 = i1+1
+ i2 = i2+1
+ dist = dist+(x(i2)-x(i1))**2
+ 10 continue
+ u(i) = u(i-1)+sqrt(dist)
+ 20 continue
+ if(u(m).le.0.) go to 90
+ do 30 i=2,m
+ u(i) = u(i)/u(m)
+ 30 continue
+ u(m) = 0.1e+01
+ 40 if(w(1).le.0.) go to 90
+ m1 = m-1
+ do 50 i=1,m1
+ if(u(i).ge.u(i+1) .or. w(i).le.0.) go to 90
+ 50 continue
+ if(iopt.ge.0) go to 70
+ if(n.le.nmin .or. n.gt.nest) go to 90
+ per = u(m)-u(1)
+ j1 = k1
+ t(j1) = u(1)
+ i1 = n-k
+ t(i1) = u(m)
+ j2 = j1
+ i2 = i1
+ do 60 i=1,k
+ i1 = i1+1
+ i2 = i2-1
+ j1 = j1+1
+ j2 = j2-1
+ t(j2) = t(i2)-per
+ t(i1) = t(j1)+per
+ 60 continue
+ call fpchep(u,m,t,n,k,ier)
+ if (ier.eq.0) go to 80
+ go to 90
+ 70 if(s.lt.0.) go to 90
+ if(s.eq.0. .and. nest.lt.(m+2*k)) go to 90
+ ier = 0
+c we partition the working space and determine the spline approximation.
+ 80 ifp = 1
+ iz = ifp+nest
+ ia1 = iz+ncc
+ ia2 = ia1+nest*k1
+ ib = ia2+nest*k
+ ig1 = ib+nest*k2
+ ig2 = ig1+nest*k2
+ iq = ig2+nest*k1
+ call fpclos(iopt,idim,m,u,mx,x,w,k,s,nest,tol,maxit,k1,k2,n,t,
+ * ncc,c,fp,wrk(ifp),wrk(iz),wrk(ia1),wrk(ia2),wrk(ib),wrk(ig1),
+ * wrk(ig2),wrk(iq),iwrk,ier)
+ 90 return
+ end
Added: branches/Interpolate1D/fitpack/cocosp.f
===================================================================
--- branches/Interpolate1D/fitpack/cocosp.f 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/fitpack/cocosp.f 2008-07-18 19:44:12 UTC (rev 4550)
@@ -0,0 +1,180 @@
+ subroutine cocosp(m,x,y,w,n,t,e,maxtr,maxbin,c,sq,sx,bind,wrk,
+ * lwrk,iwrk,kwrk,ier)
+c given the set of data points (x(i),y(i)) and the set of positive
+c numbers w(i),i=1,2,...,m, subroutine cocosp determines the weighted
+c least-squares cubic spline s(x) with given knots t(j),j=1,2,...,n
+c which satisfies the following concavity/convexity conditions
+c s''(t(j+3))*e(j) <= 0, j=1,2,...n-6
+c the fit is given in the b-spline representation( b-spline coef-
+c ficients c(j),j=1,2,...n-4) and can be evaluated by means of
+c subroutine splev.
+c
+c calling sequence:
+c call cocosp(m,x,y,w,n,t,e,maxtr,maxbin,c,sq,sx,bind,wrk,
+c * lwrk,iwrk,kwrk,ier)
+c
+c parameters:
+c m : integer. on entry m must specify the number of data points.
+c m > 3. unchanged on exit.
+c x : real array of dimension at least (m). before entry, x(i)
+c must be set to the i-th value of the independent variable x,
+c for i=1,2,...,m. these values must be supplied in strictly
+c ascending order. unchanged on exit.
+c y : real array of dimension at least (m). before entry, y(i)
+c must be set to the i-th value of the dependent variable y,
+c for i=1,2,...,m. unchanged on exit.
+c w : real array of dimension at least (m). before entry, w(i)
+c must be set to the i-th value in the set of weights. the
+c w(i) must be strictly positive. unchanged on exit.
+c n : integer. on entry n must contain the total number of knots
+c of the cubic spline. m+4>=n>=8. unchanged on exit.
+c t : real array of dimension at least (n). before entry, this
+c array must contain the knots of the spline, i.e. the position
+c of the interior knots t(5),t(6),...,t(n-4) as well as the
+c position of the boundary knots t(1),t(2),t(3),t(4) and t(n-3)
+c t(n-2),t(n-1),t(n) needed for the b-spline representation.
+c unchanged on exit. see also the restrictions (ier=10).
+c e : real array of dimension at least (n). before entry, e(j)
+c must be set to 1 if s(x) must be locally concave at t(j+3),
+c to (-1) if s(x) must be locally convex at t(j+3) and to 0
+c if no convexity constraint is imposed at t(j+3),j=1,2,..,n-6.
+c e(n-5),...,e(n) are not used. unchanged on exit.
+c maxtr : integer. on entry maxtr must contain an over-estimate of the
+c total number of records in the used tree structure, to indic-
+c ate the storage space available to the routine. maxtr >=1
+c in most practical situation maxtr=100 will be sufficient.
+c always large enough is
+c n-5 n-6
+c maxtr = ( ) + ( ) with l the greatest
+c l l+1
+c integer <= (n-6)/2 . unchanged on exit.
+c maxbin: integer. on entry maxbin must contain an over-estimate of the
+c number of knots where s(x) will have a zero second derivative
+c maxbin >=1. in most practical situation maxbin = 10 will be
+c sufficient. always large enough is maxbin=n-6.
+c unchanged on exit.
+c c : real array of dimension at least (n).
+c on succesful exit, this array will contain the coefficients
+c c(1),c(2),..,c(n-4) in the b-spline representation of s(x)
+c sq : real. on succesful exit, sq contains the weighted sum of
+c squared residuals of the spline approximation returned.
+c sx : real array of dimension at least m. on succesful exit
+c this array will contain the spline values s(x(i)),i=1,...,m
+c bind : logical array of dimension at least (n). on succesful exit
+c this array will indicate the knots where s''(x)=0, i.e.
+c s''(t(j+3)) .eq. 0 if bind(j) = .true.
+c s''(t(j+3)) .ne. 0 if bind(j) = .false., j=1,2,...,n-6
+c wrk : real array of dimension at least m*4+n*7+maxbin*(maxbin+n+1)
+c used as working space.
+c lwrk : integer. on entry,lwrk must specify the actual dimension of
+c the array wrk as declared in the calling (sub)program.lwrk
+c must not be too small (see wrk). unchanged on exit.
+c iwrk : integer array of dimension at least (maxtr*4+2*(maxbin+1))
+c used as working space.
+c kwrk : integer. on entry,kwrk must specify the actual dimension of
+c the array iwrk as declared in the calling (sub)program. kwrk
+c must not be too small (see iwrk). unchanged on exit.
+c ier : integer. error flag
+c ier=0 : succesful exit.
+c ier>0 : abnormal termination: no approximation is returned
+c ier=1 : the number of knots where s''(x)=0 exceeds maxbin.
+c probably causes : maxbin too small.
+c ier=2 : the number of records in the tree structure exceeds
+c maxtr.
+c probably causes : maxtr too small.
+c ier=3 : the algoritm finds no solution to the posed quadratic
+c programming problem.
+c probably causes : rounding errors.
+c ier=10 : on entry, the input data are controlled on validity.
+c the following restrictions must be satisfied
+c m>3, maxtr>=1, maxbin>=1, 8<=n<=m+4,w(i) > 0,
+c x(1)<x(2)<...<x(m), t(1)<=t(2)<=t(3)<=t(4)<=x(1),
+c x(1)<t(5)<t(6)<...<t(n-4)<x(m)<=t(n-3)<=...<=t(n),
+c kwrk>=maxtr*4+2*(maxbin+1),
+c lwrk>=m*4+n*7+maxbin*(maxbin+n+1),
+c the schoenberg-whitney conditions, i.e. there must
+c be a subset of data points xx(j) such that
+c t(j) < xx(j) < t(j+4), j=1,2,...,n-4
+c if one of these restrictions is found to be violated
+c control is immediately repassed to the calling program
+c
+c
+c other subroutines required:
+c fpcosp,fpbspl,fpadno,fpdeno,fpseno,fpfrno,fpchec
+c
+c references:
+c dierckx p. : an algorithm for cubic spline fitting with convexity
+c constraints, computing 24 (1980) 349-371.
+c dierckx p. : an algorithm for least-squares cubic spline fitting
+c with convexity and concavity constraints, report tw39,
+c dept. computer science, k.u.leuven, 1978.
+c dierckx p. : curve and surface fitting with splines, monographs on
+c numerical analysis, oxford university press, 1993.
+c
+c author:
+c p. dierckx
+c dept. computer science, k.u.leuven
+c celestijnenlaan 200a, b-3001 heverlee, belgium.
+c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
+c
+c creation date : march 1978
+c latest update : march 1987.
+c
+c ..
+c ..scalar arguments..
+ real*8 sq
+ integer m,n,maxtr,maxbin,lwrk,kwrk,ier
+c ..array arguments..
+ real*8 x(m),y(m),w(m),t(n),e(n),c(n),sx(m),wrk(lwrk)
+ integer iwrk(kwrk)
+ logical bind(n)
+c ..local scalars..
+ integer i,ia,ib,ic,iq,iu,iz,izz,ji,jib,jjb,jl,jr,ju,kwest,
+ * lwest,mb,nm,n6
+ real*8 one
+c ..
+c set constant
+ one = 0.1e+01
+c before starting computations a data check is made. if the input data
+c are invalid, control is immediately repassed to the calling program.
+ ier = 10
+ if(m.lt.4 .or. n.lt.8) go to 40
+ if(maxtr.lt.1 .or. maxbin.lt.1) go to 40
+ lwest = 7*n+m*4+maxbin*(1+n+maxbin)
+ kwest = 4*maxtr+2*(maxbin+1)
+ if(lwrk.lt.lwest .or. kwrk.lt.kwest) go to 40
+ if(w(1).le.0.) go to 40
+ do 10 i=2,m
+ if(x(i-1).ge.x(i) .or. w(i).le.0.) go to 40
+ 10 continue
+ call fpchec(x,m,t,n,3,ier)
+ if (ier.eq.0) go to 20
+ go to 40
+c set numbers e(i)
+ 20 n6 = n-6
+ do 30 i=1,n6
+ if(e(i).gt.0.) e(i) = one
+ if(e(i).lt.0.) e(i) = -one
+ 30 continue
+c we partition the working space and determine the spline approximation
+ nm = n+maxbin
+ mb = maxbin+1
+ ia = 1
+ ib = ia+4*n
+ ic = ib+nm*maxbin
+ iz = ic+n
+ izz = iz+n
+ iu = izz+n
+ iq = iu+maxbin
+ ji = 1
+ ju = ji+maxtr
+ jl = ju+maxtr
+ jr = jl+maxtr
+ jjb = jr+maxtr
+ jib = jjb+mb
+ call fpcosp(m,x,y,w,n,t,e,maxtr,maxbin,c,sq,sx,bind,nm,mb,wrk(ia),
+ *
+ * wrk(ib),wrk(ic),wrk(iz),wrk(izz),wrk(iu),wrk(iq),iwrk(ji),
+ * iwrk(ju),iwrk(jl),iwrk(jr),iwrk(jjb),iwrk(jib),ier)
+ 40 return
+ end
Added: branches/Interpolate1D/fitpack/concon.f
===================================================================
--- branches/Interpolate1D/fitpack/concon.f 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/fitpack/concon.f 2008-07-18 19:44:12 UTC (rev 4550)
@@ -0,0 +1,233 @@
+ subroutine concon(iopt,m,x,y,w,v,s,nest,maxtr,maxbin,n,t,c,sq,
+ * sx,bind,wrk,lwrk,iwrk,kwrk,ier)
+c given the set of data points (x(i),y(i)) and the set of positive
+c numbers w(i), i=1,2,...,m,subroutine concon determines a cubic spline
+c approximation s(x) which satisfies the following local convexity
+c constraints s''(x(i))*v(i) <= 0, i=1,2,...,m.
+c the number of knots n and the position t(j),j=1,2,...n is chosen
+c automatically by the routine in a way that
+c sq = sum((w(i)*(y(i)-s(x(i))))**2) be <= s.
+c the fit is given in the b-spline representation (b-spline coef-
+c ficients c(j),j=1,2,...n-4) and can be evaluated by means of
+c subroutine splev.
+c
+c calling sequence:
+c
+c call concon(iopt,m,x,y,w,v,s,nest,maxtr,maxbin,n,t,c,sq,
+c * sx,bind,wrk,lwrk,iwrk,kwrk,ier)
+c
+c parameters:
+c iopt: integer flag.
+c if iopt=0, the routine will start with the minimal number of
+c knots to guarantee that the convexity conditions will be
+c satisfied. if iopt=1, the routine will continue with the set
+c of knots found at the last call of the routine.
+c attention: a call with iopt=1 must always be immediately
+c preceded by another call with iopt=1 or iopt=0.
+c unchanged on exit.
+c m : integer. on entry m must specify the number of data points.
+c m > 3. unchanged on exit.
+c x : real array of dimension at least (m). before entry, x(i)
+c must be set to the i-th value of the independent variable x,
+c for i=1,2,...,m. these values must be supplied in strictly
+c ascending order. unchanged on exit.
+c y : real array of dimension at least (m). before entry, y(i)
+c must be set to the i-th value of the dependent variable y,
+c for i=1,2,...,m. unchanged on exit.
+c w : real array of dimension at least (m). before entry, w(i)
+c must be set to the i-th value in the set of weights. the
+c w(i) must be strictly positive. unchanged on exit.
+c v : real array of dimension at least (m). before entry, v(i)
+c must be set to 1 if s(x) must be locally concave at x(i),
+c to (-1) if s(x) must be locally convex at x(i) and to 0
+c if no convexity constraint is imposed at x(i).
+c s : real. on entry s must specify an over-estimate for the
+c the weighted sum of squared residuals sq of the requested
+c spline. s >=0. unchanged on exit.
+c nest : integer. on entry nest must contain an over-estimate of the
+c total number of knots of the spline returned, to indicate
+c the storage space available to the routine. nest >=8.
+c in most practical situation nest=m/2 will be sufficient.
+c always large enough is nest=m+4. unchanged on exit.
+c maxtr : integer. on entry maxtr must contain an over-estimate of the
+c total number of records in the used tree structure, to indic-
+c ate the storage space available to the routine. maxtr >=1
+c in most practical situation maxtr=100 will be sufficient.
+c always large enough is
+c nest-5 nest-6
+c maxtr = ( ) + ( ) with l the greatest
+c l l+1
+c integer <= (nest-6)/2 . unchanged on exit.
+c maxbin: integer. on entry maxbin must contain an over-estimate of the
+c number of knots where s(x) will have a zero second derivative
+c maxbin >=1. in most practical situation maxbin = 10 will be
+c sufficient. always large enough is maxbin=nest-6.
+c unchanged on exit.
+c n : integer.
+c on exit with ier <=0, n will contain the total number of
+c knots of the spline approximation returned. if the comput-
+c ation mode iopt=1 is used this value of n should be left
+c unchanged between subsequent calls.
+c t : real array of dimension at least (nest).
+c on exit with ier<=0, this array will contain the knots of the
+c spline,i.e. the position of the interior knots t(5),t(6),...,
+c t(n-4) as well as the position of the additional knots
+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)
+c needed for the the b-spline representation.
+c if the computation mode iopt=1 is used, the values of t(1),
+c t(2),...,t(n) should be left unchanged between subsequent
+c calls.
+c c : real array of dimension at least (nest).
+c on succesful exit, this array will contain the coefficients
+c c(1),c(2),..,c(n-4) in the b-spline representation of s(x)
+c sq : real. unless ier>0 , sq contains the weighted sum of
+c squared residuals of the spline approximation returned.
+c sx : real array of dimension at least m. on exit with ier<=0
+c this array will contain the spline values s(x(i)),i=1,...,m
+c if the computation mode iopt=1 is used, the values of sx(1),
+c sx(2),...,sx(m) should be left unchanged between subsequent
+c calls.
+c bind: logical array of dimension at least nest. on exit with ier<=0
+c this array will indicate the knots where s''(x)=0, i.e.
+c s''(t(j+3)) .eq. 0 if bind(j) = .true.
+c s''(t(j+3)) .ne. 0 if bind(j) = .false., j=1,2,...,n-6
+c if the computation mode iopt=1 is used, the values of bind(1)
+c ,...,bind(n-6) should be left unchanged between subsequent
+c calls.
+c wrk : real array of dimension at least (m*4+nest*8+maxbin*(maxbin+
+c nest+1)). used as working space.
+c lwrk : integer. on entry,lwrk must specify the actual dimension of
+c the array wrk as declared in the calling (sub)program.lwrk
+c must not be too small (see wrk). unchanged on exit.
+c iwrk : integer array of dimension at least (maxtr*4+2*(maxbin+1))
+c used as working space.
+c kwrk : integer. on entry,kwrk must specify the actual dimension of
+c the array iwrk as declared in the calling (sub)program. kwrk
+c must not be too small (see iwrk). unchanged on exit.
+c ier : integer. error flag
+c ier=0 : normal return, s(x) satisfies the concavity/convexity
+c constraints and sq <= s.
+c ier<0 : abnormal termination: s(x) satisfies the concavity/
+c convexity constraints but sq > s.
+c ier=-3 : the requested storage space exceeds the available
+c storage space as specified by the parameter nest.
+c probably causes: nest too small. if nest is already
+c large (say nest > m/2), it may also indicate that s
+c is too small.
+c the approximation returned is the least-squares cubic
+c spline according to the knots t(1),...,t(n) (n=nest)
+c which satisfies the convexity constraints.
+c ier=-2 : the maximal number of knots n=m+4 has been reached.
+c probably causes: s too small.
+c ier=-1 : the number of knots n is less than the maximal number
+c m+4 but concon finds that adding one or more knots
+c will not further reduce the value of sq.
+c probably causes : s too small.
+c ier>0 : abnormal termination: no approximation is returned
+c ier=1 : the number of knots where s''(x)=0 exceeds maxbin.
+c probably causes : maxbin too small.
+c ier=2 : the number of records in the tree structure exceeds
+c maxtr.
+c probably causes : maxtr too small.
+c ier=3 : the algoritm finds no solution to the posed quadratic
+c programming problem.
+c probably causes : rounding errors.
+c ier=4 : the minimum number of knots (given by n) to guarantee
+c that the concavity/convexity conditions will be
+c satisfied is greater than nest.
+c probably causes: nest too small.
+c ier=5 : the minimum number of knots (given by n) to guarantee
+c that the concavity/convexity conditions will be
+c satisfied is greater than m+4.
+c probably causes: strongly alternating convexity and
+c concavity conditions. normally the situation can be
+c coped with by adding n-m-4 extra data points (found
+c by linear interpolation e.g.) with a small weight w(i)
+c and a v(i) number equal to zero.
+c ier=10 : on entry, the input data are controlled on validity.
+c the following restrictions must be satisfied
+c 0<=iopt<=1, m>3, nest>=8, s>=0, maxtr>=1, maxbin>=1,
+c kwrk>=maxtr*4+2*(maxbin+1), w(i)>0, x(i) < x(i+1),
+c lwrk>=m*4+nest*8+maxbin*(maxbin+nest+1)
+c if one of these restrictions is found to be violated
+c control is immediately repassed to the calling program
+c
+c further comments:
+c as an example of the use of the computation mode iopt=1, the
+c following program segment will cause concon to return control
+c each time a spline with a new set of knots has been computed.
+c .............
+c iopt = 0
+c s = 0.1e+60 (s very large)
+c do 10 i=1,m
+c call concon(iopt,m,x,y,w,v,s,nest,maxtr,maxbin,n,t,c,sq,sx,
+c * bind,wrk,lwrk,iwrk,kwrk,ier)
+c ......
+c s = sq
+c iopt=1
+c 10 continue
+c .............
+c
+c other subroutines required:
+c fpcoco,fpcosp,fpbspl,fpadno,fpdeno,fpseno,fpfrno
+c
+c references:
+c dierckx p. : an algorithm for cubic spline fitting with convexity
+c constraints, computing 24 (1980) 349-371.
+c dierckx p. : an algorithm for least-squares cubic spline fitting
+c with convexity and concavity constraints, report tw39,
+c dept. computer science, k.u.leuven, 1978.
+c dierckx p. : curve and surface fitting with splines, monographs on
+c numerical analysis, oxford university press, 1993.
+c
+c author:
+c p. dierckx
+c dept. computer science, k.u.leuven
+c celestijnenlaan 200a, b-3001 heverlee, belgium.
+c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
+c
+c creation date : march 1978
+c latest update : march 1987.
+c
+c ..
+c ..scalar arguments..
+ real*8 s,sq
+ integer iopt,m,nest,maxtr,maxbin,n,lwrk,kwrk,ier
+c ..array arguments..
+ real*8 x(m),y(m),w(m),v(m),t(nest),c(nest),sx(m),wrk(lwrk)
+ integer iwrk(kwrk)
+ logical bind(nest)
+c ..local scalars..
+ integer i,lwest,kwest,ie,iw,lww
+ real*8 one
+c ..
+c set constant
+ one = 0.1e+01
+c before starting computations a data check is made. if the input data
+c are invalid, control is immediately repassed to the calling program.
+ ier = 10
+ if(iopt.lt.0 .or. iopt.gt.1) go to 30
+ if(m.lt.4 .or. nest.lt.8) go to 30
+ if(s.lt.0.) go to 30
+ if(maxtr.lt.1 .or. maxbin.lt.1) go to 30
+ lwest = 8*nest+m*4+maxbin*(1+nest+maxbin)
+ kwest = 4*maxtr+2*(maxbin+1)
+ if(lwrk.lt.lwest .or. kwrk.lt.kwest) go to 30
+ if(iopt.gt.0) go to 20
+ if(w(1).le.0.) go to 30
+ if(v(1).gt.0.) v(1) = one
+ if(v(1).lt.0.) v(1) = -one
+ do 10 i=2,m
+ if(x(i-1).ge.x(i) .or. w(i).le.0.) go to 30
+ if(v(i).gt.0.) v(i) = one
+ if(v(i).lt.0.) v(i) = -one
+ 10 continue
+ 20 ier = 0
+c we partition the working space and determine the spline approximation
+ ie = 1
+ iw = ie+nest
+ lww = lwrk-nest
+ call fpcoco(iopt,m,x,y,w,v,s,nest,maxtr,maxbin,n,t,c,sq,sx,
+ * bind,wrk(ie),wrk(iw),lww,iwrk,kwrk,ier)
+ 30 return
+ end
Added: branches/Interpolate1D/fitpack/concur.f
===================================================================
--- branches/Interpolate1D/fitpack/concur.f 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/fitpack/concur.f 2008-07-18 19:44:12 UTC (rev 4550)
@@ -0,0 +1,370 @@
+ subroutine concur(iopt,idim,m,u,mx,x,xx,w,ib,db,nb,ie,de,ne,k,s,
+ * nest,n,t,nc,c,np,cp,fp,wrk,lwrk,iwrk,ier)
+c given the ordered set of m points x(i) in the idim-dimensional space
+c and given also a corresponding set of strictly increasing values u(i)
+c and the set of positive numbers w(i),i=1,2,...,m, subroutine concur
+c determines a smooth approximating spline curve s(u), i.e.
+c x1 = s1(u)
+c x2 = s2(u) ub = u(1) <= u <= u(m) = ue
+c .........
+c xidim = sidim(u)
+c with sj(u),j=1,2,...,idim spline functions of odd degree k with
+c common knots t(j),j=1,2,...,n.
+c in addition these splines will satisfy the following boundary
+c constraints (l)
+c if ib > 0 : sj (u(1)) = db(idim*l+j) ,l=0,1,...,ib-1
+c and (l)
+c if ie > 0 : sj (u(m)) = de(idim*l+j) ,l=0,1,...,ie-1.
+c if iopt=-1 concur calculates the weighted least-squares spline curve
+c according to a given set of knots.
+c if iopt>=0 the number of knots of the splines sj(u) and the position
+c t(j),j=1,2,...,n is chosen automatically by the routine. the smooth-
+c ness of s(u) is then achieved by minimalizing the discontinuity
+c jumps of the k-th derivative of s(u) at the knots t(j),j=k+2,k+3,...,
+c n-k-1. the amount of smoothness is determined by the condition that
+c f(p)=sum((w(i)*dist(x(i),s(u(i))))**2) be <= s, with s a given non-
+c negative constant, called the smoothing factor.
+c the fit s(u) is given in the b-spline representation and can be
+c evaluated by means of subroutine curev.
+c
+c calling sequence:
+c call concur(iopt,idim,m,u,mx,x,xx,w,ib,db,nb,ie,de,ne,k,s,nest,n,
+c * t,nc,c,np,cp,fp,wrk,lwrk,iwrk,ier)
+c
+c parameters:
+c iopt : integer flag. on entry iopt must specify whether a weighted
+c least-squares spline curve (iopt=-1) or a smoothing spline
+c curve (iopt=0 or 1) must be determined.if iopt=0 the routine
+c will start with an initial set of knots t(i)=ub,t(i+k+1)=ue,
+c i=1,2,...,k+1. if iopt=1 the routine will continue with the
+c knots found at the last call of the routine.
+c attention: a call with iopt=1 must always be immediately
+c preceded by another call with iopt=1 or iopt=0.
+c unchanged on exit.
+c idim : integer. on entry idim must specify the dimension of the
+c curve. 0 < idim < 11.
+c unchanged on exit.
+c m : integer. on entry m must specify the number of data points.
+c m > k-max(ib-1,0)-max(ie-1,0). unchanged on exit.
+c u : real array of dimension at least (m). before entry,
+c u(i) must be set to the i-th value of the parameter variable
+c u for i=1,2,...,m. these values must be supplied in
+c strictly ascending order and will be unchanged on exit.
+c mx : integer. on entry mx must specify the actual dimension of
+c the arrays x and xx as declared in the calling (sub)program
+c mx must not be too small (see x). unchanged on exit.
+c x : real array of dimension at least idim*m.
+c before entry, x(idim*(i-1)+j) must contain the j-th coord-
+c inate of the i-th data point for i=1,2,...,m and j=1,2,...,
+c idim. unchanged on exit.
+c xx : real array of dimension at least idim*m.
+c used as working space. on exit xx contains the coordinates
+c of the data points to which a spline curve with zero deriv-
+c ative constraints has been determined.
+c if the computation mode iopt =1 is used xx should be left
+c unchanged between calls.
+c w : real array of dimension at least (m). before entry, w(i)
+c must be set to the i-th value in the set of weights. the
+c w(i) must be strictly positive. unchanged on exit.
+c see also further comments.
+c ib : integer. on entry ib must specify the number of derivative
+c constraints for the curve at the begin point. 0<=ib<=(k+1)/2
+c unchanged on exit.
+c db : real array of dimension nb. before entry db(idim*l+j) must
+c contain the l-th order derivative of sj(u) at u=u(1) for
+c j=1,2,...,idim and l=0,1,...,ib-1 (if ib>0).
+c unchanged on exit.
+c nb : integer, specifying the dimension of db. nb>=max(1,idim*ib)
+c unchanged on exit.
+c ie : integer. on entry ie must specify the number of derivative
+c constraints for the curve at the end point. 0<=ie<=(k+1)/2
+c unchanged on exit.
+c de : real array of dimension ne. before entry de(idim*l+j) must
+c contain the l-th order derivative of sj(u) at u=u(m) for
+c j=1,2,...,idim and l=0,1,...,ie-1 (if ie>0).
+c unchanged on exit.
+c ne : integer, specifying the dimension of de. ne>=max(1,idim*ie)
+c unchanged on exit.
+c k : integer. on entry k must specify the degree of the splines.
+c k=1,3 or 5.
+c unchanged on exit.
+c s : real.on entry (in case iopt>=0) s must specify the smoothing
+c factor. s >=0. unchanged on exit.
+c for advice on the choice of s see further comments.
+c nest : integer. on entry nest must contain an over-estimate of the
+c total number of knots of the splines returned, to indicate
+c the storage space available to the routine. nest >=2*k+2.
+c in most practical situation nest=m/2 will be sufficient.
+c always large enough is nest=m+k+1+max(0,ib-1)+max(0,ie-1),
+c the number of knots needed for interpolation (s=0).
+c unchanged on exit.
+c n : integer.
+c unless ier = 10 (in case iopt >=0), n will contain the
+c total number of knots of the smoothing spline curve returned
+c if the computation mode iopt=1 is used this value of n
+c should be left unchanged between subsequent calls.
+c in case iopt=-1, the value of n must be specified on entry.
+c t : real array of dimension at least (nest).
+c on succesful exit, this array will contain the knots of the
+c spline curve,i.e. the position of the interior knots t(k+2),
+c t(k+3),..,t(n-k-1) as well as the position of the additional
+c t(1)=t(2)=...=t(k+1)=ub and t(n-k)=...=t(n)=ue needed for
+c the b-spline representation.
+c if the computation mode iopt=1 is used, the values of t(1),
+c t(2),...,t(n) should be left unchanged between subsequent
+c calls. if the computation mode iopt=-1 is used, the values
+c t(k+2),...,t(n-k-1) must be supplied by the user, before
+c entry. see also the restrictions (ier=10).
+c nc : integer. on entry nc must specify the actual dimension of
+c the array c as declared in the calling (sub)program. nc
+c must not be too small (see c). unchanged on exit.
+c c : real array of dimension at least (nest*idim).
+c on succesful exit, this array will contain the coefficients
+c in the b-spline representation of the spline curve s(u),i.e.
+c the b-spline coefficients of the spline sj(u) will be given
+c in c(n*(j-1)+i),i=1,2,...,n-k-1 for j=1,2,...,idim.
+c cp : real array of dimension at least 2*(k+1)*idim.
+c on exit cp will contain the b-spline coefficients of a
+c polynomial curve which satisfies the boundary constraints.
+c if the computation mode iopt =1 is used cp should be left
+c unchanged between calls.
+c np : integer. on entry np must specify the actual dimension of
+c the array cp as declared in the calling (sub)program. np
+c must not be too small (see cp). unchanged on exit.
+c fp : real. unless ier = 10, fp contains the weighted sum of
+c squared residuals of the spline curve returned.
+c wrk : real array of dimension at least m*(k+1)+nest*(6+idim+3*k).
+c used as working space. if the computation mode iopt=1 is
+c used, the values wrk(1),...,wrk(n) should be left unchanged
+c between subsequent calls.
+c lwrk : integer. on entry,lwrk must specify the actual dimension of
+c the array wrk as declared in the calling (sub)program. lwrk
+c must not be too small (see wrk). unchanged on exit.
+c iwrk : integer array of dimension at least (nest).
+c used as working space. if the computation mode iopt=1 is
+c used,the values iwrk(1),...,iwrk(n) should be left unchanged
+c between subsequent calls.
+c ier : integer. unless the routine detects an error, ier contains a
+c non-positive value on exit, i.e.
+c ier=0 : normal return. the curve returned has a residual sum of
+c squares fp such that abs(fp-s)/s <= tol with tol a relat-
+c ive tolerance set to 0.001 by the program.
+c ier=-1 : normal return. the curve returned is an interpolating
+c spline curve, satisfying the constraints (fp=0).
+c ier=-2 : normal return. the curve returned is the weighted least-
+c squares polynomial curve of degree k, satisfying the
+c constraints. in this extreme case fp gives the upper
+c bound fp0 for the smoothing factor s.
+c ier=1 : error. the required storage space exceeds the available
+c storage space, as specified by the parameter nest.
+c probably causes : nest too small. if nest is already
+c large (say nest > m/2), it may also indicate that s is
+c too small
+c the approximation returned is the least-squares spline
+c curve according to the knots t(1),t(2),...,t(n). (n=nest)
+c the parameter fp gives the corresponding weighted sum of
+c squared residuals (fp>s).
+c ier=2 : error. a theoretically impossible result was found during
+c the iteration proces for finding a smoothing spline curve
+c with fp = s. probably causes : s too small.
+c there is an approximation returned but the corresponding
+c weighted sum of squared residuals does not satisfy the
+c condition abs(fp-s)/s < tol.
+c ier=3 : error. the maximal number of iterations maxit (set to 20
+c by the program) allowed for finding a smoothing curve
+c with fp=s has been reached. probably causes : s too small
+c there is an approximation returned but the corresponding
+c weighted sum of squared residuals does not satisfy the
+c condition abs(fp-s)/s < tol.
+c ier=10 : error. on entry, the input data are controlled on validity
+c the following restrictions must be satisfied.
+c -1<=iopt<=1, k = 1,3 or 5, m>k-max(0,ib-1)-max(0,ie-1),
+c nest>=2k+2, 0<idim<=10, lwrk>=(k+1)*m+nest*(6+idim+3*k),
+c nc >=nest*idim ,u(1)<u(2)<...<u(m),w(i)>0 i=1,2,...,m,
+c mx>=idim*m,0<=ib<=(k+1)/2,0<=ie<=(k+1)/2,nb>=1,ne>=1,
+c nb>=ib*idim,ne>=ib*idim,np>=2*(k+1)*idim,
+c if iopt=-1:2*k+2<=n<=min(nest,mmax) with mmax = m+k+1+
+c max(0,ib-1)+max(0,ie-1)
+c u(1)<t(k+2)<t(k+3)<...<t(n-k-1)<u(m)
+c the schoenberg-whitney conditions, i.e. there
+c must be a subset of data points uu(j) such that
+c t(j) < uu(j) < t(j+k+1), j=1+max(0,ib-1),...
+c ,n+k-1-max(0,ie-1)
+c if iopt>=0: s>=0
+c if s=0 : nest >=mmax (see above)
+c if one of these conditions is found to be violated,control
+c is immediately repassed to the calling program. in that
+c case there is no approximation returned.
+c
+c further comments:
+c by means of the parameter s, the user can control the tradeoff
+c between closeness of fit and smoothness of fit of the approximation.
+c if s is too large, the curve will be too smooth and signal will be
+c lost ; if s is too small the curve will pick up too much noise. in
+c the extreme cases the program will return an interpolating curve if
+c s=0 and the least-squares polynomial curve of degree k if s is
+c very large. between these extremes, a properly chosen s will result
+c in a good compromise between closeness of fit and smoothness of fit.
+c to decide whether an approximation, corresponding to a certain s is
+c satisfactory the user is highly recommended to inspect the fits
+c graphically.
+c recommended values for s depend on the weights w(i). if these are
+c taken as 1/d(i) with d(i) an estimate of the standard deviation of
+c x(i), a good s-value should be found in the range (m-sqrt(2*m),m+
+c sqrt(2*m)). if nothing is known about the statistical error in x(i)
+c each w(i) can be set equal to one and s determined by trial and
+c error, taking account of the comments above. the best is then to
+c start with a very large value of s ( to determine the least-squares
+c polynomial curve and the upper bound fp0 for s) and then to
+c progressively decrease the value of s ( say by a factor 10 in the
+c beginning, i.e. s=fp0/10, fp0/100,...and more carefully as the
+c approximating curve shows more detail) to obtain closer fits.
+c to economize the search for a good s-value the program provides with
+c different modes of computation. at the first call of the routine, or
+c whenever he wants to restart with the initial set of knots the user
+c must set iopt=0.
+c if iopt=1 the program will continue with the set of knots found at
+c the last call of the routine. this will save a lot of computation
+c time if concur is called repeatedly for different values of s.
+c the number of knots of the spline returned and their location will
+c depend on the value of s and on the complexity of the shape of the
+c curve underlying the data. but, if the computation mode iopt=1 is
+c used, the knots returned may also depend on the s-values at previous
+c calls (if these were smaller). therefore, if after a number of
+c trials with different s-values and iopt=1, the user can finally
+c accept a fit as satisfactory, it may be worthwhile for him to call
+c concur once more with the selected value for s but now with iopt=0.
+c indeed, concur may then return an approximation of the same quality
+c of fit but with fewer knots and therefore better if data reduction
+c is also an important objective for the user.
+c
+c the form of the approximating curve can strongly be affected by
+c the choice of the parameter values u(i). if there is no physical
+c reason for choosing a particular parameter u, often good results
+c will be obtained with the choice
+c v(1)=0, v(i)=v(i-1)+q(i), i=2,...,m, u(i)=v(i)/v(m), i=1,..,m
+c where
+c q(i)= sqrt(sum j=1,idim (xj(i)-xj(i-1))**2 )
+c other possibilities for q(i) are
+c q(i)= sum j=1,idim (xj(i)-xj(i-1))**2
+c q(i)= sum j=1,idim abs(xj(i)-xj(i-1))
+c q(i)= max j=1,idim abs(xj(i)-xj(i-1))
+c q(i)= 1
+c
+c other subroutines required:
+c fpback,fpbspl,fpched,fpcons,fpdisc,fpgivs,fpknot,fprati,fprota
+c curev,fppocu,fpadpo,fpinst
+c
+c references:
+c dierckx p. : algorithms for smoothing data with periodic and
+c parametric splines, computer graphics and image
+c processing 20 (1982) 171-184.
+c dierckx p. : algorithms for smoothing data with periodic and param-
+c etric splines, report tw55, dept. computer science,
+c k.u.leuven, 1981.
+c dierckx p. : curve and surface fitting with splines, monographs on
+c numerical analysis, oxford university press, 1993.
+c
+c author:
+c p.dierckx
+c dept. computer science, k.u. leuven
+c celestijnenlaan 200a, b-3001 heverlee, belgium.
+c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
+c
+c creation date : may 1979
+c latest update : march 1987
+c
+c ..
+c ..scalar arguments..
+ real*8 s,fp
+ integer iopt,idim,m,mx,ib,nb,ie,ne,k,nest,n,nc,np,lwrk,ier
+c ..array arguments..
+ real*8 u(m),x(mx),xx(mx),db(nb),de(ne),w(m),t(nest),c(nc),wrk(lwrk
+ *)
+ real*8 cp(np)
+ integer iwrk(nest)
+c ..local scalars..
+ real*8 tol,dist
+ integer i,ib1,ie1,ja,jb,jfp,jg,jq,jz,j,k1,k2,lwest,maxit,nmin,
+ * ncc,kk,mmin,nmax,mxx
+c ..function references
+ integer max0
+c ..
+c we set up the parameters tol and maxit
+ maxit = 20
+ tol = 0.1e-02
+c before starting computations a data check is made. if the input data
+c are invalid, control is immediately repassed to the calling program.
+ ier = 10
+ if(iopt.lt.(-1) .or. iopt.gt.1) go to 90
+ if(idim.le.0 .or. idim.gt.10) go to 90
+ if(k.le.0 .or. k.gt.5) go to 90
+ k1 = k+1
+ kk = k1/2
+ if(kk*2.ne.k1) go to 90
+ k2 = k1+1
+ if(ib.lt.0 .or. ib.gt.kk) go to 90
+ if(ie.lt.0 .or. ie.gt.kk) go to 90
+ nmin = 2*k1
+ ib1 = max0(0,ib-1)
+ ie1 = max0(0,ie-1)
+ mmin = k1-ib1-ie1
+ if(m.lt.mmin .or. nest.lt.nmin) go to 90
+ if(nb.lt.(idim*ib) .or. ne.lt.(idim*ie)) go to 90
+ if(np.lt.(2*k1*idim)) go to 90
+ mxx = m*idim
+ ncc = nest*idim
+ if(mx.lt.mxx .or. nc.lt.ncc) go to 90
+ lwest = m*k1+nest*(6+idim+3*k)
+ if(lwrk.lt.lwest) go to 90
+ if(w(1).le.0.) go to 90
+ do 10 i=2,m
+ if(u(i-1).ge.u(i) .or. w(i).le.0.) go to 90
+ 10 continue
+ if(iopt.ge.0) go to 30
+ if(n.lt.nmin .or. n.gt.nest) go to 90
+ j = n
+ do 20 i=1,k1
+ t(i) = u(1)
+ t(j) = u(m)
+ j = j-1
+ 20 continue
+ call fpched(u,m,t,n,k,ib,ie,ier)
+ if (ier.eq.0) go to 40
+ go to 90
+ 30 if(s.lt.0.) go to 90
+ nmax = m+k1+ib1+ie1
+ if(s.eq.0. .and. nest.lt.nmax) go to 90
+ ier = 0
+ if(iopt.gt.0) go to 70
+c we determine a polynomial curve satisfying the boundary constraints.
+ 40 call fppocu(idim,k,u(1),u(m),ib,db,nb,ie,de,ne,cp,np)
+c we generate new data points which will be approximated by a spline
+c with zero derivative constraints.
+ j = nmin
+ do 50 i=1,k1
+ wrk(i) = u(1)
+ wrk(j) = u(m)
+ j = j-1
+ 50 continue
+c evaluate the polynomial curve
+ call curev(idim,wrk,nmin,cp,np,k,u,m,xx,mxx,ier)
+c substract from the old data, the values of the polynomial curve
+ do 60 i=1,mxx
+ xx(i) = x(i)-xx(i)
+ 60 continue
+c we partition the working space and determine the spline curve.
+ 70 jfp = 1
+ jz = jfp+nest
+ ja = jz+ncc
+ jb = ja+nest*k1
+ jg = jb+nest*k2
+ jq = jg+nest*k2
+ call fpcons(iopt,idim,m,u,mxx,xx,w,ib,ie,k,s,nest,tol,maxit,k1,
+ * k2,n,t,ncc,c,fp,wrk(jfp),wrk(jz),wrk(ja),wrk(jb),wrk(jg),wrk(jq),
+ *
+ * iwrk,ier)
+c add the polynomial curve to the calculated spline.
+ call fpadpo(idim,t,n,c,ncc,k,cp,np,wrk(jz),wrk(ja),wrk(jb))
+ 90 return
+ end
Added: branches/Interpolate1D/fitpack/cualde.f
===================================================================
--- branches/Interpolate1D/fitpack/cualde.f 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/fitpack/cualde.f 2008-07-18 19:44:12 UTC (rev 4550)
@@ -0,0 +1,91 @@
+ subroutine cualde(idim,t,n,c,nc,k1,u,d,nd,ier)
+c subroutine cualde evaluates at the point u all the derivatives
+c (l)
+c d(idim*l+j) = sj (u) ,l=0,1,...,k, j=1,2,...,idim
+c of a spline curve s(u) of order k1 (degree k=k1-1) and dimension idim
+c given in its b-spline representation.
+c
+c calling sequence:
+c call cualde(idim,t,n,c,nc,k1,u,d,nd,ier)
+c
+c input parameters:
+c idim : integer, giving the dimension of the spline curve.
+c t : array,length n, which contains the position of the knots.
+c n : integer, giving the total number of knots of s(u).
+c c : array,length nc, which contains the b-spline coefficients.
+c nc : integer, giving the total number of coefficients of s(u).
+c k1 : integer, giving the order of s(u) (order=degree+1).
+c u : real, which contains the point where the derivatives must
+c be evaluated.
+c nd : integer, giving the dimension of the array d. nd >= k1*idim
+c
+c output parameters:
+c d : array,length nd,giving the different curve derivatives.
+c d(idim*l+j) will contain the j-th coordinate of the l-th
+c derivative of the curve at the point u.
+c ier : error flag
+c ier = 0 : normal return
+c ier =10 : invalid input data (see restrictions)
+c
+c restrictions:
+c nd >= k1*idim
+c t(k1) <= u <= t(n-k1+1)
+c
+c further comments:
+c if u coincides with a knot, right derivatives are computed
+c ( left derivatives if u = t(n-k1+1) ).
+c
+c other subroutines required: fpader.
+c
+c references :
+c de boor c : on calculating with b-splines, j. approximation theory
+c 6 (1972) 50-62.
+c cox m.g. : the numerical evaluation of b-splines, j. inst. maths
+c applics 10 (1972) 134-149.
+c dierckx p. : curve and surface fitting with splines, monographs on
+c numerical analysis, oxford university press, 1993.
+c
+c author :
+c p.dierckx
+c dept. computer science, k.u.leuven
+c celestijnenlaan 200a, b-3001 heverlee, belgium.
+c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
+c
+c latest update : march 1987
+c
+c ..scalar arguments..
+ integer idim,n,nc,k1,nd,ier
+ real*8 u
+c ..array arguments..
+ real*8 t(n),c(nc),d(nd)
+c ..local scalars..
+ integer i,j,kk,l,m,nk1
+c ..local array..
+ real*8 h(6)
+c ..
+c before starting computations a data check is made. if the input data
+c are invalid control is immediately repassed to the calling program.
+ ier = 10
+ if(nd.lt.(k1*idim)) go to 500
+ nk1 = n-k1
+ if(u.lt.t(k1) .or. u.gt.t(nk1+1)) go to 500
+c search for knot interval t(l) <= u < t(l+1)
+ l = k1
+ 100 if(u.lt.t(l+1) .or. l.eq.nk1) go to 200
+ l = l+1
+ go to 100
+ 200 if(t(l).ge.t(l+1)) go to 500
+ ier = 0
+c calculate the derivatives.
+ j = 1
+ do 400 i=1,idim
+ call fpader(t,n,c(j),k1,u,l,h)
+ m = i
+ do 300 kk=1,k1
+ d(m) = h(kk)
+ m = m+idim
+ 300 continue
+ j = j+n
+ 400 continue
+ 500 return
+ end
Added: branches/Interpolate1D/fitpack/curev.f
===================================================================
--- branches/Interpolate1D/fitpack/curev.f 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/fitpack/curev.f 2008-07-18 19:44:12 UTC (rev 4550)
@@ -0,0 +1,110 @@
+ subroutine curev(idim,t,n,c,nc,k,u,m,x,mx,ier)
+c subroutine curev evaluates in a number of points u(i),i=1,2,...,m
+c a spline curve s(u) of degree k and dimension idim, given in its
+c b-spline representation.
+c
+c calling sequence:
+c call curev(idim,t,n,c,nc,k,u,m,x,mx,ier)
+c
+c input parameters:
+c idim : integer, giving the dimension of the spline curve.
+c t : array,length n, which contains the position of the knots.
+c n : integer, giving the total number of knots of s(u).
+c c : array,length nc, which contains the b-spline coefficients.
+c nc : integer, giving the total number of coefficients of s(u).
+c k : integer, giving the degree of s(u).
+c u : array,length m, which contains the points where s(u) must
+c be evaluated.
+c m : integer, giving the number of points where s(u) must be
+c evaluated.
+c mx : integer, giving the dimension of the array x. mx >= m*idim
+c
+c output parameters:
+c x : array,length mx,giving the value of s(u) at the different
+c points. x(idim*(i-1)+j) will contain the j-th coordinate
+c of the i-th point on the curve.
+c ier : error flag
+c ier = 0 : normal return
+c ier =10 : invalid input data (see restrictions)
+c
+c restrictions:
+c m >= 1
+c mx >= m*idim
+c t(k+1) <= u(i) <= u(i+1) <= t(n-k) , i=1,2,...,m-1.
+c
+c other subroutines required: fpbspl.
+c
+c references :
+c de boor c : on calculating with b-splines, j. approximation theory
+c 6 (1972) 50-62.
+c cox m.g. : the numerical evaluation of b-splines, j. inst. maths
+c applics 10 (1972) 134-149.
+c dierckx p. : curve and surface fitting with splines, monographs on
+c numerical analysis, oxford university press, 1993.
+c
+c author :
+c p.dierckx
+c dept. computer science, k.u.leuven
+c celestijnenlaan 200a, b-3001 heverlee, belgium.
+c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
+c
+c latest update : march 1987
+c
+c ..scalar arguments..
+ integer idim,n,nc,k,m,mx,ier
+c ..array arguments..
+ real*8 t(n),c(nc),u(m),x(mx)
+c ..local scalars..
+ integer i,j,jj,j1,k1,l,ll,l1,mm,nk1
+ real*8 arg,sp,tb,te
+c ..local array..
+ real*8 h(6)
+c ..
+c before starting computations a data check is made. if the input data
+c are invalid control is immediately repassed to the calling program.
+ ier = 10
+ if (m.lt.1) go to 100
+ if (m.eq.1) go to 30
+ go to 10
+ 10 do 20 i=2,m
+ if(u(i).lt.u(i-1)) go to 100
+ 20 continue
+ 30 if(mx.lt.(m*idim)) go to 100
+ ier = 0
+c fetch tb and te, the boundaries of the approximation interval.
+ k1 = k+1
+ nk1 = n-k1
+ tb = t(k1)
+ te = t(nk1+1)
+ l = k1
+ l1 = l+1
+c main loop for the different points.
+ mm = 0
+ do 80 i=1,m
+c fetch a new u-value arg.
+ arg = u(i)
+ if(arg.lt.tb) arg = tb
+ if(arg.gt.te) arg = te
+c search for knot interval t(l) <= arg < t(l+1)
+ 40 if(arg.lt.t(l1) .or. l.eq.nk1) go to 50
+ l = l1
+ l1 = l+1
+ go to 40
+c evaluate the non-zero b-splines at arg.
+ 50 call fpbspl(t,n,k,arg,l,h)
+c find the value of s(u) at u=arg.
+ ll = l-k1
+ do 70 j1=1,idim
+ jj = ll
+ sp = 0.
+ do 60 j=1,k1
+ jj = jj+1
+ sp = sp+c(jj)*h(j)
+ 60 continue
+ mm = mm+1
+ x(mm) = sp
+ ll = ll+n
+ 70 continue
+ 80 continue
+ 100 return
+ end
Added: branches/Interpolate1D/fitpack/curfit.f
===================================================================
--- branches/Interpolate1D/fitpack/curfit.f 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/fitpack/curfit.f 2008-07-18 19:44:12 UTC (rev 4550)
@@ -0,0 +1,261 @@
+ subroutine curfit(iopt,m,x,y,w,xb,xe,k,s,nest,n,t,c,fp,
+ * wrk,lwrk,iwrk,ier)
+c given the set of data points (x(i),y(i)) and the set of positive
+c numbers w(i),i=1,2,...,m,subroutine curfit determines a smooth spline
+c approximation of degree k on the interval xb <= x <= xe.
+c if iopt=-1 curfit calculates the weighted least-squares spline
+c according to a given set of knots.
+c if iopt>=0 the number of knots of the spline s(x) and the position
+c t(j),j=1,2,...,n is chosen automatically by the routine. the smooth-
+c ness of s(x) is then achieved by minimalizing the discontinuity
+c jumps of the k-th derivative of s(x) at the knots t(j),j=k+2,k+3,...,
+c n-k-1. the amount of smoothness is determined by the condition that
+c f(p)=sum((w(i)*(y(i)-s(x(i))))**2) be <= s, with s a given non-
+c negative constant, called the smoothing factor.
+c the fit s(x) is given in the b-spline representation (b-spline coef-
+c ficients c(j),j=1,2,...,n-k-1) and can be evaluated by means of
+c subroutine splev.
+c
+c calling sequence:
+c call curfit(iopt,m,x,y,w,xb,xe,k,s,nest,n,t,c,fp,wrk,
+c * lwrk,iwrk,ier)
+c
+c parameters:
+c iopt : integer flag. on entry iopt must specify whether a weighted
+c least-squares spline (iopt=-1) or a smoothing spline (iopt=
+c 0 or 1) must be determined. if iopt=0 the routine will start
+c with an initial set of knots t(i)=xb, t(i+k+1)=xe, i=1,2,...
+c k+1. if iopt=1 the routine will continue with the knots
+c found at the last call of the routine.
+c attention: a call with iopt=1 must always be immediately
+c preceded by another call with iopt=1 or iopt=0.
+c unchanged on exit.
+c m : integer. on entry m must specify the number of data points.
+c m > k. unchanged on exit.
+c x : real array of dimension at least (m). before entry, x(i)
+c must be set to the i-th value of the independent variable x,
+c for i=1,2,...,m. these values must be supplied in strictly
+c ascending order. unchanged on exit.
+c y : real array of dimension at least (m). before entry, y(i)
+c must be set to the i-th value of the dependent variable y,
+c for i=1,2,...,m. unchanged on exit.
+c w : real array of dimension at least (m). before entry, w(i)
+c must be set to the i-th value in the set of weights. the
+c w(i) must be strictly positive. unchanged on exit.
+c see also further comments.
+c xb,xe : real values. on entry xb and xe must specify the boundaries
+c of the approximation interval. xb<=x(1), xe>=x(m).
+c unchanged on exit.
+c k : integer. on entry k must specify the degree of the spline.
+c 1<=k<=5. it is recommended to use cubic splines (k=3).
+c the user is strongly dissuaded from choosing k even,together
+c with a small s-value. unchanged on exit.
+c s : real.on entry (in case iopt>=0) s must specify the smoothing
+c factor. s >=0. unchanged on exit.
+c for advice on the choice of s see further comments.
+c nest : integer. on entry nest must contain an over-estimate of the
+c total number of knots of the spline returned, to indicate
+c the storage space available to the routine. nest >=2*k+2.
+c in most practical situation nest=m/2 will be sufficient.
+c always large enough is nest=m+k+1, the number of knots
+c needed for interpolation (s=0). unchanged on exit.
+c n : integer.
+c unless ier =10 (in case iopt >=0), n will contain the
+c total number of knots of the spline approximation returned.
+c if the computation mode iopt=1 is used this value of n
+c should be left unchanged between subsequent calls.
+c in case iopt=-1, the value of n must be specified on entry.
+c t : real array of dimension at least (nest).
+c on succesful exit, this array will contain the knots of the
+c spline,i.e. the position of the interior knots t(k+2),t(k+3)
+c ...,t(n-k-1) as well as the position of the additional knots
+c t(1)=t(2)=...=t(k+1)=xb and t(n-k)=...=t(n)=xe needed for
+c the b-spline representation.
+c if the computation mode iopt=1 is used, the values of t(1),
+c t(2),...,t(n) should be left unchanged between subsequent
+c calls. if the computation mode iopt=-1 is used, the values
+c t(k+2),...,t(n-k-1) must be supplied by the user, before
+c entry. see also the restrictions (ier=10).
+c c : real array of dimension at least (nest).
+c on succesful exit, this array will contain the coefficients
+c c(1),c(2),..,c(n-k-1) in the b-spline representation of s(x)
+c fp : real. unless ier=10, fp contains the weighted sum of
+c squared residuals of the spline approximation returned.
+c wrk : real array of dimension at least (m*(k+1)+nest*(7+3*k)).
+c used as working space. if the computation mode iopt=1 is
+c used, the values wrk(1),...,wrk(n) should be left unchanged
+c between subsequent calls.
+c lwrk : integer. on entry,lwrk must specify the actual dimension of
+c the array wrk as declared in the calling (sub)program.lwrk
+c must not be too small (see wrk). unchanged on exit.
+c iwrk : integer array of dimension at least (nest).
+c used as working space. if the computation mode iopt=1 is
+c used,the values iwrk(1),...,iwrk(n) should be left unchanged
+c between subsequent calls.
+c ier : integer. unless the routine detects an error, ier contains a
+c non-positive value on exit, i.e.
+c ier=0 : normal return. the spline returned has a residual sum of
+c squares fp such that abs(fp-s)/s <= tol with tol a relat-
+c ive tolerance set to 0.001 by the program.
+c ier=-1 : normal return. the spline returned is an interpolating
+c spline (fp=0).
+c ier=-2 : normal return. the spline returned is the weighted least-
+c squares polynomial of degree k. in this extreme case fp
+c gives the upper bound fp0 for the smoothing factor s.
+c ier=1 : error. the required storage space exceeds the available
+c storage space, as specified by the parameter nest.
+c probably causes : nest too small. if nest is already
+c large (say nest > m/2), it may also indicate that s is
+c too small
+c the approximation returned is the weighted least-squares
+c spline according to the knots t(1),t(2),...,t(n). (n=nest)
+c the parameter fp gives the corresponding weighted sum of
+c squared residuals (fp>s).
+c ier=2 : error. a theoretically impossible result was found during
+c the iteration proces for finding a smoothing spline with
+c fp = s. probably causes : s too small.
+c there is an approximation returned but the corresponding
+c weighted sum of squared residuals does not satisfy the
+c condition abs(fp-s)/s < tol.
+c ier=3 : error. the maximal number of iterations maxit (set to 20
+c by the program) allowed for finding a smoothing spline
+c with fp=s has been reached. probably causes : s too small
+c there is an approximation returned but the corresponding
+c weighted sum of squared residuals does not satisfy the
+c condition abs(fp-s)/s < tol.
+c ier=10 : error. on entry, the input data are controlled on validity
+c the following restrictions must be satisfied.
+c -1<=iopt<=1, 1<=k<=5, m>k, nest>2*k+2, w(i)>0,i=1,2,...,m
+c xb<=x(1)<x(2)<...<x(m)<=xe, lwrk>=(k+1)*m+nest*(7+3*k)
+c if iopt=-1: 2*k+2<=n<=min(nest,m+k+1)
+c xb<t(k+2)<t(k+3)<...<t(n-k-1)<xe
+c the schoenberg-whitney conditions, i.e. there
+c must be a subset of data points xx(j) such that
+c t(j) < xx(j) < t(j+k+1), j=1,2,...,n-k-1
+c if iopt>=0: s>=0
+c if s=0 : nest >= m+k+1
+c if one of these conditions is found to be violated,control
+c is immediately repassed to the calling program. in that
+c case there is no approximation returned.
+c
+c further comments:
+c by means of the parameter s, the user can control the tradeoff
+c between closeness of fit and smoothness of fit of the approximation.
+c if s is too large, the spline will be too smooth and signal will be
+c lost ; if s is too small the spline will pick up too much noise. in
+c the extreme cases the program will return an interpolating spline if
+c s=0 and the weighted least-squares polynomial of degree k if s is
+c very large. between these extremes, a properly chosen s will result
+c in a good compromise between closeness of fit and smoothness of fit.
+c to decide whether an approximation, corresponding to a certain s is
+c satisfactory the user is highly recommended to inspect the fits
+c graphically.
+c recommended values for s depend on the weights w(i). if these are
+c taken as 1/d(i) with d(i) an estimate of the standard deviation of
+c y(i), a good s-value should be found in the range (m-sqrt(2*m),m+
+c sqrt(2*m)). if nothing is known about the statistical error in y(i)
+c each w(i) can be set equal to one and s determined by trial and
+c error, taking account of the comments above. the best is then to
+c start with a very large value of s ( to determine the least-squares
+c polynomial and the corresponding upper bound fp0 for s) and then to
+c progressively decrease the value of s ( say by a factor 10 in the
+c beginning, i.e. s=fp0/10, fp0/100,...and more carefully as the
+c approximation shows more detail) to obtain closer fits.
+c to economize the search for a good s-value the program provides with
+c different modes of computation. at the first call of the routine, or
+c whenever he wants to restart with the initial set of knots the user
+c must set iopt=0.
+c if iopt=1 the program will continue with the set of knots found at
+c the last call of the routine. this will save a lot of computation
+c time if curfit is called repeatedly for different values of s.
+c the number of knots of the spline returned and their location will
+c depend on the value of s and on the complexity of the shape of the
+c function underlying the data. but, if the computation mode iopt=1
+c is used, the knots returned may also depend on the s-values at
+c previous calls (if these were smaller). therefore, if after a number
+c of trials with different s-values and iopt=1, the user can finally
+c accept a fit as satisfactory, it may be worthwhile for him to call
+c curfit once more with the selected value for s but now with iopt=0.
+c indeed, curfit may then return an approximation of the same quality
+c of fit but with fewer knots and therefore better if data reduction
+c is also an important objective for the user.
+c
+c other subroutines required:
+c fpback,fpbspl,fpchec,fpcurf,fpdisc,fpgivs,fpknot,fprati,fprota
+c
+c references:
+c dierckx p. : an algorithm for smoothing, differentiation and integ-
+c ration of experimental data using spline functions,
+c j.comp.appl.maths 1 (1975) 165-184.
+c dierckx p. : a fast algorithm for smoothing data on a rectangular
+c grid while using spline functions, siam j.numer.anal.
+c 19 (1982) 1286-1304.
+c dierckx p. : an improved algorithm for curve fitting with spline
+c functions, report tw54, dept. computer science,k.u.
+c leuven, 1981.
+c dierckx p. : curve and surface fitting with splines, monographs on
+c numerical analysis, oxford university press, 1993.
+c
+c author:
+c p.dierckx
+c dept. computer science, k.u. leuven
+c celestijnenlaan 200a, b-3001 heverlee, belgium.
+c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
+c
+c creation date : may 1979
+c latest update : march 1987
+c
+c ..
+c ..scalar arguments..
+ real*8 xb,xe,s,fp
+ integer iopt,m,k,nest,n,lwrk,ier
+c ..array arguments..
+ real*8 x(m),y(m),w(m),t(nest),c(nest),wrk(lwrk)
+ integer iwrk(nest)
+c ..local scalars..
+ real*8 tol
+ integer i,ia,ib,ifp,ig,iq,iz,j,k1,k2,lwest,maxit,nmin
+c ..
+c we set up the parameters tol and maxit
+ maxit = 20
+ tol = 0.1d-02
+c before starting computations a data check is made. if the input data
+c are invalid, control is immediately repassed to the calling program.
+ ier = 10
+ if(k.le.0 .or. k.gt.5) go to 50
+ k1 = k+1
+ k2 = k1+1
+ if(iopt.lt.(-1) .or. iopt.gt.1) go to 50
+ nmin = 2*k1
+ if(m.lt.k1 .or. nest.lt.nmin) go to 50
+ lwest = m*k1+nest*(7+3*k)
+ if(lwrk.lt.lwest) go to 50
+ if(xb.gt.x(1) .or. xe.lt.x(m) .or. w(1).le.0.) go to 50
+ do 10 i=2,m
+ if(x(i-1).ge.x(i) .or. w(i).le.0.) go to 50
+ 10 continue
+ if(iopt.ge.0) go to 30
+ if(n.lt.nmin .or. n.gt.nest) go to 50
+ j = n
+ do 20 i=1,k1
+ t(i) = xb
+ t(j) = xe
+ j = j-1
+ 20 continue
+ call fpchec(x,m,t,n,k,ier)
+ if (ier.eq.0) go to 40
+ go to 50
+ 30 if(s.lt.0.) go to 50
+ if(s.eq.0. .and. nest.lt.(m+k1)) go to 50
+ ier = 0
+c we partition the working space and determine the spline approximation.
+ 40 ifp = 1
+ iz = ifp+nest
+ ia = iz+nest
+ ib = ia+nest*k1
+ ig = ib+nest*k2
+ iq = ig+nest*k2
+ call fpcurf(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,c,fp,
+ * wrk(ifp),wrk(iz),wrk(ia),wrk(ib),wrk(ig),wrk(iq),iwrk,ier)
+ 50 return
+ end
Added: branches/Interpolate1D/fitpack/dblint.f
===================================================================
--- branches/Interpolate1D/fitpack/dblint.f 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/fitpack/dblint.f 2008-07-18 19:44:12 UTC (rev 4550)
@@ -0,0 +1,88 @@
+ real*8 function dblint(tx,nx,ty,ny,c,kx,ky,xb,xe,yb,ye,wrk)
+c function dblint calculates the double integral
+c / xe / ye
+c | | s(x,y) dx dy
+c xb / yb /
+c with s(x,y) a bivariate spline of degrees kx and ky, given in the
+c b-spline representation.
+c
+c calling sequence:
+c aint = dblint(tx,nx,ty,ny,c,kx,ky,xb,xe,yb,ye,wrk)
+c
+c input parameters:
+c tx : real array, length nx, which contains the position of the
+c knots in the x-direction.
+c nx : integer, giving the total number of knots in the x-direction
+c ty : real array, length ny, which contains the position of the
+c knots in the y-direction.
+c ny : integer, giving the total number of knots in the y-direction
+c c : real array, length (nx-kx-1)*(ny-ky-1), which contains the
+c b-spline coefficients.
+c kx,ky : integer values, giving the degrees of the spline.
+c xb,xe : real values, containing the boundaries of the integration
+c yb,ye domain. s(x,y) is considered to be identically zero out-
+c side the rectangle (tx(kx+1),tx(nx-kx))*(ty(ky+1),ty(ny-ky))
+c
+c output parameters:
+c aint : real , containing the double integral of s(x,y).
+c wrk : real array of dimension at least (nx+ny-kx-ky-2).
+c used as working space.
+c on exit, wrk(i) will contain the integral
+c / xe
+c | ni,kx+1(x) dx , i=1,2,...,nx-kx-1
+c xb /
+c with ni,kx+1(x) the normalized b-spline defined on
+c the knots tx(i),...,tx(i+kx+1)
+c wrk(j+nx-kx-1) will contain the integral
+c / ye
+c | nj,ky+1(y) dy , j=1,2,...,ny-ky-1
+c yb /
+c with nj,ky+1(y) the normalized b-spline defined on
+c the knots ty(j),...,ty(j+ky+1)
+c
+c other subroutines required: fpintb
+c
+c references :
+c gaffney p.w. : the calculation of indefinite integrals of b-splines
+c j. inst. maths applics 17 (1976) 37-41.
+c dierckx p. : curve and surface fitting with splines, monographs on
+c numerical analysis, oxford university press, 1993.
+c
+c author :
+c p.dierckx
+c dept. computer science, k.u.leuven
+c celestijnenlaan 200a, b-3001 heverlee, belgium.
+c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
+c
+c latest update : march 1989
+c
+c ..scalar arguments..
+ integer nx,ny,kx,ky
+ real*8 xb,xe,yb,ye
+c ..array arguments..
+ real*8 tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1)),wrk(nx+ny-kx-ky-2)
+c ..local scalars..
+ integer i,j,l,m,nkx1,nky1
+ real*8 res
+c ..
+ nkx1 = nx-kx-1
+ nky1 = ny-ky-1
+c we calculate the integrals of the normalized b-splines ni,kx+1(x)
+ call fpintb(tx,nx,wrk,nkx1,xb,xe)
+c we calculate the integrals of the normalized b-splines nj,ky+1(y)
+ call fpintb(ty,ny,wrk(nkx1+1),nky1,yb,ye)
+c calculate the integral of s(x,y)
+ dblint = 0.
+ do 200 i=1,nkx1
+ res = wrk(i)
+ if(res.eq.0.) go to 200
+ m = (i-1)*nky1
+ l = nkx1
+ do 100 j=1,nky1
+ m = m+1
+ l = l+1
+ dblint = dblint+res*wrk(l)*c(m)
+ 100 continue
+ 200 continue
+ return
+ end
Added: branches/Interpolate1D/fitpack/evapol.f
===================================================================
--- branches/Interpolate1D/fitpack/evapol.f 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/fitpack/evapol.f 2008-07-18 19:44:12 UTC (rev 4550)
@@ -0,0 +1,82 @@
+ real*8 function evapol(tu,nu,tv,nv,c,rad,x,y)
+c function program evacir evaluates the function f(x,y) = s(u,v),
+c defined through the transformation
+c x = u*rad(v)*cos(v) y = u*rad(v)*sin(v)
+c and where s(u,v) is a bicubic spline ( 0<=u<=1 , -pi<=v<=pi ), given
+c in its standard b-spline representation.
+c
+c calling sequence:
+c f = evapol(tu,nu,tv,nv,c,rad,x,y)
+c
+c input parameters:
+c tu : real array, length nu, which contains the position of the
+c knots in the u-direction.
+c nu : integer, giving the total number of knots in the u-direction
+c tv : real array, length nv, which contains the position of the
+c knots in the v-direction.
+c nv : integer, giving the total number of knots in the v-direction
+c c : real array, length (nu-4)*(nv-4), which contains the
+c b-spline coefficients.
+c rad : real function subprogram, defining the boundary of the
+c approximation domain. must be declared external in the
+c calling (sub)-program
+c x,y : real values.
+c before entry x and y must be set to the co-ordinates of
+c the point where f(x,y) must be evaluated.
+c
+c output parameter:
+c f : real
+c on exit f contains the value of f(x,y)
+c
+c other subroutines required:
+c bispev,fpbisp,fpbspl
+c
+c references :
+c de boor c : on calculating with b-splines, j. approximation theory
+c 6 (1972) 50-62.
+c cox m.g. : the numerical evaluation of b-splines, j. inst. maths
+c applics 10 (1972) 134-149.
+c dierckx p. : curve and surface fitting with splines, monographs on
+c numerical analysis, oxford university press, 1993.
+c
+c author :
+c p.dierckx
+c dept. computer science, k.u.leuven
+c celestijnenlaan 200a, b-3001 heverlee, belgium.
+c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
+c
+c latest update : march 1989
+c
+c ..scalar arguments..
+ integer nu,nv
+ real*8 x,y
+c ..array arguments..
+ real*8 tu(nu),tv(nv),c((nu-4)*(nv-4))
+c ..user specified function
+ real*8 rad
+c ..local scalars..
+ integer ier
+ real*8 u,v,r,f,one,dist
+c ..local arrays
+ real*8 wrk(8)
+ integer iwrk(2)
+c ..function references
+ real*8 atan2,sqrt
+c ..
+c calculate the (u,v)-coordinates of the given point.
+ one = 1
+ u = 0.
+ v = 0.
+ dist = x**2+y**2
+ if(dist.le.0.) go to 10
+ v = atan2(y,x)
+ r = rad(v)
+ if(r.le.0.) go to 10
+ u = sqrt(dist)/r
+ if(u.gt.one) u = one
+c evaluate s(u,v)
+ 10 call bispev(tu,nu,tv,nv,c,3,3,u,1,v,1,f,wrk,8,iwrk,2,ier)
+ evapol = f
+ return
+ end
+
Added: branches/Interpolate1D/fitpack/fourco.f
===================================================================
--- branches/Interpolate1D/fitpack/fourco.f 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/fitpack/fourco.f 2008-07-18 19:44:12 UTC (rev 4550)
@@ -0,0 +1,96 @@
+ subroutine fourco(t,n,c,alfa,m,ress,resc,wrk1,wrk2,ier)
+c subroutine fourco calculates the integrals
+c /t(n-3)
+c ress(i) = ! s(x)*sin(alfa(i)*x) dx and
+c t(4)/
+c /t(n-3)
+c resc(i) = ! s(x)*cos(alfa(i)*x) dx, i=1,...,m,
+c t(4)/
+c where s(x) denotes a cubic spline which is given in its
+c b-spline representation.
+c
+c calling sequence:
+c call fourco(t,n,c,alfa,m,ress,resc,wrk1,wrk2,ier)
+c
+c input parameters:
+c t : real array,length n, containing the knots of s(x).
+c n : integer, containing the total number of knots. n>=10.
+c c : real array,length n, containing the b-spline coefficients.
+c alfa : real array,length m, containing the parameters alfa(i).
+c m : integer, specifying the number of integrals to be computed.
+c wrk1 : real array,length n. used as working space
+c wrk2 : real array,length n. used as working space
+c
+c output parameters:
+c ress : real array,length m, containing the integrals ress(i).
+c resc : real array,length m, containing the integrals resc(i).
+c ier : error flag:
+c ier=0 : normal return.
+c ier=10: invalid input data (see restrictions).
+c
+c restrictions:
+c n >= 10
+c t(4) < t(5) < ... < t(n-4) < t(n-3).
+c t(1) <= t(2) <= t(3) <= t(4).
+c t(n-3) <= t(n-2) <= t(n-1) <= t(n).
+c
+c other subroutines required: fpbfou,fpcsin
+c
+c references :
+c dierckx p. : calculation of fouriercoefficients of discrete
+c functions using cubic splines. j. computational
+c and applied mathematics 3 (1977) 207-209.
+c dierckx p. : curve and surface fitting with splines, monographs on
+c numerical analysis, oxford university press, 1993.
+c
+c author :
+c p.dierckx
+c dept. computer science, k.u.leuven
+c celestijnenlaan 200a, b-3001 heverlee, belgium.
+c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
+c
+c latest update : march 1987
+c
+c ..scalar arguments..
+ integer n,m,ier
+c ..array arguments..
+ real*8 t(n),c(n),wrk1(n),wrk2(n),alfa(m),ress(m),resc(m)
+c ..local scalars..
+ integer i,j,n4
+ real*8 rs,rc
+c ..
+ n4 = n-4
+c before starting computations a data check is made. in the input data
+c are invalid, control is immediately repassed to the calling program.
+ ier = 10
+ if(n.lt.10) go to 50
+ j = n
+ do 10 i=1,3
+ if(t(i).gt.t(i+1)) go to 50
+ if(t(j).lt.t(j-1)) go to 50
+ j = j-1
+ 10 continue
+ do 20 i=4,n4
+ if(t(i).ge.t(i+1)) go to 50
+ 20 continue
+ ier = 0
+c main loop for the different alfa(i).
+ do 40 i=1,m
+c calculate the integrals
+c wrk1(j) = integral(nj,4(x)*sin(alfa*x)) and
+c wrk2(j) = integral(nj,4(x)*cos(alfa*x)), j=1,2,...,n-4,
+c where nj,4(x) denotes the normalised cubic b-spline defined on the
+c knots t(j),t(j+1),...,t(j+4).
+ call fpbfou(t,n,alfa(i),wrk1,wrk2)
+c calculate the integrals ress(i) and resc(i).
+ rs = 0.
+ rc = 0.
+ do 30 j=1,n4
+ rs = rs+c(j)*wrk1(j)
+ rc = rc+c(j)*wrk2(j)
+ 30 continue
+ ress(i) = rs
+ resc(i) = rc
+ 40 continue
+ 50 return
+ end
Added: branches/Interpolate1D/fitpack/fpader.f
===================================================================
--- branches/Interpolate1D/fitpack/fpader.f 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/fitpack/fpader.f 2008-07-18 19:44:12 UTC (rev 4550)
@@ -0,0 +1,57 @@
+ subroutine fpader(t,n,c,k1,x,l,d)
+c subroutine fpader calculates the derivatives
+c (j-1)
+c d(j) = s (x) , j=1,2,...,k1
+c of a spline of order k1 at the point t(l)<=x<t(l+1), using the
+c stable recurrence scheme of de boor
+c ..
+c ..scalar arguments..
+ real*8 x
+ integer n,k1,l
+c ..array arguments..
+ real*8 t(n),c(n),d(k1)
+c ..local scalars..
+ integer i,ik,j,jj,j1,j2,ki,kj,li,lj,lk
+ real*8 ak,fac,one
+c ..local array..
+ real*8 h(20)
+c ..
+ one = 0.1d+01
+ lk = l-k1
+ do 100 i=1,k1
+ ik = i+lk
+ h(i) = c(ik)
+ 100 continue
+ kj = k1
+ fac = one
+ do 700 j=1,k1
+ ki = kj
+ j1 = j+1
+ if(j.eq.1) go to 300
+ i = k1
+ do 200 jj=j,k1
+ li = i+lk
+ lj = li+kj
+ h(i) = (h(i)-h(i-1))/(t(lj)-t(li))
+ i = i-1
+ 200 continue
+ 300 do 400 i=j,k1
+ d(i) = h(i)
+ 400 continue
+ if(j.eq.k1) go to 600
+ do 500 jj=j1,k1
+ ki = ki-1
+ i = k1
+ do 500 j2=jj,k1
+ li = i+lk
+ lj = li+ki
+ d(i) = ((x-t(li))*d(i)+(t(lj)-x)*d(i-1))/(t(lj)-t(li))
+ i = i-1
+ 500 continue
+ 600 d(j) = d(k1)*fac
+ ak = k1-j
+ fac = fac*ak
+ kj = kj-1
+ 700 continue
+ return
+ end
Added: branches/Interpolate1D/fitpack/fpadno.f
===================================================================
--- branches/Interpolate1D/fitpack/fpadno.f 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/fitpack/fpadno.f 2008-07-18 19:44:12 UTC (rev 4550)
@@ -0,0 +1,59 @@
+ subroutine fpadno(maxtr,up,left,right,info,count,merk,jbind,
+ * n1,ier)
+c subroutine fpadno adds a branch of length n1 to the triply linked
+c tree,the information of which is kept in the arrays up,left,right
+c and info. the information field of the nodes of this new branch is
+c given in the array jbind. in linking the new branch fpadno takes
+c account of the property of the tree that
+c info(k) < info(right(k)) ; info(k) < info(left(k))
+c if necessary the subroutine calls subroutine fpfrno to collect the
+c free nodes of the tree. if no computer words are available at that
+c moment, the error parameter ier is set to 1.
+c ..
+c ..scalar arguments..
+ integer maxtr,count,merk,n1,ier
+c ..array arguments..
+ integer up(maxtr),left(maxtr),right(maxtr),info(maxtr),jbind(n1)
+c ..local scalars..
+ integer k,niveau,point
+ logical bool
+c ..subroutine references..
+c fpfrno
+c ..
+ point = 1
+ niveau = 1
+ 10 k = left(point)
+ bool = .true.
+ 20 if(k.eq.0) go to 50
+ if (info(k)-jbind(niveau).lt.0) go to 30
+ if (info(k)-jbind(niveau).eq.0) go to 40
+ go to 50
+ 30 point = k
+ k = right(point)
+ bool = .false.
+ go to 20
+ 40 point = k
+ niveau = niveau+1
+ go to 10
+ 50 if(niveau.gt.n1) go to 90
+ count = count+1
+ if(count.le.maxtr) go to 60
+ call fpfrno(maxtr,up,left,right,info,point,merk,n1,count,ier)
+ if(ier.ne.0) go to 100
+ 60 info(count) = jbind(niveau)
+ left(count) = 0
+ right(count) = k
+ if(bool) go to 70
+ bool = .true.
+ right(point) = count
+ up(count) = up(point)
+ go to 80
+ 70 up(count) = point
+ left(point) = count
+ 80 point = count
+ niveau = niveau+1
+ k = 0
+ go to 50
+ 90 ier = 0
+ 100 return
+ end
Added: branches/Interpolate1D/fitpack/fpadpo.f
===================================================================
--- branches/Interpolate1D/fitpack/fpadpo.f 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/fitpack/fpadpo.f 2008-07-18 19:44:12 UTC (rev 4550)
@@ -0,0 +1,70 @@
+ subroutine fpadpo(idim,t,n,c,nc,k,cp,np,cc,t1,t2)
+c given a idim-dimensional spline curve of degree k, in its b-spline
+c representation ( knots t(j),j=1,...,n , b-spline coefficients c(j),
+c j=1,...,nc) and given also a polynomial curve in its b-spline
+c representation ( coefficients cp(j), j=1,...,np), subroutine fpadpo
+c calculates the b-spline representation (coefficients c(j),j=1,...,nc)
+c of the sum of the two curves.
+c
+c other subroutine required : fpinst
+c
+c ..
+c ..scalar arguments..
+ integer idim,k,n,nc,np
+c ..array arguments..
+ real*8 t(n),c(nc),cp(np),cc(nc),t1(n),t2(n)
+c ..local scalars..
+ integer i,ii,j,jj,k1,l,l1,n1,n2,nk1,nk2
+c ..
+ k1 = k+1
+ nk1 = n-k1
+c initialization
+ j = 1
+ l = 1
+ do 20 jj=1,idim
+ l1 = j
+ do 10 ii=1,k1
+ cc(l1) = cp(l)
+ l1 = l1+1
+ l = l+1
+ 10 continue
+ j = j+n
+ l = l+k1
+ 20 continue
+ if(nk1.eq.k1) go to 70
+ n1 = k1*2
+ j = n
+ l = n1
+ do 30 i=1,k1
+ t1(i) = t(i)
+ t1(l) = t(j)
+ l = l-1
+ j = j-1
+ 30 continue
+c find the b-spline representation of the given polynomial curve
+c according to the given set of knots.
+ nk2 = nk1-1
+ do 60 l=k1,nk2
+ l1 = l+1
+ j = 1
+ do 40 i=1,idim
+ call fpinst(0,t1,n1,cc(j),k,t(l1),l,t2,n2,cc(j),n)
+ j = j+n
+ 40 continue
+ do 50 i=1,n2
+ t1(i) = t2(i)
+ 50 continue
+ n1 = n2
+ 60 continue
+c find the b-spline representation of the resulting curve.
+ 70 j = 1
+ do 90 jj=1,idim
+ l = j
+ do 80 i=1,nk1
+ c(l) = cc(l)+c(l)
+ l = l+1
+ 80 continue
+ j = j+n
+ 90 continue
+ return
+ end
Added: branches/Interpolate1D/fitpack/fpback.f
===================================================================
--- branches/Interpolate1D/fitpack/fpback.f 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/fitpack/fpback.f 2008-07-18 19:44:12 UTC (rev 4550)
@@ -0,0 +1,31 @@
+ subroutine fpback(a,z,n,k,c,nest)
+c subroutine fpback calculates the solution of the system of
+c equations a*c = z with a a n x n upper triangular matrix
+c of bandwidth k.
+c ..
+c ..scalar arguments..
+ integer n,k,nest
+c ..array arguments..
+ real*8 a(nest,k),z(n),c(n)
+c ..local scalars..
+ real*8 store
+ integer i,i1,j,k1,l,m
+c ..
+ k1 = k-1
+ c(n) = z(n)/a(n,1)
+ i = n-1
+ if(i.eq.0) go to 30
+ do 20 j=2,n
+ store = z(i)
+ i1 = k1
+ if(j.le.k1) i1 = j-1
+ m = i
+ do 10 l=1,i1
+ m = m+1
+ store = store-c(m)*a(i,l+1)
+ 10 continue
+ c(i) = store/a(i,1)
+ i = i-1
+ 20 continue
+ 30 return
+ end
Added: branches/Interpolate1D/fitpack/fpbacp.f
===================================================================
--- branches/Interpolate1D/fitpack/fpbacp.f 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/fitpack/fpbacp.f 2008-07-18 19:44:12 UTC (rev 4550)
@@ -0,0 +1,58 @@
+ subroutine fpbacp(a,b,z,n,k,c,k1,nest)
+c subroutine fpbacp calculates the solution of the system of equations
+c g * c = z with g a n x n upper triangular matrix of the form
+c ! a ' !
+c g = ! ' b !
+c ! 0 ' !
+c with b a n x k matrix and a a (n-k) x (n-k) upper triangular
+c matrix of bandwidth k1.
+c ..
+c ..scalar arguments..
+ integer n,k,k1,nest
+c ..array arguments..
+ real*8 a(nest,k1),b(nest,k),z(n),c(n)
+c ..local scalars..
+ integer i,i1,j,l,l0,l1,n2
+ real*8 store
+c ..
+ n2 = n-k
+ l = n
+ do 30 i=1,k
+ store = z(l)
+ j = k+2-i
+ if(i.eq.1) go to 20
+ l0 = l
+ do 10 l1=j,k
+ l0 = l0+1
+ store = store-c(l0)*b(l,l1)
+ 10 continue
+ 20 c(l) = store/b(l,j-1)
+ l = l-1
+ if(l.eq.0) go to 80
+ 30 continue
+ do 50 i=1,n2
+ store = z(i)
+ l = n2
+ do 40 j=1,k
+ l = l+1
+ store = store-c(l)*b(i,j)
+ 40 continue
+ c(i) = store
+ 50 continue
+ i = n2
+ c(i) = c(i)/a(i,1)
+ if(i.eq.1) go to 80
+ do 70 j=2,n2
+ i = i-1
+ store = c(i)
+ i1 = k
+ if(j.le.k) i1=j-1
+ l = i
+ do 60 l0=1,i1
+ l = l+1
+ store = store-c(l)*a(i,l0+1)
+ 60 continue
+ c(i) = store/a(i,1)
+ 70 continue
+ 80 return
+ end
Added: branches/Interpolate1D/fitpack/fpbfout.f
===================================================================
--- branches/Interpolate1D/fitpack/fpbfout.f 2008-07-18 13:19:37 UTC (rev 4549)
+++ branches/Interpolate1D/fitpack/fpbfout.f 2008-07-18 19:44:12 UTC (rev 4550)
@@ -0,0 +1,197 @@
+ subroutine fpbfou(t,n,par,ress,resc)
+c subroutine fpbfou calculates the integrals
+c /t(n-3)
+c ress(j) = ! nj,4(x)*sin(par*x) dx and
+c t(4)/
+c /t(n-3)
+c resc(j) = ! nj,4(x)*cos(par*x) dx , j=1,2,...n-4
+c t(4)/
+c where nj,4(x) denotes the cubic b-spline defined on the knots
+c t(j),t(j+1),...,t(j+4).
+c
+c calling sequence:
+c call fpbfou(t,n,par,ress,resc)
+c
+c input parameters:
+c t : real array,length n, containing the knots.
+c n : integer, containing the number of knots.
+c par : real, containing the value of the parameter par.
+c
+c output parameters:
+c ress : real array,length n, containing the integrals ress(j).
+c resc : real array,length n, containing the integrals resc(j).
+c
+c restrictions:
+c n >= 10, t(4) < t(5) < ... < t(n-4) < t(n-3).
+c ..
+c ..scalar arguments..
+ integer n
+ real*8 par
+c ..array arguments..
+ real*8 t(n),ress(n),resc(n)
+c ..local scalars..
+ integer i,ic,ipj,is,j,jj,jp1,jp4,k,li,lj,ll,nmj,nm3,nm7
+ real*8 ak,beta,con1,con2,c1,c2,delta,eps,fac,f1,f2,f3,one,quart,
+ * sign,six,s1,s2,term
+c ..local arrays..
+ real*8 co(5),si(5),hs(5),hc(5),rs(3),rc(3)
+c ..function references..
+ real*8 cos,sin,abs
+c ..
+c initialization.
+ one = 0.1e+01
+ six = 0.6e+01
+ eps = 0.1e-07
+ quart = 0.25e0
+ con1 = 0.5e-01
+ con2 = 0.12e+03
+ nm3 = n-3
+ nm7 = n-7
+ if(par.ne.0.) term = six/par
+ beta = par*t(4)
+ co(1) = cos(beta)
+ si(1) = sin(beta)
+c calculate the integrals ress(j) and resc(j), j=1,2,3 by setting up
+c a divided difference table.
+ do 30 j=1,3
+ jp1 = j+1
+ jp4 = j+4
+ beta = par*t(jp4)
+ co(jp1) = cos(beta)
+ si(jp1) = sin(beta)
+ call fpcsin(t(4),t(jp4),par,si(1),co(1),si(jp1),co(jp1),
+ * rs(j),rc(j))
+ i = 5-j
+ hs(i) = 0.
+ hc(i) = 0.
+ do 10 jj=1,j
+ ipj = i+jj
+ hs(ipj) = rs(jj)
+ hc(ipj) = rc(jj)
+ 10 continue
+ do 20 jj=1,3
+ if(i.lt.jj) i = jj
+ k = 5
+ li = jp4
+ do 20 ll=i,4
+ lj = li-jj
+ fac = t(li)-t(lj)
+ hs(k) = (hs(k)-hs(k-1))/fac
+ hc(k) = (hc(k)-hc(k-1))/fac
+ k = k-1
+ li = li-1
+ 20 continue
+ ress(j) = hs(5)-hs(4)
+ resc(j) = hc(5)-hc(4)
+ 30 continue
+ if(nm7.lt.4) go to 160
+c calculate the integrals ress(j) and resc(j),j=4,5,...,n-7.
+ do 150 j=4,nm7
+ jp4 = j+4
+ beta = par*t(jp4)
+ co(5) = cos(beta)
+ si(5) = sin(beta)
+ delta = t(jp4)-t(j)
+c the way of computing ress(j) and resc(j) depends on the value of
+c beta = par*(t(j+4)-t(j)).
+ beta = delta*par
+ if(abs(beta).le.one) go to 60
+c if !beta! > 1 the integrals are calculated by setting up a divided
+c difference table.
+ do 40 k=1,5
+ hs(k) = si(k)
+ hc(k) = co(k)
+ 40 continue
+ do 50 jj=1,3
+ k = 5
+ li = jp4
+ do 50 ll=jj,4
+ lj = li-jj
+ fac = par*(t(li)-t(lj))
+ hs(k) = (hs(k)-hs(k-1))/fac
+ hc(k) = (hc(k)-hc(k-1))/fac
+ k = k-1
+ li = li-1
+ 50 continue
+ s2 = (hs(5)-hs(4))*term
+ c2 = (hc(5)-hc(4))*term
+ go to 130
+c if !beta! <= 1 the integrals are calculated by evaluating a series
+c expansion.
+ 60 f3 = 0.
+ do 70 i=1,4
+ ipj = i+j
+ hs(i) = par*(t(ipj)-t(j))
+ hc(i) = hs(i)
+ f3 = f3+hs(i)
+ 70 continue
+ f3 = f3*con1
+ c1 = quart
+ s1 = f3
+ if(abs(f3).le.eps) go to 120
+ sign = one
+ fac = con2
+ k = 5
+ is = 0
+ do 110 ic=1,20
+ k = k+1
+ ak = k
+ fac = fac*ak
+ f1 = 0.
+ f3 = 0.
+ do 80 i=1,4
+ f1 = f1+hc(i)
+ f2 = f1*hs(i)
+ hc(i) = f2
+ f3 = f3+f2
+ 80 continue
+ f3 = f3*six/fac
+ if(is.eq.0) go to 90
+ is = 0
+ s1 = s1+f3*sign
+ go to 100
+ 90 sign = -sign
+ is = 1
+ c1 = c1+f3*sign
+ 100 if(abs(f3).le.eps) go to 120
+ 110 continue
+ 120 s2 = delta*(co(1)*s1+si(1)*c1)
+ c2 = delta*(co(1)*c1-si(1)*s1)
+ 130 ress(j) = s2
+ resc(j) = c2
+ do 140 i=1,4
+ co(i) = co(i+1)
+ si(i) = si(i+1)
+ 140 continue
+ 150 continue
+c calculate the integrals ress(j) and resc(j),j=n-6,n-5,n-4 by setting
+c up a divided difference table.
+ 160 do 190 j=1,3
+ nmj = nm3-j
+ i = 5-j
+ call fpcsin(t(nm3),t(nmj),par,si(4),co(4),si(i-1),co(i-1),
+ * rs(j),rc(j))
+ hs(i) = 0.
+ hc(i) = 0.