[Scipy-svn] r4610 - in branches/Interpolate1D: . extensions extensions/fitpack extensions/fitpack/fitpack
scipy-svn@scip...
scipy-svn@scip...
Thu Aug 7 16:25:34 CDT 2008
Author: fcady
Date: 2008-08-07 16:25:16 -0500 (Thu, 07 Aug 2008)
New Revision: 4610
Added:
branches/Interpolate1D/extensions/
branches/Interpolate1D/extensions/_fitpack.pyf
branches/Interpolate1D/extensions/_interpolate.cpp
branches/Interpolate1D/extensions/fitpack/
branches/Interpolate1D/extensions/fitpack/fitpack/
branches/Interpolate1D/extensions/interpolate.h
branches/Interpolate1D/extensions/multipack.h
branches/Interpolate1D/extensions/ndimage/
Removed:
branches/Interpolate1D/_fitpack.pyf
branches/Interpolate1D/_interpolate.cpp
branches/Interpolate1D/extensions/fitpack/fitpack/Makefile
branches/Interpolate1D/extensions/fitpack/fitpack/README
branches/Interpolate1D/extensions/fitpack/fitpack/bispev.f
branches/Interpolate1D/extensions/fitpack/fitpack/clocur.f
branches/Interpolate1D/extensions/fitpack/fitpack/cocosp.f
branches/Interpolate1D/extensions/fitpack/fitpack/concon.f
branches/Interpolate1D/extensions/fitpack/fitpack/concur.f
branches/Interpolate1D/extensions/fitpack/fitpack/cualde.f
branches/Interpolate1D/extensions/fitpack/fitpack/curev.f
branches/Interpolate1D/extensions/fitpack/fitpack/curfit.f
branches/Interpolate1D/extensions/fitpack/fitpack/dblint.f
branches/Interpolate1D/extensions/fitpack/fitpack/evapol.f
branches/Interpolate1D/extensions/fitpack/fitpack/fourco.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpader.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpadno.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpadpo.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpback.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpbacp.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpbfout.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpbisp.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpbspl.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpchec.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpched.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpchep.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpclos.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpcoco.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpcons.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpcosp.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpcsin.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpcurf.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpcuro.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpcyt1.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpcyt2.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpdeno.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpdisc.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpfrno.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpgivs.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpgrdi.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpgrpa.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpgrre.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpgrsp.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpinst.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpintb.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpknot.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpopdi.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpopsp.f
branches/Interpolate1D/extensions/fitpack/fitpack/fporde.f
branches/Interpolate1D/extensions/fitpack/fitpack/fppara.f
branches/Interpolate1D/extensions/fitpack/fitpack/fppasu.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpperi.f
branches/Interpolate1D/extensions/fitpack/fitpack/fppocu.f
branches/Interpolate1D/extensions/fitpack/fitpack/fppogr.f
branches/Interpolate1D/extensions/fitpack/fitpack/fppola.f
branches/Interpolate1D/extensions/fitpack/fitpack/fprank.f
branches/Interpolate1D/extensions/fitpack/fitpack/fprati.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpregr.f
branches/Interpolate1D/extensions/fitpack/fitpack/fprota.f
branches/Interpolate1D/extensions/fitpack/fitpack/fprppo.f
branches/Interpolate1D/extensions/fitpack/fitpack/fprpsp.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpseno.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpspgr.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpsphe.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpsuev.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpsurf.f
branches/Interpolate1D/extensions/fitpack/fitpack/fpsysy.f
branches/Interpolate1D/extensions/fitpack/fitpack/fptrnp.f
branches/Interpolate1D/extensions/fitpack/fitpack/fptrpe.f
branches/Interpolate1D/extensions/fitpack/fitpack/insert.f
branches/Interpolate1D/extensions/fitpack/fitpack/parcur.f
branches/Interpolate1D/extensions/fitpack/fitpack/parder.f
branches/Interpolate1D/extensions/fitpack/fitpack/parsur.f
branches/Interpolate1D/extensions/fitpack/fitpack/percur.f
branches/Interpolate1D/extensions/fitpack/fitpack/pogrid.f
branches/Interpolate1D/extensions/fitpack/fitpack/polar.f
branches/Interpolate1D/extensions/fitpack/fitpack/profil.f
branches/Interpolate1D/extensions/fitpack/fitpack/regrid.f
branches/Interpolate1D/extensions/fitpack/fitpack/spalde.f
branches/Interpolate1D/extensions/fitpack/fitpack/spgrid.f
branches/Interpolate1D/extensions/fitpack/fitpack/sphere.f
branches/Interpolate1D/extensions/fitpack/fitpack/splder.f
branches/Interpolate1D/extensions/fitpack/fitpack/splev.f
branches/Interpolate1D/extensions/fitpack/fitpack/splint.f
branches/Interpolate1D/extensions/fitpack/fitpack/sproot.f
branches/Interpolate1D/extensions/fitpack/fitpack/surev.f
branches/Interpolate1D/extensions/fitpack/fitpack/surfit.f
branches/Interpolate1D/fitpack/
branches/Interpolate1D/interpolate.h
branches/Interpolate1D/multipack.h
branches/Interpolate1D/ndimage/
Modified:
branches/Interpolate1D/setup.py
Log:
for cleanness of file, put all under-the-hood extensions into a single directory
Deleted: branches/Interpolate1D/_fitpack.pyf
===================================================================
--- branches/Interpolate1D/_fitpack.pyf 2008-08-07 21:05:37 UTC (rev 4609)
+++ branches/Interpolate1D/_fitpack.pyf 2008-08-07 21:25:16 UTC (rev 4610)
@@ -1,479 +0,0 @@
-! -*- f90 -*-
-! Author: Pearu Peterson <pearu@cens.ioc.ee>
-!
-python module _dfitpack ! in
-
- usercode '''
-
-static double dmax(double* seq,int len) {
- double val;
- int i;
- if (len<1)
- return -1e308;
- val = seq[0];
- for(i=1;i<len;++i)
- if (seq[i]>val) val = seq[i];
- return val;
-}
-static double dmin(double* seq,int len) {
- double val;
- int i;
- if (len<1)
- return 1e308;
- val = seq[0];
- for(i=1;i<len;++i)
- if (seq[i]<val) val = seq[i];
- return val;
-}
-static double calc_b(double* x,int m,double* tx,int nx) {
- double val1 = dmin(x,m);
- double val2 = dmin(tx,nx);
- if (val2>val1) return val1;
- val1 = dmax(tx,nx);
- return val2 - (val1-val2)/nx;
-}
-static double calc_e(double* x,int m,double* tx,int nx) {
- double val1 = dmax(x,m);
- double val2 = dmax(tx,nx);
- if (val2<val1) return val1;
- val1 = dmin(tx,nx);
- return val2 + (val2-val1)/nx;
-}
-static int imax(int i1,int i2) {
- return MAX(i1,i2);
-}
-
-static int calc_surfit_lwrk1(int m, int kx, int ky, int nxest, int nyest) {
- int u = nxest-kx-1;
- int v = nyest-ky-1;
- int km = MAX(kx,ky)+1;
- int ne = MAX(nxest,nyest);
- int bx = kx*v+ky+1;
- int by = ky*u+kx+1;
- int b1,b2;
- if (bx<=by) {b1=bx;b2=bx+v-ky;}
- else {b1=by;b2=by+u-kx;}
- return u*v*(2+b1+b2)+2*(u+v+km*(m+ne)+ne-kx-ky)+b2+1;
-}
-static int calc_surfit_lwrk2(int m, int kx, int ky, int nxest, int nyest) {
- int u = nxest-kx-1;
- int v = nyest-ky-1;
- int bx = kx*v+ky+1;
- int by = ky*u+kx+1;
- int b2 = (bx<=by?bx+v-ky:by+u-kx);
- return u*v*(b2+1)+b2;
-}
-
-static int calc_regrid_lwrk(int mx, int my, int kx, int ky,
- int nxest, int nyest) {
- int u = MAX(my, nxest);
- return 4+nxest*(my+2*kx+5)+nyest*(2*ky+5)+mx*(kx+1)+my*(ky+1)+u;
-}
-'''
-
- interface
-
- !!!!!!!!!! Univariate spline !!!!!!!!!!!
-
- subroutine splev(t,n,c,k,x,y,m,ier)
- ! y = splev(t,c,k,x)
- real*8 dimension(n),intent(in) :: t
- integer intent(hide),depend(t) :: n=len(t)
- real*8 dimension(n),depend(n,k),check(len(c)==n),intent(in) :: c
- integer :: k
- real*8 dimension(m),intent(in) :: x
- real*8 dimension(m),depend(m),intent(out) :: y
- integer intent(hide),depend(x) :: m=len(x)
- integer intent(hide) :: ier
- end subroutine splev
-
- subroutine splder(t,n,c,k,nu,x,y,m,wrk,ier)
- ! dy = splder(t,c,k,x,[nu])
- real*8 dimension(n) :: t
- integer depend(t),intent(hide) :: n=len(t)
- real*8 dimension(n),depend(n,k),check(len(c)==n),intent(in) :: c
- integer :: k
- integer depend(k),check(0<=nu && nu<=k) :: nu = 1
- real*8 dimension(m) :: x
- real*8 dimension(m),depend(m),intent(out) :: y
- integer depend(x),intent(hide) :: m=len(x)
- real*8 dimension(n),depend(n),intent(cache,hide) :: wrk
- integer intent(hide) :: ier
- end subroutine splder
-
- function splint(t,n,c,k,a,b,wrk)
- ! iy = splint(t,c,k,a,b)
- real*8 dimension(n),intent(in) :: t
- integer intent(hide),depend(t) :: n=len(t)
- real*8 dimension(n),depend(n),check(len(c)==n) :: c
- integer intent(in) :: k
- real*8 intent(in) :: a
- real*8 intent(in) :: b
- real*8 dimension(n),depend(n),intent(cache,hide) :: wrk
- real*8 :: splint
- end function splint
-
- subroutine sproot(t,n,c,zero,mest,m,ier)
- ! zero,m,ier = sproot(t,c[,mest])
- real*8 dimension(n),intent(in) :: t
- integer intent(hide),depend(t),check(n>=8) :: n=len(t)
- real*8 dimension(n),depend(n),check(len(c)==n) :: c
- real*8 dimension(mest),intent(out),depend(mest) :: zero
- integer optional,intent(in),depend(n) :: mest=3*(n-7)
- integer intent(out) :: m
- integer intent(out) :: ier
- end subroutine sproot
-
- subroutine spalde(t,n,c,k,x,d,ier)
- ! d,ier = spalde(t,c,k,x)
-
- callprotoargument double*,int*,double*,int*,double*,double*,int*
- callstatement {int k1=k+1; (*f2py_func)(t,&n,c,&k1,&x,d,&ier); }
-
- real*8 dimension(n) :: t
- integer intent(hide),depend(t) :: n=len(t)
- real*8 dimension(n),depend(n),check(len(c)==n) :: c
- integer intent(in) :: k
- real*8 intent(in) :: x
- real*8 dimension(k+1),intent(out),depend(k) :: d
- integer intent(out) :: ier
- end subroutine spalde
-
- subroutine curfit(iopt,m,x,y,w,xb,xe,k,s,nest,n,t,c,fp,wrk,lwrk,iwrk,ier)
- ! in curfit.f
- integer :: iopt
- integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
- real*8 dimension(m) :: x
- real*8 dimension(m),depend(m),check(len(y)==m) :: y
- real*8 dimension(m),depend(m),check(len(w)==m) :: w
- real*8 optional,depend(x),check(xb<=x[0]) :: xb = x[0]
- real*8 optional,depend(x,m),check(xe>=x[m-1]) :: xe = x[m-1]
- integer optional,check(1<=k && k <=5),intent(in) :: k=3
- real*8 optional,check(s>=0.0) :: s = 0.0
- integer intent(hide),depend(t) :: nest=len(t)
- integer intent(out), depend(nest) :: n=nest
- real*8 dimension(nest),intent(inout) :: t
- real*8 dimension(n),intent(out) :: c
- real*8 intent(out) :: fp
- real*8 dimension(lwrk),intent(inout) :: wrk
- integer intent(hide),depend(wrk) :: lwrk=len(wrk)
- integer dimension(nest),intent(inout) :: iwrk
- integer intent(out) :: ier
- end subroutine curfit
-
- subroutine percur(iopt,m,x,y,w,k,s,nest,n,t,c,fp,wrk,lwrk,iwrk,ier)
- ! in percur.f
- integer :: iopt
- integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
- real*8 dimension(m) :: x
- real*8 dimension(m),depend(m),check(len(y)==m) :: y
- real*8 dimension(m),depend(m),check(len(w)==m) :: w
- integer optional,check(1<=k && k <=5),intent(in) :: k=3
- real*8 optional,check(s>=0.0) :: s = 0.0
- integer intent(hide),depend(t) :: nest=len(t)
- integer intent(out), depend(nest) :: n=nest
- real*8 dimension(nest),intent(inout) :: t
- real*8 dimension(n),intent(out) :: c
- real*8 intent(out) :: fp
- real*8 dimension(lwrk),intent(inout) :: wrk
- integer intent(hide),depend(wrk) :: lwrk=len(wrk)
- integer dimension(nest),intent(inout) :: iwrk
- integer intent(out) :: ier
- end subroutine percur
-
-
- subroutine parcur(iopt,ipar,idim,m,u,mx,x,w,ub,ue,k,s,nest,n,t,nc,c,fp,wrk,lwrk,iwrk,ier)
- ! in parcur.f
- integer check(iopt>=-1 && iopt <= 1):: iopt
- integer check(ipar == 1 || ipar == 0) :: ipar
- integer check(idim > 0 && idim < 11) :: idim
- integer intent(hide),depend(u,k),check(m>k) :: m=len(u)
- real*8 dimension(m), intent(inout) :: u
- integer intent(hide),depend(x,idim,m),check(mx>=idim*m) :: mx=len(x)
- real*8 dimension(mx) :: x
- real*8 dimension(m) :: w
- real*8 :: ub
- real*8 :: ue
- integer optional, check(1<=k && k<=5) :: k=3.0
- real*8 optional, check(s>=0.0) :: s = 0.0
- integer intent(hide), depend(t) :: nest=len(t)
- integer intent(out), depend(nest) :: n=nest
- real*8 dimension(nest), intent(inout) :: t
- integer intent(hide), depend(c,nest,idim), check(nc>=idim*nest) :: nc=len(c)
- real*8 dimension(nc), intent(out) :: c
- real*8 intent(out) :: fp
- real*8 dimension(lwrk), intent(inout) :: wrk
- integer intent(hide),depend(wrk) :: lwrk=len(wrk)
- integer dimension(nest), intent(inout) :: iwrk
- integer intent(out) :: ier
- end subroutine parcur
-
-
- subroutine fpcurf0(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,c,fp,fpint,wrk,nrdata,ier)
- ! x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier = \
- ! fpcurf0(x,y,k,[w,xb,xe,s,nest])
-
- fortranname fpcurf
- callprotoargument int*,double*,double*,double*,int*,double*,double*,int*,double*,int*,double*,int*,int*,int*,int*,double*,double*,double*,double*,double*,double*,double*,double*,double*,int*,int*
- callstatement (*f2py_func)(&iopt,x,y,w,&m,&xb,&xe,&k,&s,&nest,&tol,&maxit,&k1,&k2,&n,t,c,&fp,fpint,wrk,wrk+nest,wrk+nest*k2,wrk+nest*2*k2,wrk+nest*3*k2,nrdata,&ier)
-
- integer intent(hide) :: iopt = 0
- real*8 dimension(m),intent(in,out) :: x
- real*8 dimension(m),depend(m),check(len(y)==m),intent(in,out) :: y
- real*8 dimension(m),depend(m),check(len(w)==m),intent(in,out) :: w = 1.0
- integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
- real*8 intent(in,out),depend(x),check(xb<=x[0]) :: xb = x[0]
- real*8 intent(in,out),depend(x,m),check(xe>=x[m-1]) :: xe = x[m-1]
- integer check(1<=k && k<=5),intent(in,out) :: k
- real*8 check(s>=0.0),depend(m),intent(in,out) :: s = m
- integer intent(in),depend(m,s,k,k1),check(nest>=2*k1) :: nest = (s==0.0?m+k+1:MAX(m/2,2*k1))
- real*8 intent(hide) :: tol = 0.001
- integer intent(hide) :: maxit = 20
- integer intent(hide),depend(k) :: k1=k+1
- integer intent(hide),depend(k) :: k2=k+2
- integer intent(out) :: n
- real*8 dimension(nest),intent(out),depend(nest) :: t
- real*8 dimension(nest),depend(nest),intent(out) :: c
- real*8 intent(out) :: fp
- real*8 dimension(nest),depend(nest),intent(out,cache) :: fpint
- real*8 dimension(nest*3*k2+m*k1),intent(cache,hide),depend(nest,k1,k2,m) :: wrk
- integer dimension(nest),depend(nest),intent(out,cache) :: nrdata
- integer intent(out) :: ier
- end subroutine fpcurf0
-
- subroutine fpcurf1(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,c,fp,fpint,wrk,nrdata,ier)
- ! x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier = \
- ! fpcurf1(x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier)
-
- fortranname fpcurf
- callprotoargument int*,double*,double*,double*,int*,double*,double*,int*,double*,int*,double*,int*,int*,int*,int*,double*,double*,double*,double*,double*,double*,double*,double*,double*,int*,int*
- callstatement (*f2py_func)(&iopt,x,y,w,&m,&xb,&xe,&k,&s,&nest,&tol,&maxit,&k1,&k2,&n,t,c,&fp,fpint,wrk,wrk+nest,wrk+nest*k2,wrk+nest*2*k2,wrk+nest*3*k2,nrdata,&ier)
-
- integer intent(hide) :: iopt = 1
- real*8 dimension(m),intent(in,out,overwrite) :: x
- real*8 dimension(m),depend(m),check(len(y)==m),intent(in,out,overwrite) :: y
- real*8 dimension(m),depend(m),check(len(w)==m),intent(in,out,overwrite) :: w
- integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
- real*8 intent(in,out) :: xb
- real*8 intent(in,out) :: xe
- integer check(1<=k && k<=5),intent(in,out) :: k
- real*8 check(s>=0.0),intent(in,out) :: s
- integer intent(hide),depend(t) :: nest = len(t)
- real*8 intent(hide) :: tol = 0.001
- integer intent(hide) :: maxit = 20
- integer intent(hide),depend(k) :: k1=k+1
- integer intent(hide),depend(k) :: k2=k+2
- integer intent(in,out) :: n
- real*8 dimension(nest),intent(in,out,overwrite) :: t
- real*8 dimension(nest),depend(nest),check(len(c)==nest),intent(in,out,overwrite) :: c
- real*8 intent(in,out) :: fp
- real*8 dimension(nest),depend(nest),check(len(fpint)==nest),intent(in,out,cache,overwrite) :: fpint
- real*8 dimension(nest*3*k2+m*k1),intent(cache,hide),depend(nest,k1,k2,m) :: wrk
- integer dimension(nest),depend(nest),check(len(nrdata)==nest),intent(in,out,cache,overwrite) :: nrdata
- integer intent(in,out) :: ier
- end subroutine fpcurf1
-
- subroutine fpcurfm1(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,c,fp,fpint,wrk,nrdata,ier)
- ! x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier = \
- ! fpcurfm1(x,y,k,t,[w,xb,xe])
-
- fortranname fpcurf
- callprotoargument int*,double*,double*,double*,int*,double*,double*,int*,double*,int*,double*,int*,int*,int*,int*,double*,double*,double*,double*,double*,double*,double*,double*,double*,int*,int*
- callstatement (*f2py_func)(&iopt,x,y,w,&m,&xb,&xe,&k,&s,&nest,&tol,&maxit,&k1,&k2,&n,t,c,&fp,fpint,wrk,wrk+nest,wrk+nest*k2,wrk+nest*2*k2,wrk+nest*3*k2,nrdata,&ier)
-
- integer intent(hide) :: iopt = -1
- real*8 dimension(m),intent(in,out) :: x
- real*8 dimension(m),depend(m),check(len(y)==m),intent(in,out) :: y
- real*8 dimension(m),depend(m),check(len(w)==m),intent(in,out) :: w = 1.0
- integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
- real*8 intent(in,out),depend(x),check(xb<=x[0]) :: xb = x[0]
- real*8 intent(in,out),depend(x,m),check(xe>=x[m-1]) :: xe = x[m-1]
- integer check(1<=k && k<=5),intent(in,out) :: k
- real*8 intent(out) :: s = -1
- integer intent(hide),depend(n) :: nest = n
- real*8 intent(hide) :: tol = 0.001
- integer intent(hide) :: maxit = 20
- integer intent(hide),depend(k) :: k1=k+1
- integer intent(hide),depend(k) :: k2=k+2
- integer intent(out),depend(t) :: n = len(t)
- real*8 dimension(n),intent(in,out,overwrite) :: t
- real*8 dimension(nest),depend(nest),intent(out) :: c
- real*8 intent(out) :: fp
- real*8 dimension(nest),depend(nest),intent(out,cache) :: fpint
- real*8 dimension(nest*3*k2+m*k1),intent(cache,hide),depend(nest,k1,k2,m) :: wrk
- integer dimension(nest),depend(nest),intent(out,cache) :: nrdata
- integer intent(out) :: ier
- end subroutine fpcurfm1
-
- !!!!!!!!!! Bivariate spline !!!!!!!!!!!
-
- subroutine bispev(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk,lwrk,iwrk,kwrk,ier)
- ! z,ier = bispev(tx,ty,c,kx,ky,x,y)
- real*8 dimension(nx),intent(in) :: tx
- integer intent(hide),depend(tx) :: nx=len(tx)
- real*8 dimension(ny),intent(in) :: ty
- integer intent(hide),depend(ty) :: ny=len(ty)
- real*8 intent(in),dimension((nx-kx-1)*(ny-ky-1)),depend(nx,ny,kx,ky),&
- check(len(c)==(nx-kx-1)*(ny-ky-1)):: c
- integer :: kx
- integer :: ky
- real*8 intent(in),dimension(mx) :: x
- integer intent(hide),depend(x) :: mx=len(x)
- real*8 intent(in),dimension(my) :: y
- integer intent(hide),depend(y) :: my=len(y)
- real*8 dimension(mx,my),depend(mx,my),intent(out,c) :: z
- real*8 dimension(lwrk),depend(lwrk),intent(hide,cache) :: wrk
- integer intent(hide),depend(mx,kx,my,ky) :: lwrk=mx*(kx+1)+my*(ky+1)
- integer dimension(kwrk),depend(kwrk),intent(hide,cache) :: iwrk
- integer intent(hide),depend(mx,my) :: kwrk=mx+my
- integer intent(out) :: ier
- end subroutine bispev
-
- subroutine surfit_smth(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,nyest,&
- nmax,eps,nx,tx,ny,ty,c,fp,wrk1,lwrk1,wrk2,lwrk2,&
- iwrk,kwrk,ier)
- ! nx,tx,ny,ty,c,fp,ier = surfit_smth(x,y,z,[w,xb,xe,yb,ye,kx,ky,s,eps,lwrk2])
-
- fortranname surfit
-
- integer intent(hide) :: iopt=0
- integer intent(hide),depend(x,kx,ky),check(m>=(kx+1)*(ky+1)) &
- :: m=len(x)
- real*8 dimension(m) :: x
- real*8 dimension(m),depend(m),check(len(y)==m) :: y
- real*8 dimension(m),depend(m),check(len(z)==m) :: z
- real*8 optional,dimension(m),depend(m),check(len(w)==m) :: w = 1.0
- real*8 optional,depend(x,m) :: xb=dmin(x,m)
- real*8 optional,depend(x,m) :: xe=dmax(x,m)
- real*8 optional,depend(y,m) :: yb=dmin(y,m)
- real*8 optional,depend(y,m) :: ye=dmax(y,m)
- integer check(1<=kx && kx<=5) :: kx = 3
- integer check(1<=ky && ky<=5) :: ky = 3
- real*8 optional,check(0.0<=s) :: s = m
- integer optional,depend(kx,m),check(nxest>=2*(kx+1)) &
- :: nxest = imax(kx+1+sqrt(m/2),2*(kx+1))
- integer optional,depend(ky,m),check(nyest>=2*(ky+1)) &
- :: nyest = imax(ky+1+sqrt(m/2),2*(ky+1))
- integer intent(hide),depend(nxest,nyest) :: nmax=MAX(nxest,nyest)
- real*8 optional,check(0.0<eps && eps<1.0) :: eps = 1e-16
- integer intent(out) :: nx
- real*8 dimension(nmax),intent(out),depend(nmax) :: tx
- integer intent(out) :: ny
- real*8 dimension(nmax),intent(out),depend(nmax) :: ty
- real*8 dimension((nxest-kx-1)*(nyest-ky-1)), &
- depend(kx,ky,nxest,nyest),intent(out) :: c
- real*8 intent(out) :: fp
- real*8 dimension(lwrk1),intent(cache,out),depend(lwrk1) :: wrk1
- integer intent(hide),depend(m,kx,ky,nxest,nyest) &
- :: lwrk1=calc_surfit_lwrk1(m,kx,ky,nxest,nyest)
- real*8 dimension(lwrk2),intent(cache,hide),depend(lwrk2) :: wrk2
- integer optional,intent(in),depend(kx,ky,nxest,nyest) &
- :: lwrk2=calc_surfit_lwrk2(m,kx,ky,nxest,nyest)
- integer dimension(kwrk),depend(kwrk),intent(cache,hide) :: iwrk
- integer intent(hide),depend(m,nxest,nyest,kx,ky) &
- :: kwrk=m+(nxest-2*kx-1)*(nyest-2*ky-1)
- integer intent(out) :: ier
- end subroutine surfit_smth
-
- subroutine surfit_lsq(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,nyest,&
- nmax,eps,nx,tx,ny,ty,c,fp,wrk1,lwrk1,wrk2,lwrk2,&
- iwrk,kwrk,ier)
- ! tx,ty,c,fp,ier = surfit_lsq(x,y,z,tx,ty,[w,xb,xe,yb,ye,kx,ky,eps,lwrk2])
-
- fortranname surfit
-
- integer intent(hide) :: iopt=-1
- integer intent(hide),depend(x,kx,ky),check(m>=(kx+1)*(ky+1)) &
- :: m=len(x)
- real*8 dimension(m) :: x
- real*8 dimension(m),depend(m),check(len(y)==m) :: y
- real*8 dimension(m),depend(m),check(len(z)==m) :: z
- real*8 optional,dimension(m),depend(m),check(len(w)==m) :: w = 1.0
- real*8 optional,depend(x,tx,m,nx) :: xb=calc_b(x,m,tx,nx)
- real*8 optional,depend(x,tx,m,nx) :: xe=calc_e(x,m,tx,nx)
- real*8 optional,depend(y,ty,m,ny) :: yb=calc_b(y,m,ty,ny)
- real*8 optional,depend(y,ty,m,ny) :: ye=calc_e(y,m,ty,ny)
- integer check(1<=kx && kx<=5) :: kx = 3
- integer check(1<=ky && ky<=5) :: ky = 3
- real*8 intent(hide) :: s = 0.0
- integer intent(hide),depend(nx) :: nxest = nx
- integer intent(hide),depend(ny) :: nyest = ny
- integer intent(hide),depend(nx,ny) :: nmax=MAX(nx,ny)
- real*8 optional,check(0.0<eps && eps<1.0) :: eps = 1e-16
- integer intent(hide),depend(tx,kx),check(2*kx+2<=nx) :: nx = len(tx)
- real*8 dimension(nx),intent(in,out,overwrite) :: tx
- integer intent(hide),depend(ty,ky),check(2*ky+2<=ny) :: ny = len(ty)
- real*8 dimension(ny),intent(in,out,overwrite) :: ty
- real*8 dimension((nx-kx-1)*(ny-ky-1)),depend(kx,ky,nx,ny),intent(out) :: c
- real*8 intent(out) :: fp
- real*8 dimension(lwrk1),intent(cache,hide),depend(lwrk1) :: wrk1
- integer intent(hide),depend(m,kx,ky,nxest,nyest) &
- :: lwrk1=calc_surfit_lwrk1(m,kx,ky,nxest,nyest)
- real*8 dimension(lwrk2),intent(cache,hide),depend(lwrk2) :: wrk2
- integer optional,intent(in),depend(kx,ky,nxest,nyest) &
- :: lwrk2=calc_surfit_lwrk2(m,kx,ky,nxest,nyest)
- integer dimension(kwrk),depend(kwrk),intent(cache,hide) :: iwrk
- integer intent(hide),depend(m,nx,ny,kx,ky) &
- :: kwrk=m+(nx-2*kx-1)*(ny-2*ky-1)
- integer intent(out) :: ier
- end subroutine surfit_lsq
-
- subroutine regrid_smth(iopt,mx,x,my,y,z,xb,xe,yb,ye,kx,ky,s,&
- nxest,nyest,nx,tx,ny,ty,c,fp,wrk,lwrk,iwrk,kwrk,ier)
- ! nx,tx,ny,ty,c,fp,ier = regrid_smth(x,y,z,[xb,xe,yb,ye,kx,ky,s])
-
- fortranname regrid
-
- integer intent(hide) :: iopt=0
- integer intent(hide),depend(x,kx),check(mx>kx) :: mx=len(x)
- real*8 dimension(mx) :: x
- integer intent(hide),depend(y,ky),check(my>ky) :: my=len(y)
- real*8 dimension(my) :: y
- real*8 dimension(mx*my),depend(mx,my),check(len(z)==mx*my) :: z
- real*8 optional,depend(x,mx) :: xb=dmin(x,mx)
- real*8 optional,depend(x,mx) :: xe=dmax(x,mx)
- real*8 optional,depend(y,my) :: yb=dmin(y,my)
- real*8 optional,depend(y,my) :: ye=dmax(y,my)
- integer optional,check(1<=kx && kx<=5) :: kx = 3
- integer optional,check(1<=ky && ky<=5) :: ky = 3
- real*8 optional,check(0.0<=s) :: s = 0.0
- integer intent(hide),depend(kx,mx),check(nxest>=2*(kx+1)) &
- :: nxest = mx+kx+1
- integer intent(hide),depend(ky,my),check(nyest>=2*(ky+1)) &
- :: nyest = my+ky+1
- integer intent(out) :: nx
- real*8 dimension(nxest),intent(out),depend(nxest) :: tx
- integer intent(out) :: ny
- real*8 dimension(nyest),intent(out),depend(nyest) :: ty
- real*8 dimension((nxest-kx-1)*(nyest-ky-1)), &
- depend(kx,ky,nxest,nyest),intent(out) :: c
- real*8 intent(out) :: fp
- real*8 dimension(lwrk),intent(cache,hide),depend(lwrk) :: wrk
- integer intent(hide),depend(mx,my,kx,ky,nxest,nyest) &
- :: lwrk=calc_regrid_lwrk(mx,my,kx,ky,nxest,nyest)
- integer dimension(kwrk),depend(kwrk),intent(cache,hide) :: iwrk
- integer intent(hide),depend(mx,my,nxest,nyest) &
- :: kwrk=3+mx+my+nxest+nyest
- integer intent(out) :: ier
- end subroutine regrid_smth
-
- function dblint(tx,nx,ty,ny,c,kx,ky,xb,xe,yb,ye,wrk)
- ! iy = dblint(tx,ty,c,kx,ky,xb,xe,yb,ye)
- real*8 dimension(nx),intent(in) :: tx
- integer intent(hide),depend(tx) :: nx=len(tx)
- real*8 dimension(ny),intent(in) :: ty
- integer intent(hide),depend(ty) :: ny=len(ty)
- real*8 intent(in),dimension((nx-kx-1)*(ny-ky-1)),depend(nx,ny,kx,ky),&
- check(len(c)==(nx-kx-1)*(ny-ky-1)):: c
- integer :: kx
- integer :: ky
- real*8 intent(in) :: xb
- real*8 intent(in) :: xe
- real*8 intent(in) :: yb
- real*8 intent(in) :: ye
- real*8 dimension(nx+ny-kx-ky-2),depend(nx,ny,kx,ky),intent(cache,hide) :: wrk
- real*8 :: dblint
- end function dblint
- end interface
-end python module _dfitpack
-
Deleted: branches/Interpolate1D/_interpolate.cpp
===================================================================
--- branches/Interpolate1D/_interpolate.cpp 2008-08-07 21:05:37 UTC (rev 4609)
+++ branches/Interpolate1D/_interpolate.cpp 2008-08-07 21:25:16 UTC (rev 4610)
@@ -1,236 +0,0 @@
-#include "Python.h"
-#include <stdlib.h>
-
-#include "interpolate.h"
-#include "numpy/arrayobject.h"
-
-using namespace std;
-
-extern "C" {
-
-static PyObject* linear_method(PyObject*self, PyObject* args, PyObject* kywds)
-{
- static char *kwlist[] = {"x","y","new_x","new_y", NULL};
- PyObject *py_x, *py_y, *py_new_x, *py_new_y;
- py_x = py_y = py_new_x = py_new_y = NULL;
- PyObject *arr_x, *arr_y, *arr_new_x, *arr_new_y;
- arr_x = arr_y = arr_new_x = arr_new_y = NULL;
-
- if(!PyArg_ParseTupleAndKeywords(args,kywds,"OOOO:linear_dddd",kwlist,&py_x, &py_y, &py_new_x, &py_new_y))
- return NULL;
- arr_x = PyArray_FROMANY(py_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
- if (!arr_x) {
- PyErr_SetString(PyExc_ValueError, "x must be a 1-D array of floats");
- goto fail;
- }
- arr_y = PyArray_FROMANY(py_y, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
- if (!arr_y) {
- PyErr_SetString(PyExc_ValueError, "y must be a 1-D array of floats");
- goto fail;
- }
- arr_new_x = PyArray_FROMANY(py_new_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
- if (!arr_new_x) {
- PyErr_SetString(PyExc_ValueError, "new_x must be a 1-D array of floats");
- goto fail;
- }
- arr_new_y = PyArray_FROMANY(py_new_y, PyArray_DOUBLE, 1, 1, NPY_INOUT_ARRAY);
- if (!arr_new_y) {
- PyErr_SetString(PyExc_ValueError, "new_y must be a 1-D array of floats");
- goto fail;
- }
-
- linear((double*)PyArray_DATA(arr_x), (double*)PyArray_DATA(arr_y),
- PyArray_DIM(arr_x,0), (double*)PyArray_DATA(arr_new_x),
- (double*)PyArray_DATA(arr_new_y), PyArray_DIM(arr_new_x,0));
-
- Py_DECREF(arr_x);
- Py_DECREF(arr_y);
- Py_DECREF(arr_new_x);
- Py_DECREF(arr_new_y);
-
- Py_RETURN_NONE;
-
-fail:
- Py_XDECREF(arr_x);
- Py_XDECREF(arr_y);
- Py_XDECREF(arr_new_x);
- Py_XDECREF(arr_new_y);
- return NULL;
-}
-
-static PyObject* loginterp_method(PyObject*self, PyObject* args, PyObject* kywds)
-{
- static char *kwlist[] = {"x","y","new_x","new_y", NULL};
- PyObject *py_x, *py_y, *py_new_x, *py_new_y;
- py_x = py_y = py_new_x = py_new_y = NULL;
- PyObject *arr_x, *arr_y, *arr_new_x, *arr_new_y;
- arr_x = arr_y = arr_new_x = arr_new_y = NULL;
-
- if(!PyArg_ParseTupleAndKeywords(args,kywds,"OOOO:loginterp_dddd",kwlist,&py_x, &py_y, &py_new_x, &py_new_y))
- return NULL;
- arr_x = PyArray_FROMANY(py_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
- if (!arr_x) {
- PyErr_SetString(PyExc_ValueError, "x must be a 1-D array of floats");
- goto fail;
- }
- arr_y = PyArray_FROMANY(py_y, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
- if (!arr_y) {
- PyErr_SetString(PyExc_ValueError, "y must be a 1-D array of floats");
- goto fail;
- }
- arr_new_x = PyArray_FROMANY(py_new_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
- if (!arr_new_x) {
- PyErr_SetString(PyExc_ValueError, "new_x must be a 1-D array of floats");
- goto fail;
- }
- arr_new_y = PyArray_FROMANY(py_new_y, PyArray_DOUBLE, 1, 1, NPY_INOUT_ARRAY);
- if (!arr_new_y) {
- PyErr_SetString(PyExc_ValueError, "new_y must be a 1-D array of floats");
- goto fail;
- }
-
- loginterp((double*)PyArray_DATA(arr_x), (double*)PyArray_DATA(arr_y),
- PyArray_DIM(arr_x,0), (double*)PyArray_DATA(arr_new_x),
- (double*)PyArray_DATA(arr_new_y), PyArray_DIM(arr_new_x,0));
-
- Py_DECREF(arr_x);
- Py_DECREF(arr_y);
- Py_DECREF(arr_new_x);
- Py_DECREF(arr_new_y);
-
- Py_RETURN_NONE;
-
-fail:
- Py_XDECREF(arr_x);
- Py_XDECREF(arr_y);
- Py_XDECREF(arr_new_x);
- Py_XDECREF(arr_new_y);
- return NULL;
-}
-
-static PyObject* window_average_method(PyObject*self, PyObject* args, PyObject* kywds)
-{
- static char *kwlist[] = {"x","y","new_x","new_y", NULL};
- PyObject *py_x, *py_y, *py_new_x, *py_new_y;
- py_x = py_y = py_new_x = py_new_y = NULL;
- PyObject *arr_x, *arr_y, *arr_new_x, *arr_new_y;
- arr_x = arr_y = arr_new_x = arr_new_y = NULL;
- double width;
-
- if(!PyArg_ParseTupleAndKeywords(args,kywds,"OOOOd:loginterp_dddd",kwlist,&py_x, &py_y, &py_new_x, &py_new_y, &width))
- return NULL;
- arr_x = PyArray_FROMANY(py_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
- if (!arr_x) {
- PyErr_SetString(PyExc_ValueError, "x must be a 1-D array of floats");
- goto fail;
- }
- arr_y = PyArray_FROMANY(py_y, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
- if (!arr_y) {
- PyErr_SetString(PyExc_ValueError, "y must be a 1-D array of floats");
- goto fail;
- }
- arr_new_x = PyArray_FROMANY(py_new_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
- if (!arr_new_x) {
- PyErr_SetString(PyExc_ValueError, "new_x must be a 1-D array of floats");
- goto fail;
- }
- arr_new_y = PyArray_FROMANY(py_new_y, PyArray_DOUBLE, 1, 1, NPY_INOUT_ARRAY);
- if (!arr_new_y) {
- PyErr_SetString(PyExc_ValueError, "new_y must be a 1-D array of floats");
- goto fail;
- }
-
- window_average((double*)PyArray_DATA(arr_x), (double*)PyArray_DATA(arr_y),
- PyArray_DIM(arr_x,0), (double*)PyArray_DATA(arr_new_x),
- (double*)PyArray_DATA(arr_new_y), PyArray_DIM(arr_new_x,0), width);
-
- Py_DECREF(arr_x);
- Py_DECREF(arr_y);
- Py_DECREF(arr_new_x);
- Py_DECREF(arr_new_y);
-
- Py_RETURN_NONE;
-
-fail:
- Py_XDECREF(arr_x);
- Py_XDECREF(arr_y);
- Py_XDECREF(arr_new_x);
- Py_XDECREF(arr_new_y);
- return NULL;
-}
-
-static PyObject* block_average_above_method(PyObject*self, PyObject* args, PyObject* kywds)
-{
- static char *kwlist[] = {"x","y","new_x","new_y", NULL};
- PyObject *py_x, *py_y, *py_new_x, *py_new_y;
- py_x = py_y = py_new_x = py_new_y = NULL;
- PyObject *arr_x, *arr_y, *arr_new_x, *arr_new_y;
- arr_x = arr_y = arr_new_x = arr_new_y = NULL;
-
- if(!PyArg_ParseTupleAndKeywords(args,kywds,"OOOO:loginterp_dddd",kwlist,&py_x, &py_y, &py_new_x, &py_new_y))
- return NULL;
- arr_x = PyArray_FROMANY(py_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
- if (!arr_x) {
- PyErr_SetString(PyExc_ValueError, "x must be a 1-D array of floats");
- goto fail;
- }
- arr_y = PyArray_FROMANY(py_y, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
- if (!arr_y) {
- PyErr_SetString(PyExc_ValueError, "y must be a 1-D array of floats");
- goto fail;
- }
- arr_new_x = PyArray_FROMANY(py_new_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
- if (!arr_new_x) {
- PyErr_SetString(PyExc_ValueError, "new_x must be a 1-D array of floats");
- goto fail;
- }
- arr_new_y = PyArray_FROMANY(py_new_y, PyArray_DOUBLE, 1, 1, NPY_INOUT_ARRAY);
- if (!arr_new_y) {
- PyErr_SetString(PyExc_ValueError, "new_y must be a 1-D array of floats");
- goto fail;
- }
-
- block_average_above((double*)PyArray_DATA(arr_x), (double*)PyArray_DATA(arr_y),
- PyArray_DIM(arr_x,0), (double*)PyArray_DATA(arr_new_x),
- (double*)PyArray_DATA(arr_new_y), PyArray_DIM(arr_new_x,0));
-
- Py_DECREF(arr_x);
- Py_DECREF(arr_y);
- Py_DECREF(arr_new_x);
- Py_DECREF(arr_new_y);
-
- Py_RETURN_NONE;
-
-fail:
- Py_XDECREF(arr_x);
- Py_XDECREF(arr_y);
- Py_XDECREF(arr_new_x);
- Py_XDECREF(arr_new_y);
- return NULL;
-}
-
-static PyMethodDef interpolate_methods[] = {
- {"linear_dddd", (PyCFunction)linear_method, METH_VARARGS|METH_KEYWORDS,
- ""},
- {"loginterp_dddd", (PyCFunction)loginterp_method, METH_VARARGS|METH_KEYWORDS,
- ""},
- {"window_average_ddddd", (PyCFunction)window_average_method, METH_VARARGS|METH_KEYWORDS,
- ""},
- {"block_average_above_dddd", (PyCFunction)block_average_above_method, METH_VARARGS|METH_KEYWORDS,
- ""},
- {NULL, NULL, 0, NULL}
-};
-
-
-PyMODINIT_FUNC init_interpolate(void)
-{
- PyObject* m;
- m = Py_InitModule3("_interpolate", interpolate_methods,
- "A few interpolation routines.\n"
- );
- if (m == NULL)
- return;
- import_array();
-}
-
-} // extern "C"
Copied: branches/Interpolate1D/extensions/_fitpack.pyf (from rev 4591, branches/Interpolate1D/_fitpack.pyf)
Copied: branches/Interpolate1D/extensions/_interpolate.cpp (from rev 4587, branches/Interpolate1D/_interpolate.cpp)
Copied: branches/Interpolate1D/extensions/fitpack (from rev 4587, branches/Interpolate1D/fitpack)
Copied: branches/Interpolate1D/extensions/fitpack/fitpack (from rev 4587, branches/Interpolate1D/fitpack)
Deleted: branches/Interpolate1D/extensions/fitpack/fitpack/Makefile
===================================================================
--- branches/Interpolate1D/fitpack/Makefile 2008-07-31 19:09:00 UTC (rev 4587)
+++ branches/Interpolate1D/extensions/fitpack/fitpack/Makefile 2008-08-07 21:25:16 UTC (rev 4610)
@@ -1,19 +0,0 @@
-# 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
-
-
-
-
-
-
Deleted: branches/Interpolate1D/extensions/fitpack/fitpack/README
===================================================================
--- branches/Interpolate1D/fitpack/README 2008-07-31 19:09:00 UTC (rev 4587)
+++ branches/Interpolate1D/extensions/fitpack/fitpack/README 2008-08-07 21:25:16 UTC (rev 4610)
@@ -1,3 +0,0 @@
-- ddierckx is a 'real*8' version of dierckx
- generated by Pearu Peterson <pearu@ioc.ee>.
-- dierckx (in netlib) is fitpack by P. Dierckx
Deleted: branches/Interpolate1D/extensions/fitpack/fitpack/bispev.f
===================================================================
--- branches/Interpolate1D/fitpack/bispev.f 2008-07-31 19:09:00 UTC (rev 4587)
+++ branches/Interpolate1D/extensions/fitpack/fitpack/bispev.f 2008-08-07 21:25:16 UTC (rev 4610)
@@ -1,103 +0,0 @@
- 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
Deleted: branches/Interpolate1D/extensions/fitpack/fitpack/clocur.f
===================================================================
--- branches/Interpolate1D/fitpack/clocur.f 2008-07-31 19:09:00 UTC (rev 4587)
+++ branches/Interpolate1D/extensions/fitpack/fitpack/clocur.f 2008-08-07 21:25:16 UTC (rev 4610)
@@ -1,352 +0,0 @@
- 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
Deleted: branches/Interpolate1D/extensions/fitpack/fitpack/cocosp.f
===================================================================
--- branches/Interpolate1D/fitpack/cocosp.f 2008-07-31 19:09:00 UTC (rev 4587)
+++ branches/Interpolate1D/extensions/fitpack/fitpack/cocosp.f 2008-08-07 21:25:16 UTC (rev 4610)
@@ -1,180 +0,0 @@
- 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
Deleted: branches/Interpolate1D/extensions/fitpack/fitpack/concon.f
===================================================================
--- branches/Interpolate1D/fitpack/concon.f 2008-07-31 19:09:00 UTC (rev 4587)
+++ branches/Interpolate1D/extensions/fitpack/fitpack/concon.f 2008-08-07 21:25:16 UTC (rev 4610)
@@ -1,233 +0,0 @@
- 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
Deleted: branches/Interpolate1D/extensions/fitpack/fitpack/concur.f
===================================================================
--- branches/Interpolate1D/fitpack/concur.f 2008-07-31 19:09:00 UTC (rev 4587)
+++ branches/Interpolate1D/extensions/fitpack/fitpack/concur.f 2008-08-07 21:25:16 UTC (rev 4610)
@@ -1,370 +0,0 @@
- 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
Deleted: branches/Interpolate1D/extensions/fitpack/fitpack/cualde.f
===================================================================
--- branches/Interpolate1D/fitpack/cualde.f 2008-07-31 19:09:00 UTC (rev 4587)
+++ branches/Interpolate1D/extensions/fitpack/fitpack/cualde.f 2008-08-07 21:25:16 UTC (rev 4610)
@@ -1,91 +0,0 @@
- 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
Deleted: branches/Interpolate1D/extensions/fitpack/fitpack/curev.f
===================================================================
--- branches/Interpolate1D/fitpack/curev.f 2008-07-31 19:09:00 UTC (rev 4587)
+++ branches/Interpolate1D/extensions/fitpack/fitpack/curev.f 2008-08-07 21:25:16 UTC (rev 4610)
@@ -1,110 +0,0 @@
- 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
Deleted: branches/Interpolate1D/extensions/fitpack/fitpack/curfit.f
===================================================================
--- branches/Interpolate1D/fitpack/curfit.f 2008-07-31 19:09:00 UTC (rev 4587)
+++ branches/Interpolate1D/extensions/fitpack/fitpack/curfit.f 2008-08-07 21:25:16 UTC (rev 4610)
@@ -1,261 +0,0 @@
- 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
Deleted: branches/Interpolate1D/extensions/fitpack/fitpack/dblint.f
===================================================================
--- branches/Interpolate1D/fitpack/dblint.f 2008-07-31 19:09:00 UTC (rev 4587)
+++ branches/Interpolate1D/extensions/fitpack/fitpack/dblint.f 2008-08-07 21:25:16 UTC (rev 4610)
@@ -1,88 +0,0 @@
- 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 integra