b1accf0ae3
This adds a new `splines` module which implements a high-level interface for creating and evaluating splines and rewrite almost all modules to use it. Also, notably: 1. both `simplespline` and DIERCKX splines can now used with a uniform interface 2. most complexity due to handling working space arrays is gone 3. memory management has been significantly simplified too
4610 lines
157 KiB
Fortran
4610 lines
157 KiB
Fortran
module dierckx
|
|
|
|
use const_and_precisions, only : wp_
|
|
implicit none
|
|
|
|
contains
|
|
|
|
subroutine bispev(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk,lwrk, &
|
|
iwrk,kwrk,ier)
|
|
! subroutine bispev evaluates on a grid (x(i),y(j)),i=1,...,mx; j=1,...
|
|
! ,my a bivariate spline s(x,y) of degrees kx and ky, given in the
|
|
! b-spline representation.
|
|
!
|
|
! calling sequence:
|
|
! call bispev(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk,lwrk,
|
|
! * iwrk,kwrk,ier)
|
|
!
|
|
! input parameters:
|
|
! tx : real array, length nx, which contains the position of the
|
|
! knots in the x-direction.
|
|
! nx : integer, giving the total number of knots in the x-direction
|
|
! ty : real array, length ny, which contains the position of the
|
|
! knots in the y-direction.
|
|
! ny : integer, giving the total number of knots in the y-direction
|
|
! c : real array, length (nx-kx-1)*(ny-ky-1), which contains the
|
|
! b-spline coefficients.
|
|
! kx,ky : integer values, giving the degrees of the spline.
|
|
! x : real array of dimension (mx).
|
|
! before entry x(i) must be set to the x co-ordinate of the
|
|
! i-th grid point along the x-axis.
|
|
! tx(kx+1)<=x(i-1)<=x(i)<=tx(nx-kx), i=2,...,mx.
|
|
! mx : on entry mx must specify the number of grid points along
|
|
! the x-axis. mx >=1.
|
|
! y : real array of dimension (my).
|
|
! before entry y(j) must be set to the y co-ordinate of the
|
|
! j-th grid point along the y-axis.
|
|
! ty(ky+1)<=y(j-1)<=y(j)<=ty(ny-ky), j=2,...,my.
|
|
! my : on entry my must specify the number of grid points along
|
|
! the y-axis. my >=1.
|
|
! wrk : real array of dimension lwrk. used as workspace.
|
|
! lwrk : integer, specifying the dimension of wrk.
|
|
! lwrk >= mx*(kx+1)+my*(ky+1)
|
|
! iwrk : integer array of dimension kwrk. used as workspace.
|
|
! kwrk : integer, specifying the dimension of iwrk. kwrk >= mx+my.
|
|
!
|
|
! output parameters:
|
|
! z : real array of dimension (mx*my).
|
|
! on succesful exit z(my*(i-1)+j) contains the value of s(x,y)
|
|
! at the point (x(i),y(j)),i=1,...,mx;j=1,...,my.
|
|
! ier : integer error flag
|
|
! ier=0 : normal return
|
|
! ier=10: invalid input data (see restrictions)
|
|
!
|
|
! restrictions:
|
|
! mx >=1, my >=1, lwrk>=mx*(kx+1)+my*(ky+1), kwrk>=mx+my
|
|
! tx(kx+1) <= x(i-1) <= x(i) <= tx(nx-kx), i=2,...,mx
|
|
! ty(ky+1) <= y(j-1) <= y(j) <= ty(ny-ky), j=2,...,my
|
|
!
|
|
! other subroutines required:
|
|
! fpbisp,fpbspl
|
|
!
|
|
! references :
|
|
! de boor c : on calculating with b-splines, j. approximation theory
|
|
! 6 (1972) 50-62.
|
|
! cox m.g. : the numerical evaluation of b-splines, j. inst. maths
|
|
! applics 10 (1972) 134-149.
|
|
! dierckx p. : curve and surface fitting with splines, monographs on
|
|
! numerical analysis, oxford university press, 1993.
|
|
!
|
|
! author :
|
|
! p.dierckx
|
|
! dept. computer science, k.u.leuven
|
|
! celestijnenlaan 200a, b-3001 heverlee, belgium.
|
|
! e-mail : Paul.Dierckx@cs.kuleuven.ac.be
|
|
!
|
|
! latest update : march 1987
|
|
!
|
|
implicit none
|
|
! arguments
|
|
integer, intent(in) :: nx, ny, kx, ky, mx, my, lwrk, kwrk
|
|
integer, intent(out) :: ier
|
|
integer, intent(inout) :: iwrk(kwrk)
|
|
real(wp_), intent(in) :: tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1)),x(mx), y(my)
|
|
real(wp_), intent(out) :: z(mx*my)
|
|
real(wp_), intent(inout) :: wrk(lwrk)
|
|
! local variables
|
|
integer :: i, iw, lwest
|
|
! ..
|
|
! before starting computations a data check is made. if the input data
|
|
! are invalid control is immediately repassed to the calling program.
|
|
ier = 10
|
|
lwest = (kx+1)*mx+(ky+1)*my
|
|
if(lwrk<lwest .or. kwrk<(mx+my) .or. mx<1 .or. my<1) return
|
|
do i=2,mx
|
|
if(x(i)<x(i-1)) return
|
|
end do
|
|
do i=2,my
|
|
if(y(i)<y(i-1)) return
|
|
end do
|
|
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))
|
|
end subroutine bispev
|
|
|
|
subroutine surfit(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)
|
|
! given the set of data points (x(i),y(i),z(i)) and the set of positive
|
|
! numbers w(i),i=1,...,m, subroutine surfit determines a smooth bivar-
|
|
! iate spline approximation s(x,y) of degrees kx and ky on the rect-
|
|
! angle xb <= x <= xe, yb <= y <= ye.
|
|
! if iopt = -1 surfit calculates the weighted least-squares spline
|
|
! according to a given set of knots.
|
|
! if iopt >= 0 the total numbers nx and ny of these knots and their
|
|
! position tx(j),j=1,...,nx and ty(j),j=1,...,ny are chosen automatic-
|
|
! ally by the routine. the smoothness of s(x,y) is then achieved by
|
|
! minimalizing the discontinuity jumps in the derivatives of s(x,y)
|
|
! across the boundaries of the subpanels (tx(i),tx(i+1))*(ty(j),ty(j+1).
|
|
! the amounth of smoothness is determined by the condition that f(p) =
|
|
! sum ((w(i)*(z(i)-s(x(i),y(i))))**2) be <= s, with s a given non-neg-
|
|
! ative constant, called the smoothing factor.
|
|
! the fit is given in the b-spline representation (b-spline coefficients
|
|
! c((ny-ky-1)*(i-1)+j),i=1,...,nx-kx-1;j=1,...,ny-ky-1) and can be eval-
|
|
! uated by means of subroutine bispev.
|
|
!
|
|
! calling sequence:
|
|
! call surfit(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)
|
|
!
|
|
! parameters:
|
|
! iopt : integer flag. on entry iopt must specify whether a weighted
|
|
! least-squares spline (iopt=-1) or a smoothing spline (iopt=0
|
|
! or 1) must be determined.
|
|
! if iopt=0 the routine will start with an initial set of knots
|
|
! tx(i)=xb,tx(i+kx+1)=xe,i=1,...,kx+1;ty(i)=yb,ty(i+ky+1)=ye,i=
|
|
! 1,...,ky+1. if iopt=1 the routine will continue with the set
|
|
! of knots found at the last call of the routine.
|
|
! attention: a call with iopt=1 must always be immediately pre-
|
|
! ceded by another call with iopt=1 or iopt=0.
|
|
! unchanged on exit.
|
|
! m : integer. on entry m must specify the number of data points.
|
|
! m >= (kx+1)*(ky+1). unchanged on exit.
|
|
! x : real array of dimension at least (m).
|
|
! y : real array of dimension at least (m).
|
|
! z : real array of dimension at least (m).
|
|
! before entry, x(i),y(i),z(i) must be set to the co-ordinates
|
|
! of the i-th data point, for i=1,...,m. the order of the data
|
|
! points is immaterial. unchanged on exit.
|
|
! w : real array of dimension at least (m). before entry, w(i) must
|
|
! be set to the i-th value in the set of weights. the w(i) must
|
|
! be strictly positive. unchanged on exit.
|
|
! xb,xe : real values. on entry xb,xe,yb and ye must specify the bound-
|
|
! yb,ye aries of the rectangular approximation domain.
|
|
! xb<=x(i)<=xe,yb<=y(i)<=ye,i=1,...,m. unchanged on exit.
|
|
! kx,ky : integer values. on entry kx and ky must specify the degrees
|
|
! of the spline. 1<=kx,ky<=5. it is recommended to use bicubic
|
|
! (kx=ky=3) splines. unchanged on exit.
|
|
! s : real. on entry (in case iopt>=0) s must specify the smoothing
|
|
! factor. s >=0. unchanged on exit.
|
|
! for advice on the choice of s see further comments
|
|
! nxest : integer. unchanged on exit.
|
|
! nyest : integer. unchanged on exit.
|
|
! on entry, nxest and nyest must specify an upper bound for the
|
|
! number of knots required in the x- and y-directions respect.
|
|
! these numbers will also determine the storage space needed by
|
|
! the routine. nxest >= 2*(kx+1), nyest >= 2*(ky+1).
|
|
! in most practical situation nxest = kx+1+sqrt(m/2), nyest =
|
|
! ky+1+sqrt(m/2) will be sufficient. see also further comments.
|
|
! nmax : integer. on entry nmax must specify the actual dimension of
|
|
! the arrays tx and ty. nmax >= nxest, nmax >=nyest.
|
|
! unchanged on exit.
|
|
! eps : real.
|
|
! on entry, eps must specify a threshold for determining the
|
|
! effective rank of an over-determined linear system of equat-
|
|
! ions. 0 < eps < 1. if the number of decimal digits in the
|
|
! computer representation of a real number is q, then 10**(-q)
|
|
! is a suitable value for eps in most practical applications.
|
|
! unchanged on exit.
|
|
! nx : integer.
|
|
! unless ier=10 (in case iopt >=0), nx will contain the total
|
|
! number of knots with respect to the x-variable, of the spline
|
|
! approximation returned. if the computation mode iopt=1 is
|
|
! used, the value of nx should be left unchanged between sub-
|
|
! sequent calls.
|
|
! in case iopt=-1, the value of nx should be specified on entry
|
|
! tx : real array of dimension nmax.
|
|
! on succesful exit, this array will contain the knots of the
|
|
! spline with respect to the x-variable, i.e. the position of
|
|
! the interior knots tx(kx+2),...,tx(nx-kx-1) as well as the
|
|
! position of the additional knots tx(1)=...=tx(kx+1)=xb and
|
|
! tx(nx-kx)=...=tx(nx)=xe needed for the b-spline representat.
|
|
! if the computation mode iopt=1 is used, the values of tx(1),
|
|
! ...,tx(nx) should be left unchanged between subsequent calls.
|
|
! if the computation mode iopt=-1 is used, the values tx(kx+2),
|
|
! ...tx(nx-kx-1) must be supplied by the user, before entry.
|
|
! see also the restrictions (ier=10).
|
|
! ny : integer.
|
|
! unless ier=10 (in case iopt >=0), ny will contain the total
|
|
! number of knots with respect to the y-variable, of the spline
|
|
! approximation returned. if the computation mode iopt=1 is
|
|
! used, the value of ny should be left unchanged between sub-
|
|
! sequent calls.
|
|
! in case iopt=-1, the value of ny should be specified on entry
|
|
! ty : real array of dimension nmax.
|
|
! on succesful exit, this array will contain the knots of the
|
|
! spline with respect to the y-variable, i.e. the position of
|
|
! the interior knots ty(ky+2),...,ty(ny-ky-1) as well as the
|
|
! position of the additional knots ty(1)=...=ty(ky+1)=yb and
|
|
! ty(ny-ky)=...=ty(ny)=ye needed for the b-spline representat.
|
|
! if the computation mode iopt=1 is used, the values of ty(1),
|
|
! ...,ty(ny) should be left unchanged between subsequent calls.
|
|
! if the computation mode iopt=-1 is used, the values ty(ky+2),
|
|
! ...ty(ny-ky-1) must be supplied by the user, before entry.
|
|
! see also the restrictions (ier=10).
|
|
! c : real array of dimension at least (nxest-kx-1)*(nyest-ky-1).
|
|
! on succesful exit, c contains the coefficients of the spline
|
|
! approximation s(x,y)
|
|
! fp : real. unless ier=10, fp contains the weighted sum of
|
|
! squared residuals of the spline approximation returned.
|
|
! wrk1 : real array of dimension (lwrk1). used as workspace.
|
|
! if the computation mode iopt=1 is used the value of wrk1(1)
|
|
! should be left unchanged between subsequent calls.
|
|
! on exit wrk1(2),wrk1(3),...,wrk1(1+(nx-kx-1)*(ny-ky-1)) will
|
|
! contain the values d(i)/max(d(i)),i=1,...,(nx-kx-1)*(ny-ky-1)
|
|
! with d(i) the i-th diagonal element of the reduced triangular
|
|
! matrix for calculating the b-spline coefficients. it includes
|
|
! those elements whose square is less than eps,which are treat-
|
|
! ed as 0 in the case of presumed rank deficiency (ier<-2).
|
|
! lwrk1 : integer. on entry lwrk1 must specify the actual dimension of
|
|
! the array wrk1 as declared in the calling (sub)program.
|
|
! lwrk1 must not be too small. let
|
|
! u = nxest-kx-1, v = nyest-ky-1, km = max(kx,ky)+1,
|
|
! ne = max(nxest,nyest), bx = kx*v+ky+1, by = ky*u+kx+1,
|
|
! if(bx<=by) b1 = bx, b2 = b1+v-ky
|
|
! if(bx>by) b1 = by, b2 = b1+u-kx then
|
|
! lwrk1 >= u*v*(2+b1+b2)+2*(u+v+km*(m+ne)+ne-kx-ky)+b2+1
|
|
! wrk2 : real array of dimension (lwrk2). used as workspace, but
|
|
! only in the case a rank deficient system is encountered.
|
|
! lwrk2 : integer. on entry lwrk2 must specify the actual dimension of
|
|
! the array wrk2 as declared in the calling (sub)program.
|
|
! lwrk2 > 0 . a save upper boundfor lwrk2 = u*v*(b2+1)+b2
|
|
! where u,v and b2 are as above. if there are enough data
|
|
! points, scattered uniformly over the approximation domain
|
|
! and if the smoothing factor s is not too small, there is a
|
|
! good chance that this extra workspace is not needed. a lot
|
|
! of memory might therefore be saved by setting lwrk2=1.
|
|
! (see also ier > 10)
|
|
! iwrk : integer array of dimension (kwrk). used as workspace.
|
|
! kwrk : integer. on entry kwrk must specify the actual dimension of
|
|
! the array iwrk as declared in the calling (sub)program.
|
|
! kwrk >= m+(nxest-2*kx-1)*(nyest-2*ky-1).
|
|
! ier : integer. unless the routine detects an error, ier contains a
|
|
! non-positive value on exit, i.e.
|
|
! ier=0 : normal return. the spline returned has a residual sum of
|
|
! squares fp such that abs(fp-s)/s <= tol with tol a relat-
|
|
! ive tolerance set to 0.001 by the program.
|
|
! ier=-1 : normal return. the spline returned is an interpolating
|
|
! spline (fp=0).
|
|
! ier=-2 : normal return. the spline returned is the weighted least-
|
|
! squares polynomial of degrees kx and ky. in this extreme
|
|
! case fp gives the upper bound for the smoothing factor s.
|
|
! ier<-2 : warning. the coefficients of the spline returned have been
|
|
! computed as the minimal norm least-squares solution of a
|
|
! (numerically) rank deficient system. (-ier) gives the rank.
|
|
! especially if the rank deficiency which can be computed as
|
|
! (nx-kx-1)*(ny-ky-1)+ier, is large the results may be inac-
|
|
! curate. they could also seriously depend on the value of
|
|
! eps.
|
|
! ier=1 : error. the required storage space exceeds the available
|
|
! storage space, as specified by the parameters nxest and
|
|
! nyest.
|
|
! probably causes : nxest or nyest too small. if these param-
|
|
! eters are already large, it may also indicate that s is
|
|
! too small
|
|
! the approximation returned is the weighted least-squares
|
|
! spline according to the current set of knots.
|
|
! the parameter fp gives the corresponding weighted sum of
|
|
! squared residuals (fp>s).
|
|
! ier=2 : error. a theoretically impossible result was found during
|
|
! the iteration proces for finding a smoothing spline with
|
|
! fp = s. probably causes : s too small or badly chosen eps.
|
|
! there is an approximation returned but the corresponding
|
|
! weighted sum of squared residuals does not satisfy the
|
|
! condition abs(fp-s)/s < tol.
|
|
! ier=3 : error. the maximal number of iterations maxit (set to 20
|
|
! by the program) allowed for finding a smoothing spline
|
|
! with fp=s has been reached. probably causes : s too small
|
|
! there is an approximation returned but the corresponding
|
|
! weighted sum of squared residuals does not satisfy the
|
|
! condition abs(fp-s)/s < tol.
|
|
! ier=4 : error. no more knots can be added because the number of
|
|
! b-spline coefficients (nx-kx-1)*(ny-ky-1) already exceeds
|
|
! the number of data points m.
|
|
! probably causes : either s or m too small.
|
|
! the approximation returned is the weighted least-squares
|
|
! spline according to the current set of knots.
|
|
! the parameter fp gives the corresponding weighted sum of
|
|
! squared residuals (fp>s).
|
|
! ier=5 : error. no more knots can be added because the additional
|
|
! knot would (quasi) coincide with an old one.
|
|
! probably causes : s too small or too large a weight to an
|
|
! inaccurate data point.
|
|
! the approximation returned is the weighted least-squares
|
|
! spline according to the current set of knots.
|
|
! the parameter fp gives the corresponding weighted sum of
|
|
! squared residuals (fp>s).
|
|
! ier=10 : error. on entry, the input data are controlled on validity
|
|
! the following restrictions must be satisfied.
|
|
! -1<=iopt<=1, 1<=kx,ky<=5, m>=(kx+1)*(ky+1), nxest>=2*kx+2,
|
|
! nyest>=2*ky+2, 0<eps<1, nmax>=nxest, nmax>=nyest,
|
|
! xb<=x(i)<=xe, yb<=y(i)<=ye, w(i)>0, i=1,...,m
|
|
! lwrk1 >= u*v*(2+b1+b2)+2*(u+v+km*(m+ne)+ne-kx-ky)+b2+1
|
|
! kwrk >= m+(nxest-2*kx-1)*(nyest-2*ky-1)
|
|
! if iopt=-1: 2*kx+2<=nx<=nxest
|
|
! xb<tx(kx+2)<tx(kx+3)<...<tx(nx-kx-1)<xe
|
|
! 2*ky+2<=ny<=nyest
|
|
! yb<ty(ky+2)<ty(ky+3)<...<ty(ny-ky-1)<ye
|
|
! if iopt>=0: s>=0
|
|
! if one of these conditions is found to be violated,control
|
|
! is immediately repassed to the calling program. in that
|
|
! case there is no approximation returned.
|
|
! ier>10 : error. lwrk2 is too small, i.e. there is not enough work-
|
|
! space for computing the minimal least-squares solution of
|
|
! a rank deficient system of linear equations. ier gives the
|
|
! requested value for lwrk2. there is no approximation re-
|
|
! turned but, having saved the information contained in nx,
|
|
! ny,tx,ty,wrk1, and having adjusted the value of lwrk2 and
|
|
! the dimension of the array wrk2 accordingly, the user can
|
|
! continue at the point the program was left, by calling
|
|
! surfit with iopt=1.
|
|
!
|
|
! further comments:
|
|
! by means of the parameter s, the user can control the tradeoff
|
|
! between closeness of fit and smoothness of fit of the approximation.
|
|
! if s is too large, the spline will be too smooth and signal will be
|
|
! lost ; if s is too small the spline will pick up too much noise. in
|
|
! the extreme cases the program will return an interpolating spline if
|
|
! s=0 and the weighted least-squares polynomial (degrees kx,ky)if s is
|
|
! very large. between these extremes, a properly chosen s will result
|
|
! in a good compromise between closeness of fit and smoothness of fit.
|
|
! to decide whether an approximation, corresponding to a certain s is
|
|
! satisfactory the user is highly recommended to inspect the fits
|
|
! graphically.
|
|
! recommended values for s depend on the weights w(i). if these are
|
|
! taken as 1/d(i) with d(i) an estimate of the standard deviation of
|
|
! z(i), a good s-value should be found in the range (m-sqrt(2*m),m+
|
|
! sqrt(2*m)). if nothing is known about the statistical error in z(i)
|
|
! each w(i) can be set equal to one and s determined by trial and
|
|
! error, taking account of the comments above. the best is then to
|
|
! start with a very large value of s ( to determine the least-squares
|
|
! polynomial and the corresponding upper bound fp0 for s) and then to
|
|
! progressively decrease the value of s ( say by a factor 10 in the
|
|
! beginning, i.e. s=fp0/10, fp0/100,...and more carefully as the
|
|
! approximation shows more detail) to obtain closer fits.
|
|
! to choose s very small is strongly discouraged. this considerably
|
|
! increases computation time and memory requirements. it may also
|
|
! cause rank-deficiency (ier<-2) and endager numerical stability.
|
|
! to economize the search for a good s-value the program provides with
|
|
! different modes of computation. at the first call of the routine, or
|
|
! whenever he wants to restart with the initial set of knots the user
|
|
! must set iopt=0.
|
|
! if iopt=1 the program will continue with the set of knots found at
|
|
! the last call of the routine. this will save a lot of computation
|
|
! time if surfit is called repeatedly for different values of s.
|
|
! the number of knots of the spline returned and their location will
|
|
! depend on the value of s and on the complexity of the shape of the
|
|
! function underlying the data. if the computation mode iopt=1
|
|
! is used, the knots returned may also depend on the s-values at
|
|
! previous calls (if these were smaller). therefore, if after a number
|
|
! of trials with different s-values and iopt=1, the user can finally
|
|
! accept a fit as satisfactory, it may be worthwhile for him to call
|
|
! surfit once more with the selected value for s but now with iopt=0.
|
|
! indeed, surfit may then return an approximation of the same quality
|
|
! of fit but with fewer knots and therefore better if data reduction
|
|
! is also an important objective for the user.
|
|
! the number of knots may also depend on the upper bounds nxest and
|
|
! nyest. indeed, if at a certain stage in surfit the number of knots
|
|
! in one direction (say nx) has reached the value of its upper bound
|
|
! (nxest), then from that moment on all subsequent knots are added
|
|
! in the other (y) direction. this may indicate that the value of
|
|
! nxest is too small. on the other hand, it gives the user the option
|
|
! of limiting the number of knots the routine locates in any direction
|
|
! for example, by setting nxest=2*kx+2 (the lowest allowable value for
|
|
! nxest), the user can indicate that he wants an approximation which
|
|
! is a simple polynomial of degree kx in the variable x.
|
|
!
|
|
! other subroutines required:
|
|
! fpback,fpbspl,fpsurf,fpdisc,fpgivs,fprank,fprati,fprota,fporde
|
|
!
|
|
! references:
|
|
! dierckx p. : an algorithm for surface fitting with spline functions
|
|
! ima j. numer. anal. 1 (1981) 267-283.
|
|
! dierckx p. : an algorithm for surface fitting with spline functions
|
|
! report tw50, dept. computer science,k.u.leuven, 1980.
|
|
! dierckx p. : curve and surface fitting with splines, monographs on
|
|
! numerical analysis, oxford university press, 1993.
|
|
!
|
|
! author:
|
|
! p.dierckx
|
|
! dept. computer science, k.u. leuven
|
|
! celestijnenlaan 200a, b-3001 heverlee, belgium.
|
|
! e-mail : Paul.Dierckx@cs.kuleuven.ac.be
|
|
!
|
|
! creation date : may 1979
|
|
! latest update : march 1987
|
|
!
|
|
! ..
|
|
implicit none
|
|
! ..scalar arguments..
|
|
real(wp_) xb,xe,yb,ye,s,eps,fp
|
|
integer iopt,m,kx,ky,nxest,nyest,nmax,nx,ny,lwrk1,lwrk2,kwrk,ier
|
|
! ..array arguments..
|
|
real(wp_) x(m),y(m),z(m),w(m),tx(nmax),ty(nmax), &
|
|
c((nxest-kx-1)*(nyest-ky-1)),wrk1(lwrk1),wrk2(lwrk2)
|
|
integer iwrk(kwrk)
|
|
! ..local scalars..
|
|
real(wp_) tol
|
|
integer i,ib1,ib3,jb1,ki,kmax,km1,km2,kn,kwest,kx1,ky1,la,lbx, &
|
|
lby,lco,lf,lff,lfp,lh,lq,lsx,lsy,lwest,maxit,ncest,nest,nek, &
|
|
nminx,nminy,nmx,nmy,nreg,nrint,nxk,nyk
|
|
! ..function references..
|
|
integer max0
|
|
! ..subroutine references..
|
|
! fpsurf
|
|
! ..
|
|
! we set up the parameters tol and maxit.
|
|
maxit = 20
|
|
tol = 0.1e-02
|
|
! before starting computations a data check is made. if the input data
|
|
! are invalid,control is immediately repassed to the calling program.
|
|
ier = 10
|
|
if(eps<=0. .or. eps>=1.) return
|
|
if(kx<=0 .or. kx>5) return
|
|
kx1 = kx+1
|
|
if(ky<=0 .or. ky>5) return
|
|
ky1 = ky+1
|
|
kmax = max0(kx,ky)
|
|
km1 = kmax+1
|
|
km2 = km1+1
|
|
if(iopt<(-1) .or. iopt>1) return
|
|
if(m<(kx1*ky1)) return
|
|
nminx = 2*kx1
|
|
if(nxest<nminx .or. nxest>nmax) return
|
|
nminy = 2*ky1
|
|
if(nyest<nminy .or. nyest>nmax) return
|
|
nest = max0(nxest,nyest)
|
|
nxk = nxest-kx1
|
|
nyk = nyest-ky1
|
|
ncest = nxk*nyk
|
|
nmx = nxest-nminx+1
|
|
nmy = nyest-nminy+1
|
|
nrint = nmx+nmy
|
|
nreg = nmx*nmy
|
|
ib1 = kx*nyk+ky1
|
|
jb1 = ky*nxk+kx1
|
|
ib3 = kx1*nyk+1
|
|
if(ib1>jb1) then
|
|
ib1 = jb1
|
|
ib3 = ky1*nxk+1
|
|
end if
|
|
lwest = ncest*(2+ib1+ib3)+2*(nrint+nest*km2+m*km1)+ib3
|
|
kwest = m+nreg
|
|
if(lwrk1<lwest .or. kwrk<kwest) return
|
|
if(xb>=xe .or. yb>=ye) return
|
|
do i=1,m
|
|
if(w(i)<=0.) return
|
|
if(x(i)<xb .or. x(i)>xe) return
|
|
if(y(i)<yb .or. y(i)>ye) return
|
|
end do
|
|
if(iopt<0) then
|
|
if(nx<nminx .or. nx>nxest) return
|
|
nxk = nx-kx1
|
|
tx(kx1) = xb
|
|
tx(nxk+1) = xe
|
|
do i=kx1,nxk
|
|
if(tx(i+1)<=tx(i)) return
|
|
end do
|
|
if(ny<nminy .or. ny>nyest) return
|
|
nyk = ny-ky1
|
|
ty(ky1) = yb
|
|
ty(nyk+1) = ye
|
|
do i=ky1,nyk
|
|
if(ty(i+1)<=ty(i)) return
|
|
end do
|
|
else
|
|
if(s<0.) return
|
|
end if
|
|
ier = 0
|
|
! we partition the working space and determine the spline approximation
|
|
kn = 1
|
|
ki = kn+m
|
|
lq = 2
|
|
la = lq+ncest*ib3
|
|
lf = la+ncest*ib1
|
|
lff = lf+ncest
|
|
lfp = lff+ncest
|
|
lco = lfp+nrint
|
|
lh = lco+nrint
|
|
lbx = lh+ib3
|
|
nek = nest*km2
|
|
lby = lbx+nek
|
|
lsx = lby+nek
|
|
lsy = lsx+m*km1
|
|
call fpsurf(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,nyest, &
|
|
eps,tol,maxit,nest,km1,km2,ib1,ib3,ncest,nrint,nreg,nx,tx, &
|
|
ny,ty,c,fp,wrk1(1),wrk1(lfp),wrk1(lco),wrk1(lf),wrk1(lff), &
|
|
wrk1(la),wrk1(lq),wrk1(lbx),wrk1(lby),wrk1(lsx),wrk1(lsy), &
|
|
wrk1(lh),iwrk(ki),iwrk(kn),wrk2,lwrk2,ier)
|
|
end subroutine surfit
|
|
|
|
subroutine fpsurf(iopt,m,x,y,z,w,xb,xe,yb,ye,kxx,kyy,s,nxest, &
|
|
nyest,eta,tol,maxit,nmax,km1,km2,ib1,ib3,nc,intest,nrest, &
|
|
nx0,tx,ny0,ty,c,fp,fp0,fpint,coord,f,ff,a,q,bx,by,spx,spy,h, &
|
|
idx,nummer,wrk,lwrk,ier)
|
|
! ..
|
|
implicit none
|
|
! ..scalar arguments..
|
|
real(wp_) xb,xe,yb,ye,s,eta,tol,fp,fp0
|
|
integer iopt,m,kxx,kyy,nxest,nyest,maxit,nmax,km1,km2,ib1,ib3, &
|
|
nc,intest,nrest,nx0,ny0,lwrk,ier
|
|
! ..array arguments..
|
|
real(wp_) x(m),y(m),z(m),w(m),tx(nmax),ty(nmax),c(nc),fpint(intest), &
|
|
coord(intest),f(nc),ff(nc),a(nc,ib1),q(nc,ib3),bx(nmax,km2), &
|
|
by(nmax,km2),spx(m,km1),spy(m,km1),h(ib3),wrk(lwrk)
|
|
integer idx(nrest),nummer(m)
|
|
! ..local scalars..
|
|
real(wp_) acc,arg,cs,dmax,fac1,fac2,fpmax,fpms,f1,f2,f3,hxi,p,pinv, &
|
|
piv,p1,p2,p3,sigma,sn,sq,store,wi,x0,x1,y0,y1,zi,eps, &
|
|
rn,one,con1,con9,con4,ten
|
|
integer i,iband,iband1,iband3,iband4,ibb,ichang,ich1,ich3,ii, &
|
|
in,irot,iter,i1,i2,i3,j,jrot,jxy,j1,kx,kx1,kx2,ky,ky1,ky2,l, &
|
|
la,lf,lh,lwest,lx,ly,l1,l2,n,ncof,nk1x,nk1y,nminx,nminy,nreg, &
|
|
nrint,num,num1,nx,nxe,nxx,ny,nye,nyy,n1,rank
|
|
! ..local arrays..
|
|
real(wp_) hx(6),hy(6)
|
|
! ..function references..
|
|
! real(8) abs,sqrt
|
|
! integer min0
|
|
! ..subroutine references..
|
|
! fpback,fpbspl,fpgivs,fpdisc,fporde,fprank,fprota
|
|
! ..
|
|
! set constants
|
|
one = 0.1e+01
|
|
con1 = 0.1e0
|
|
con9 = 0.9e0
|
|
con4 = 0.4e-01
|
|
ten = 0.1e+02
|
|
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
! part 1: determination of the number of knots and their position. c
|
|
! **************************************************************** c
|
|
! given a set of knots we compute the least-squares spline sinf(x,y), c
|
|
! and the corresponding weighted sum of squared residuals fp=f(p=inf). c
|
|
! if iopt=-1 sinf(x,y) is the requested approximation. c
|
|
! if iopt=0 or iopt=1 we check whether we can accept the knots: c
|
|
! if fp <=s we will continue with the current set of knots. c
|
|
! if fp > s we will increase the number of knots and compute the c
|
|
! corresponding least-squares spline until finally fp<=s. c
|
|
! the initial choice of knots depends on the value of s and iopt. c
|
|
! if iopt=0 we first compute the least-squares polynomial of degree c
|
|
! kx in x and ky in y; nx=nminx=2*kx+2 and ny=nminy=2*ky+2. c
|
|
! fp0=f(0) denotes the corresponding weighted sum of squared c
|
|
! residuals c
|
|
! if iopt=1 we start with the knots found at the last call of the c
|
|
! routine, except for the case that s>=fp0; then we can compute c
|
|
! the least-squares polynomial directly. c
|
|
! eventually the independent variables x and y (and the corresponding c
|
|
! parameters) will be switched if this can reduce the bandwidth of the c
|
|
! system to be solved. c
|
|
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
! ichang denotes whether(1) or not(-1) the directions have been inter-
|
|
! changed.
|
|
ichang = -1
|
|
x0 = xb
|
|
x1 = xe
|
|
y0 = yb
|
|
y1 = ye
|
|
kx = kxx
|
|
ky = kyy
|
|
kx1 = kx+1
|
|
ky1 = ky+1
|
|
nxe = nxest
|
|
nye = nyest
|
|
eps = sqrt(eta)
|
|
if(iopt<0) go to 20
|
|
! calculation of acc, the absolute tolerance for the root of f(p)=s.
|
|
acc = tol*s
|
|
if(iopt==0) go to 10
|
|
if(fp0>s) go to 20
|
|
! initialization for the least-squares polynomial.
|
|
10 continue
|
|
nminx = 2*kx1
|
|
nminy = 2*ky1
|
|
nx = nminx
|
|
ny = nminy
|
|
ier = -2
|
|
go to 30
|
|
20 continue
|
|
nx = nx0
|
|
ny = ny0
|
|
! main loop for the different sets of knots. m is a save upper bound
|
|
! for the number of trials.
|
|
30 continue
|
|
do iter=1,m
|
|
! find the position of the additional knots which are needed for the
|
|
! b-spline representation of s(x,y).
|
|
l = nx
|
|
do i=1,kx1
|
|
tx(i) = x0
|
|
tx(l) = x1
|
|
l = l-1
|
|
end do
|
|
l = ny
|
|
do i=1,ky1
|
|
ty(i) = y0
|
|
ty(l) = y1
|
|
l = l-1
|
|
end do
|
|
! find nrint, the total number of knot intervals and nreg, the number
|
|
! of panels in which the approximation domain is subdivided by the
|
|
! intersection of knots.
|
|
nxx = nx-2*kx1+1
|
|
nyy = ny-2*ky1+1
|
|
nrint = nxx+nyy
|
|
nreg = nxx*nyy
|
|
! find the bandwidth of the observation matrix a.
|
|
! if necessary, interchange the variables x and y, in order to obtain
|
|
! a minimal bandwidth.
|
|
iband1 = kx*(ny-ky1)+ky
|
|
l = ky*(nx-kx1)+kx
|
|
if(iband1>l) then
|
|
iband1 = l
|
|
ichang = -ichang
|
|
do i=1,m
|
|
store = x(i)
|
|
x(i) = y(i)
|
|
y(i) = store
|
|
end do
|
|
store = x0
|
|
x0 = y0
|
|
y0 = store
|
|
store = x1
|
|
x1 = y1
|
|
y1 = store
|
|
n = min0(nx,ny)
|
|
do i=1,n
|
|
store = tx(i)
|
|
tx(i) = ty(i)
|
|
ty(i) = store
|
|
end do
|
|
n1 = n+1
|
|
if(nx<ny) then
|
|
do i=n1,ny
|
|
tx(i) = ty(i)
|
|
end do
|
|
else if (nx>ny) then
|
|
do i=n1,nx
|
|
ty(i) = tx(i)
|
|
end do
|
|
end if
|
|
l = nx
|
|
nx = ny
|
|
ny = l
|
|
l = nxe
|
|
nxe = nye
|
|
nye = l
|
|
l = nxx
|
|
nxx = nyy
|
|
nyy = l
|
|
l = kx
|
|
kx = ky
|
|
ky = l
|
|
kx1 = kx+1
|
|
ky1 = ky+1
|
|
end if
|
|
iband = iband1+1
|
|
! arrange the data points according to the panel they belong to.
|
|
call fporde(x,y,m,kx,ky,tx,nx,ty,ny,nummer,idx,nreg)
|
|
! find ncof, the number of b-spline coefficients.
|
|
nk1x = nx-kx1
|
|
nk1y = ny-ky1
|
|
ncof = nk1x*nk1y
|
|
! initialize the observation matrix a.
|
|
do i=1,ncof
|
|
f(i) = 0.
|
|
do j=1,iband
|
|
a(i,j) = 0.
|
|
end do
|
|
end do
|
|
! initialize the sum of squared residuals.
|
|
fp = 0.
|
|
! fetch the data points in the new order. main loop for the
|
|
! different panels.
|
|
do num=1,nreg
|
|
! fix certain constants for the current panel; jrot records the column
|
|
! number of the first non-zero element in a row of the observation
|
|
! matrix according to a data point of the panel.
|
|
num1 = num-1
|
|
lx = num1/nyy
|
|
l1 = lx+kx1
|
|
ly = num1-lx*nyy
|
|
l2 = ly+ky1
|
|
jrot = lx*nk1y+ly
|
|
! test whether there are still data points in the panel.
|
|
in = idx(num)
|
|
do
|
|
if(in==0) exit
|
|
! fetch a new data point.
|
|
wi = w(in)
|
|
zi = z(in)*wi
|
|
! evaluate for the x-direction, the (kx+1) non-zero b-splines at x(in).
|
|
call fpbspl(tx,nx,kx,x(in),l1,hx)
|
|
! evaluate for the y-direction, the (ky+1) non-zero b-splines at y(in).
|
|
call fpbspl(ty,ny,ky,y(in),l2,hy)
|
|
! store the value of these b-splines in spx and spy respectively.
|
|
do i=1,kx1
|
|
spx(in,i) = hx(i)
|
|
end do
|
|
do i=1,ky1
|
|
spy(in,i) = hy(i)
|
|
end do
|
|
! initialize the new row of observation matrix.
|
|
do i=1,iband
|
|
h(i) = 0.
|
|
end do
|
|
! calculate the non-zero elements of the new row by making the cross
|
|
! products of the non-zero b-splines in x- and y-direction.
|
|
i1 = 0
|
|
do i=1,kx1
|
|
hxi = hx(i)
|
|
j1 = i1
|
|
do j=1,ky1
|
|
j1 = j1+1
|
|
h(j1) = hxi*hy(j)*wi
|
|
end do
|
|
i1 = i1+nk1y
|
|
end do
|
|
! rotate the row into triangle by givens transformations .
|
|
irot = jrot
|
|
do i=1,iband
|
|
irot = irot+1
|
|
piv = h(i)
|
|
if(piv==0.) cycle
|
|
! calculate the parameters of the givens transformation.
|
|
call fpgivs(piv,a(irot,1),cs,sn)
|
|
! apply that transformation to the right hand side.
|
|
call fprota(cs,sn,zi,f(irot))
|
|
if(i==iband) exit
|
|
! apply that transformation to the left hand side.
|
|
i2 = 1
|
|
i3 = i+1
|
|
do j=i3,iband
|
|
i2 = i2+1
|
|
call fprota(cs,sn,h(j),a(irot,i2))
|
|
end do
|
|
end do
|
|
! add the contribution of the row to the sum of squares of residual
|
|
! right hand sides.
|
|
fp = fp+zi**2
|
|
! find the number of the next data point in the panel.
|
|
in = nummer(in)
|
|
end do
|
|
end do
|
|
! find dmax, the maximum value for the diagonal elements in the reduced
|
|
! triangle.
|
|
dmax = 0.
|
|
do i=1,ncof
|
|
if(a(i,1)<=dmax) cycle
|
|
dmax = a(i,1)
|
|
end do
|
|
! check whether the observation matrix is rank deficient.
|
|
sigma = eps*dmax
|
|
do i=1,ncof
|
|
if(a(i,1)<=sigma) go to 280
|
|
end do
|
|
! backward substitution in case of full rank.
|
|
call fpback(a,f,ncof,iband,c,nc)
|
|
rank = ncof
|
|
do i=1,ncof
|
|
q(i,1) = a(i,1)/dmax
|
|
end do
|
|
go to 300
|
|
! in case of rank deficiency, find the minimum norm solution.
|
|
! check whether there is sufficient working space
|
|
280 continue
|
|
lwest = ncof*iband+ncof+iband
|
|
if(lwrk<lwest) go to 780
|
|
do i=1,ncof
|
|
ff(i) = f(i)
|
|
do j=1,iband
|
|
q(i,j) = a(i,j)
|
|
end do
|
|
end do
|
|
lf =1
|
|
lh = lf+ncof
|
|
la = lh+iband
|
|
call fprank(q,ff,ncof,iband,nc,sigma,c,sq,rank,wrk(la), &
|
|
wrk(lf),wrk(lh))
|
|
do i=1,ncof
|
|
q(i,1) = q(i,1)/dmax
|
|
end do
|
|
! add to the sum of squared residuals, the contribution of reducing
|
|
! the rank.
|
|
fp = fp+sq
|
|
300 continue
|
|
if(ier==(-2)) fp0 = fp
|
|
! test whether the least-squares spline is an acceptable solution.
|
|
if(iopt<0) go to 820
|
|
fpms = fp-s
|
|
if(abs(fpms)<=acc) then
|
|
if(fp>0.0_wp_) then
|
|
go to 820
|
|
else
|
|
go to 815
|
|
end if
|
|
end if
|
|
! test whether we can accept the choice of knots.
|
|
if(fpms<0.) exit
|
|
! test whether we cannot further increase the number of knots.
|
|
if(ncof>m) go to 790
|
|
ier = 0
|
|
! search where to add a new knot.
|
|
! find for each interval the sum of squared residuals fpint for the
|
|
! data points having the coordinate belonging to that knot interval.
|
|
! calculate also coord which is the same sum, weighted by the position
|
|
! of the data points considered.
|
|
do i=1,nrint
|
|
fpint(i) = 0.
|
|
coord(i) = 0.
|
|
end do
|
|
do num=1,nreg
|
|
num1 = num-1
|
|
lx = num1/nyy
|
|
l1 = lx+1
|
|
ly = num1-lx*nyy
|
|
l2 = ly+1+nxx
|
|
jrot = lx*nk1y+ly
|
|
in = idx(num)
|
|
do
|
|
if(in==0) exit
|
|
store = 0.
|
|
i1 = jrot
|
|
do i=1,kx1
|
|
hxi = spx(in,i)
|
|
j1 = i1
|
|
do j=1,ky1
|
|
j1 = j1+1
|
|
store = store+hxi*spy(in,j)*c(j1)
|
|
end do
|
|
i1 = i1+nk1y
|
|
end do
|
|
store = (w(in)*(z(in)-store))**2
|
|
fpint(l1) = fpint(l1)+store
|
|
coord(l1) = coord(l1)+store*x(in)
|
|
fpint(l2) = fpint(l2)+store
|
|
coord(l2) = coord(l2)+store*y(in)
|
|
in = nummer(in)
|
|
end do
|
|
end do
|
|
! find the interval for which fpint is maximal on the condition that
|
|
! there still can be added a knot.
|
|
do
|
|
l = 0
|
|
fpmax = 0.
|
|
l1 = 1
|
|
l2 = nrint
|
|
if(nx==nxe) l1 = nxx+1
|
|
if(ny==nye) l2 = nxx
|
|
if(l1>l2) go to 810
|
|
do i=l1,l2
|
|
if(fpmax>=fpint(i)) cycle
|
|
l = i
|
|
fpmax = fpint(i)
|
|
end do
|
|
! test whether we cannot further increase the number of knots.
|
|
if(l==0) go to 785
|
|
! calculate the position of the new knot.
|
|
arg = coord(l)/fpint(l)
|
|
! test in what direction the new knot is going to be added.
|
|
if(l<=nxx) then
|
|
! addition in the x-direction.
|
|
jxy = l+kx1
|
|
fpint(l) = 0.
|
|
fac1 = tx(jxy)-arg
|
|
fac2 = arg-tx(jxy-1)
|
|
if(fac1>(ten*fac2) .or. fac2>(ten*fac1)) cycle
|
|
j = nx
|
|
do i=jxy,nx
|
|
tx(j+1) = tx(j)
|
|
j = j-1
|
|
end do
|
|
tx(jxy) = arg
|
|
nx = nx+1
|
|
else
|
|
! addition in the y-direction.
|
|
jxy = l+ky1-nxx
|
|
fpint(l) = 0.
|
|
fac1 = ty(jxy)-arg
|
|
fac2 = arg-ty(jxy-1)
|
|
if(fac1>(ten*fac2) .or. fac2>(ten*fac1)) cycle
|
|
j = ny
|
|
do i=jxy,ny
|
|
ty(j+1) = ty(j)
|
|
j = j-1
|
|
end do
|
|
ty(jxy) = arg
|
|
ny = ny+1
|
|
end if
|
|
exit
|
|
end do
|
|
! restart the computations with the new set of knots.
|
|
end do
|
|
! test whether the least-squares polynomial is a solution of our
|
|
! approximation problem.
|
|
if(ier==(-2)) go to 830
|
|
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
! part 2: determination of the smoothing spline sp(x,y) c
|
|
! ***************************************************** c
|
|
! we have determined the number of knots and their position. we now c
|
|
! compute the b-spline coefficients of the smoothing spline sp(x,y). c
|
|
! the observation matrix a is extended by the rows of a matrix, c
|
|
! expressing that sp(x,y) must be a polynomial of degree kx in x and c
|
|
! ky in y. the corresponding weights of these additional rows are set c
|
|
! to 1./p. iteratively we than have to determine the value of p c
|
|
! such that f(p)=sum((w(i)*(z(i)-sp(x(i),y(i))))**2) be = s. c
|
|
! we already know that the least-squares polynomial corresponds to c
|
|
! p=0 and that the least-squares spline corresponds to p=infinity. c
|
|
! the iteration process which is proposed here makes use of rational c
|
|
! interpolation. since f(p) is a convex and strictly decreasing c
|
|
! function of p, it can be approximated by a rational function r(p)= c
|
|
! (u*p+v)/(p+w). three values of p(p1,p2,p3) with corresponding values c
|
|
! of f(p) (f1=f(p1)-s,f2=f(p2)-s,f3=f(p3)-s) are used to calculate the c
|
|
! new value of p such that r(p)=s. convergence is guaranteed by taking c
|
|
! f1 > 0 and f3 < 0. c
|
|
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
kx2 = kx1+1
|
|
! test whether there are interior knots in the x-direction.
|
|
! and
|
|
! evaluate the discotinuity jumps of the kx-th order derivative of
|
|
! the b-splines at the knots tx(l),l=kx+2,...,nx-kx-1.
|
|
if(nk1x/=kx1) call fpdisc(tx,nx,kx2,bx,nmax)
|
|
ky2 = ky1 + 1
|
|
! test whether there are interior knots in the y-direction.
|
|
! and
|
|
! evaluate the discontinuity jumps of the ky-th order derivative of
|
|
! the b-splines at the knots ty(l),l=ky+2,...,ny-ky-1.
|
|
if(nk1y/=ky1) call fpdisc(ty,ny,ky2,by,nmax)
|
|
! initial value for p.
|
|
p1 = 0.
|
|
f1 = fp0-s
|
|
p3 = -one
|
|
f3 = fpms
|
|
p = 0.
|
|
do i=1,ncof
|
|
p = p+a(i,1)
|
|
end do
|
|
rn = ncof
|
|
p = rn/p
|
|
! find the bandwidth of the extended observation matrix.
|
|
iband3 = kx1*nk1y
|
|
iband4 = iband3 +1
|
|
ich1 = 0
|
|
ich3 = 0
|
|
! iteration process to find the root of f(p)=s.
|
|
do 770 iter=1,maxit
|
|
pinv = one/p
|
|
! store the triangularized observation matrix into q.
|
|
do i=1,ncof
|
|
ff(i) = f(i)
|
|
do j=1,iband
|
|
q(i,j) = a(i,j)
|
|
end do
|
|
ibb = iband+1
|
|
do j=ibb,iband4
|
|
q(i,j) = 0.
|
|
end do
|
|
end do
|
|
if(nk1y/=ky1) then
|
|
! extend the observation matrix with the rows of a matrix, expressing
|
|
! that for x=cst. sp(x,y) must be a polynomial in y of degree ky.
|
|
do i=ky2,nk1y
|
|
ii = i-ky1
|
|
do j=1,nk1x
|
|
! initialize the new row.
|
|
do l=1,iband
|
|
h(l) = 0.
|
|
end do
|
|
! fill in the non-zero elements of the row. jrot records the column
|
|
! number of the first non-zero element in the row.
|
|
do l=1,ky2
|
|
h(l) = by(ii,l)*pinv
|
|
end do
|
|
zi = 0.
|
|
jrot = (j-1)*nk1y+ii
|
|
! rotate the new row into triangle by givens transformations without
|
|
! square roots.
|
|
do irot=jrot,ncof
|
|
piv = h(1)
|
|
i2 = min0(iband1,ncof-irot)
|
|
if(piv/=0.) then
|
|
! calculate the parameters of the givens transformation.
|
|
call fpgivs(piv,q(irot,1),cs,sn)
|
|
! apply that givens transformation to the right hand side.
|
|
call fprota(cs,sn,zi,ff(irot))
|
|
if(i2==0) exit
|
|
! apply that givens transformation to the left hand side.
|
|
do l=1,i2
|
|
l1 = l+1
|
|
call fprota(cs,sn,h(l1),q(irot,l1))
|
|
end do
|
|
else
|
|
if(i2<=0) exit
|
|
end if
|
|
do l=1,i2
|
|
h(l) = h(l+1)
|
|
end do
|
|
h(i2+1) = 0.
|
|
end do
|
|
end do
|
|
end do
|
|
end if
|
|
if(nk1x/=kx1) then
|
|
! extend the observation matrix with the rows of a matrix expressing
|
|
! that for y=cst. sp(x,y) must be a polynomial in x of degree kx.
|
|
do i=kx2,nk1x
|
|
ii = i-kx1
|
|
do j=1,nk1y
|
|
! initialize the new row
|
|
do l=1,iband4
|
|
h(l) = 0.
|
|
end do
|
|
! fill in the non-zero elements of the row. jrot records the column
|
|
! number of the first non-zero element in the row.
|
|
j1 = 1
|
|
do l=1,kx2
|
|
h(j1) = bx(ii,l)*pinv
|
|
j1 = j1+nk1y
|
|
end do
|
|
zi = 0.
|
|
jrot = (i-kx2)*nk1y+j
|
|
! rotate the new row into triangle by givens transformations .
|
|
do irot=jrot,ncof
|
|
piv = h(1)
|
|
i2 = min0(iband3,ncof-irot)
|
|
if(piv/=0.) then
|
|
! calculate the parameters of the givens transformation.
|
|
call fpgivs(piv,q(irot,1),cs,sn)
|
|
! apply that givens transformation to the right hand side.
|
|
call fprota(cs,sn,zi,ff(irot))
|
|
if(i2==0) exit
|
|
! apply that givens transformation to the left hand side.
|
|
do l=1,i2
|
|
l1 = l+1
|
|
call fprota(cs,sn,h(l1),q(irot,l1))
|
|
end do
|
|
else
|
|
if(i2<=0) exit
|
|
end if
|
|
do l=1,i2
|
|
h(l) = h(l+1)
|
|
end do
|
|
h(i2+1) = 0.
|
|
end do
|
|
end do
|
|
end do
|
|
end if
|
|
! find dmax, the maximum value for the diagonal elements in the
|
|
! reduced triangle.
|
|
dmax = 0.
|
|
do i=1,ncof
|
|
if(q(i,1)<=dmax) cycle
|
|
dmax = q(i,1)
|
|
end do
|
|
! check whether the matrix is rank deficient.
|
|
sigma = eps*dmax
|
|
do i=1,ncof
|
|
if(q(i,1)<=sigma) go to 670
|
|
end do
|
|
! backward substitution in case of full rank.
|
|
call fpback(q,ff,ncof,iband4,c,nc)
|
|
rank = ncof
|
|
go to 675
|
|
! in case of rank deficiency, find the minimum norm solution.
|
|
670 continue
|
|
lwest = ncof*iband4+ncof+iband4
|
|
if(lwrk<lwest) go to 780
|
|
lf = 1
|
|
lh = lf+ncof
|
|
la = lh+iband4
|
|
call fprank(q,ff,ncof,iband4,nc,sigma,c,sq,rank,wrk(la), &
|
|
wrk(lf),wrk(lh))
|
|
675 continue
|
|
do i=1,ncof
|
|
q(i,1) = q(i,1)/dmax
|
|
end do
|
|
! compute f(p).
|
|
fp = 0.
|
|
do num = 1,nreg
|
|
num1 = num-1
|
|
lx = num1/nyy
|
|
ly = num1-lx*nyy
|
|
jrot = lx*nk1y+ly
|
|
in = idx(num)
|
|
do
|
|
if(in==0) exit
|
|
store = 0.
|
|
i1 = jrot
|
|
do i=1,kx1
|
|
hxi = spx(in,i)
|
|
j1 = i1
|
|
do j=1,ky1
|
|
j1 = j1+1
|
|
store = store+hxi*spy(in,j)*c(j1)
|
|
end do
|
|
i1 = i1+nk1y
|
|
end do
|
|
fp = fp+(w(in)*(z(in)-store))**2
|
|
in = nummer(in)
|
|
end do
|
|
end do
|
|
! test whether the approximation sp(x,y) is an acceptable solution.
|
|
fpms = fp-s
|
|
if(abs(fpms)<=acc) go to 820
|
|
! test whether the maximum allowable number of iterations has been
|
|
! reached.
|
|
if(iter==maxit) go to 795
|
|
! carry out one more step of the iteration process.
|
|
p2 = p
|
|
f2 = fpms
|
|
if(ich3/=0) go to 740
|
|
if((f2-f3)>acc) go to 730
|
|
! our initial choice of p is too large.
|
|
p3 = p2
|
|
f3 = f2
|
|
p = p*con4
|
|
if(p<=p1) p = p1*con9 + p2*con1
|
|
go to 770
|
|
730 if(f2<0.) ich3 = 1
|
|
740 if(ich1/=0) go to 760
|
|
if((f1-f2)>acc) go to 750
|
|
! our initial choice of p is too small
|
|
p1 = p2
|
|
f1 = f2
|
|
p = p/con4
|
|
if(p3<0.) go to 770
|
|
if(p>=p3) p = p2*con1 + p3*con9
|
|
go to 770
|
|
750 if(f2>0.) ich1 = 1
|
|
! test whether the iteration process proceeds as theoretically
|
|
! expected.
|
|
760 if(f2>=f1 .or. f2<=f3) go to 800
|
|
! find the new value of p.
|
|
p = fprati(p1,f1,p2,f2,p3,f3)
|
|
770 continue
|
|
! error codes and messages.
|
|
780 ier = lwest
|
|
go to 830
|
|
785 ier = 5
|
|
go to 830
|
|
790 ier = 4
|
|
go to 830
|
|
795 ier = 3
|
|
go to 830
|
|
800 ier = 2
|
|
go to 830
|
|
810 ier = 1
|
|
go to 830
|
|
815 ier = -1
|
|
fp = 0.
|
|
820 if(ncof/=rank) ier = -rank
|
|
! test whether x and y are in the original order.
|
|
830 if(ichang<0) go to 930
|
|
! if not, interchange x and y once more.
|
|
l1 = 1
|
|
do i=1,nk1x
|
|
l2 = i
|
|
do j=1,nk1y
|
|
f(l2) = c(l1)
|
|
l1 = l1+1
|
|
l2 = l2+nk1x
|
|
end do
|
|
end do
|
|
do i=1,ncof
|
|
c(i) = f(i)
|
|
end do
|
|
do i=1,m
|
|
store = x(i)
|
|
x(i) = y(i)
|
|
y(i) = store
|
|
end do
|
|
n = min0(nx,ny)
|
|
do i=1,n
|
|
store = tx(i)
|
|
tx(i) = ty(i)
|
|
ty(i) = store
|
|
end do
|
|
n1 = n+1
|
|
if(nx<ny) then
|
|
do i=n1,ny
|
|
tx(i) = ty(i)
|
|
end do
|
|
else if(nx>ny) then
|
|
do i=n1,nx
|
|
ty(i) = tx(i)
|
|
end do
|
|
end if
|
|
l = nx
|
|
nx = ny
|
|
ny = l
|
|
930 continue
|
|
if(iopt>=0) then
|
|
nx0 = nx
|
|
ny0 = ny
|
|
end if
|
|
end subroutine fpsurf
|
|
|
|
subroutine fpback(a,z,n,k,c,nest)
|
|
! subroutine fpback calculates the solution of the system of
|
|
! equations a*c = z with a a n x n upper triangular matrix
|
|
! of bandwidth k.
|
|
! ..
|
|
implicit none
|
|
! arguments
|
|
integer, intent(in) :: n, k, nest
|
|
real(wp_), intent(in) :: a(nest,k), z(n)
|
|
real(wp_), intent(inout) :: c(n)
|
|
! local variables
|
|
real(wp_) :: store
|
|
integer :: i, i1, j, k1, l, m
|
|
! ..
|
|
k1 = k-1
|
|
c(n) = z(n)/a(n,1)
|
|
i = n-1
|
|
if(i==0) return
|
|
do j=2,n
|
|
store = z(i)
|
|
i1 = k1
|
|
if(j<=k1) i1 = j-1
|
|
m = i
|
|
do l=1,i1
|
|
m = m+1
|
|
store = store-c(m)*a(i,l+1)
|
|
end do
|
|
c(i) = store/a(i,1)
|
|
i = i-1
|
|
end do
|
|
end subroutine fpback
|
|
|
|
subroutine fpbisp(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wx,wy,lx,ly)
|
|
implicit none
|
|
! arguments
|
|
integer, intent(in) :: nx, ny, kx, ky, mx, my
|
|
integer, intent(out) :: lx(mx), ly(my)
|
|
real(wp_), intent(in) :: tx(nx), ty(ny), c((nx-kx-1)*(ny-ky-1)), &
|
|
x(mx), y(my)
|
|
real(wp_), intent(out) :: z(mx*my)
|
|
real(wp_), intent(out) :: wx(mx,kx+1), wy(my,ky+1)
|
|
! local variables
|
|
integer :: kx1, ky1, l, l1, l2, m, nkx1, nky1, i, j, i1, j1
|
|
real(wp_) :: arg, sp, tb, te, h(6)
|
|
! ..subroutine references..
|
|
! fpbspl
|
|
! ..
|
|
kx1 = kx+1
|
|
nkx1 = nx-kx1
|
|
tb = tx(kx1)
|
|
te = tx(nkx1+1)
|
|
l = kx1
|
|
l1 = l+1
|
|
do i=1,mx
|
|
arg = x(i)
|
|
if(arg<tb) arg = tb
|
|
if(arg>te) arg = te
|
|
do
|
|
if(arg<tx(l1) .or. l==nkx1) exit
|
|
l = l1
|
|
l1 = l+1
|
|
end do
|
|
call fpbspl(tx,nx,kx,arg,l,h)
|
|
lx(i) = l-kx1
|
|
do j=1,kx1
|
|
wx(i,j) = h(j)
|
|
end do
|
|
end do
|
|
ky1 = ky+1
|
|
nky1 = ny-ky1
|
|
tb = ty(ky1)
|
|
te = ty(nky1+1)
|
|
l = ky1
|
|
l1 = l+1
|
|
do i=1,my
|
|
arg = y(i)
|
|
if(arg<tb) arg = tb
|
|
if(arg>te) arg = te
|
|
do
|
|
if(arg<ty(l1) .or. l==nky1) exit
|
|
l = l1
|
|
l1 = l+1
|
|
end do
|
|
call fpbspl(ty,ny,ky,arg,l,h)
|
|
ly(i) = l-ky1
|
|
do j=1,ky1
|
|
wy(i,j) = h(j)
|
|
end do
|
|
end do
|
|
m = 0
|
|
do i=1,mx
|
|
l = lx(i)*nky1
|
|
do i1=1,kx1
|
|
h(i1) = wx(i,i1)
|
|
end do
|
|
do j=1,my
|
|
l1 = l+ly(j)
|
|
sp = 0.0_wp_
|
|
do i1=1,kx1
|
|
l2 = l1
|
|
do j1=1,ky1
|
|
l2 = l2+1
|
|
sp = sp+c(l2)*h(i1)*wy(j,j1)
|
|
end do
|
|
l1 = l1+nky1
|
|
end do
|
|
m = m+1
|
|
z(m) = sp
|
|
end do
|
|
end do
|
|
end subroutine fpbisp
|
|
|
|
subroutine fpbspl(t,n,k,x,l,h)
|
|
! subroutine fpbspl evaluates the (k+1) non-zero b-splines of
|
|
! degree k at t(l) <= x < t(l+1) using the stable recurrence
|
|
! relation of de boor and cox.
|
|
! ..
|
|
implicit none
|
|
! arguments
|
|
integer, intent(in) :: n, k, l
|
|
real(wp_), intent(in) :: x, t(n)
|
|
real(wp_), intent(out) :: h(6)
|
|
! local variables
|
|
integer :: i, j, li, lj
|
|
real(wp_) :: f, hh(5)
|
|
! parameters
|
|
real(wp_), parameter :: one = 0.1e1_wp_
|
|
! ..
|
|
h(1) = one
|
|
do j=1,k
|
|
do i=1,j
|
|
hh(i) = h(i)
|
|
end do
|
|
h(1) = 0.0_wp_
|
|
do i=1,j
|
|
li = l+i
|
|
lj = li-j
|
|
f = hh(i)/(t(li)-t(lj))
|
|
h(i) = h(i)+f*(t(li)-x)
|
|
h(i+1) = f*(x-t(lj))
|
|
end do
|
|
end do
|
|
end subroutine fpbspl
|
|
|
|
subroutine fpchec(x,m,t,n,k,ier)
|
|
! subroutine fpchec verifies the number and the position of the knots
|
|
! t(j),j=1,2,...,n of a spline of degree k, in relation to the number
|
|
! and the position of the data points x(i),i=1,2,...,m. if all of the
|
|
! following conditions are fulfilled, the error parameter ier is set
|
|
! to zero. if one of the conditions is violated ier is set to ten.
|
|
! 1) k+1 <= n-k-1 <= m
|
|
! 2) t(1) <= t(2) <= ... <= t(k+1)
|
|
! t(n-k) <= t(n-k+1) <= ... <= t(n)
|
|
! 3) t(k+1) < t(k+2) < ... < t(n-k)
|
|
! 4) t(k+1) <= x(i) <= t(n-k)
|
|
! 5) the conditions specified by schoenberg and whitney must hold
|
|
! for at least one subset of data points, i.e. there must be a
|
|
! subset of data points y(j) such that
|
|
! t(j) < y(j) < t(j+k+1), j=1,2,...,n-k-1
|
|
! ..
|
|
implicit none
|
|
! arguments
|
|
integer, intent(in) :: m, n, k
|
|
real(wp_), intent(in) :: x(m), t(n)
|
|
integer, intent(out) :: ier
|
|
! local variables
|
|
integer :: i, j, k1, k2, l, nk1, nk2, nk3
|
|
real(wp_) :: tj, tl
|
|
! ..
|
|
k1 = k+1
|
|
k2 = k1+1
|
|
nk1 = n-k1
|
|
nk2 = nk1+1
|
|
ier = 10
|
|
! check condition no 1
|
|
if(nk1<k1 .or. nk1>m) return
|
|
! check condition no 2
|
|
j = n
|
|
do i=1,k
|
|
if(t(i)>t(i+1) .or. t(j)<t(j-1)) return
|
|
j = j-1
|
|
end do
|
|
! check condition no 3
|
|
do i=k2,nk2
|
|
if(t(i)<=t(i-1)) return
|
|
end do
|
|
! check condition no 4
|
|
if(x(1)<t(k1) .or. x(m)>t(nk2)) return
|
|
! check condition no 5
|
|
if(x(1)>=t(k2) .or. x(m)<=t(nk1)) return
|
|
i = 1
|
|
l = k2
|
|
nk3 = nk1-1
|
|
do j=2,nk3
|
|
tj = t(j)
|
|
l = l+1
|
|
tl = t(l)
|
|
do
|
|
i = i+1
|
|
if(i>=m) return
|
|
if(x(i)>tj) exit
|
|
end do
|
|
if(x(i)>=tl) return
|
|
end do
|
|
ier = 0
|
|
end subroutine fpchec
|
|
|
|
subroutine fpdisc(t,n,k2,b,nest)
|
|
! subroutine fpdisc calculates the discontinuity jumps of the kth
|
|
! derivative of the b-splines of degree k at the knots t(k+2)..t(n-k-1)
|
|
implicit none
|
|
! arguments
|
|
integer, intent(in) :: n, k2, nest
|
|
real(wp_), intent(in) :: t(n)
|
|
real(wp_), intent(out) :: b(nest,k2)
|
|
! local variables
|
|
real(wp_) :: an, fac, prod
|
|
integer :: i, ik, j, jk, k, k1, l, lj, lk, lmk, lp, nk1, nrint
|
|
real(wp_), dimension(12) :: h
|
|
! ..
|
|
k1 = k2-1
|
|
k = k1-1
|
|
nk1 = n-k1
|
|
nrint = nk1-k
|
|
an = nrint
|
|
fac = an/(t(nk1+1)-t(k1))
|
|
do l=k2,nk1
|
|
lmk = l-k1
|
|
do j=1,k1
|
|
ik = j+k1
|
|
lj = l+j
|
|
lk = lj-k2
|
|
h(j) = t(l)-t(lk)
|
|
h(ik) = t(l)-t(lj)
|
|
end do
|
|
lp = lmk
|
|
do j=1,k2
|
|
jk = j
|
|
prod = h(j)
|
|
do i=1,k
|
|
jk = jk+1
|
|
prod = prod*h(jk)*fac
|
|
end do
|
|
lk = lp+k1
|
|
b(lmk,j) = (t(lk)-t(lp))/prod
|
|
lp = lp+1
|
|
end do
|
|
end do
|
|
end subroutine fpdisc
|
|
|
|
subroutine fprank(a,f,n,m,na,tol,c,sq,rank,aa,ff,h)
|
|
! subroutine fprank finds the minimum norm solution of a least-
|
|
! squares problem in case of rank deficiency.
|
|
!
|
|
! input parameters:
|
|
! a : array, which contains the non-zero elements of the observation
|
|
! matrix after triangularization by givens transformations.
|
|
! f : array, which contains the transformed right hand side.
|
|
! n : integer,wich contains the dimension of a.
|
|
! m : integer, which denotes the bandwidth of a.
|
|
! tol : real value, giving a threshold to determine the rank of a.
|
|
!
|
|
! output parameters:
|
|
! c : array, which contains the minimum norm solution.
|
|
! sq : real value, giving the contribution of reducing the rank
|
|
! to the sum of squared residuals.
|
|
! rank : integer, which contains the rank of matrix a.
|
|
!
|
|
implicit none
|
|
! ..scalar arguments..
|
|
integer n,m,na,rank
|
|
real(wp_) tol,sq
|
|
! ..array arguments..
|
|
real(wp_) a(na,m),f(n),c(n),aa(n,m),ff(n),h(m)
|
|
! ..local scalars..
|
|
integer i,ii,ij,i1,i2,j,jj,j1,j2,j3,k,kk,m1,nl
|
|
real(wp_) cs,fac,piv,sn,yi
|
|
real(wp_) store,stor1,stor2,stor3
|
|
! ..function references..
|
|
integer min0
|
|
! ..subroutine references..
|
|
! fpgivs,fprota
|
|
! ..
|
|
m1 = m-1
|
|
! the rank deficiency nl is considered to be the number of sufficient
|
|
! small diagonal elements of a.
|
|
nl = 0
|
|
sq = 0.
|
|
do i=1,n
|
|
if(a(i,1)>tol) cycle
|
|
! if a sufficient small diagonal element is found, we put it to
|
|
! zero. the remainder of the row corresponding to that zero diagonal
|
|
! element is then rotated into triangle by givens rotations .
|
|
! the rank deficiency is increased by one.
|
|
nl = nl+1
|
|
if(i==n) cycle
|
|
yi = f(i)
|
|
do j=1,m1
|
|
h(j) = a(i,j+1)
|
|
end do
|
|
h(m) = 0.
|
|
i1 = i+1
|
|
do ii=i1,n
|
|
i2 = min0(n-ii,m1)
|
|
piv = h(1)
|
|
if(piv/=0.) then
|
|
call fpgivs(piv,a(ii,1),cs,sn)
|
|
call fprota(cs,sn,yi,f(ii))
|
|
if(i2==0) exit
|
|
do j=1,i2
|
|
j1 = j+1
|
|
call fprota(cs,sn,h(j1),a(ii,j1))
|
|
h(j) = h(j1)
|
|
end do
|
|
else
|
|
if(i2==0) exit
|
|
do j=1,i2
|
|
h(j) = h(j+1)
|
|
end do
|
|
end if
|
|
h(i2+1) = 0.
|
|
end do
|
|
! add to the sum of squared residuals the contribution of deleting
|
|
! the row with small diagonal element.
|
|
sq = sq+yi**2
|
|
end do
|
|
! rank denotes the rank of a.
|
|
rank = n-nl
|
|
! let b denote the (rank*n) upper trapezoidal matrix which can be
|
|
! obtained from the (n*n) upper triangular matrix a by deleting
|
|
! the rows and interchanging the columns corresponding to a zero
|
|
! diagonal element. if this matrix is factorized using givens
|
|
! transformations as b = (r) (u) where
|
|
! r is a (rank*rank) upper triangular matrix,
|
|
! u is a (rank*n) orthonormal matrix
|
|
! then the minimal least-squares solution c is given by c = b' v,
|
|
! where v is the solution of the system (r) (r)' v = g and
|
|
! g denotes the vector obtained from the old right hand side f, by
|
|
! removing the elements corresponding to a zero diagonal element of a.
|
|
! initialization.
|
|
do i=1,rank
|
|
do j=1,m
|
|
aa(i,j) = 0.
|
|
end do
|
|
end do
|
|
! form in aa the upper triangular matrix obtained from a by
|
|
! removing rows and columns with zero diagonal elements. form in ff
|
|
! the new right hand side by removing the elements of the old right
|
|
! hand side corresponding to a deleted row.
|
|
ii = 0
|
|
do i=1,n
|
|
if(a(i,1)<=tol) cycle
|
|
ii = ii+1
|
|
ff(ii) = f(i)
|
|
aa(ii,1) = a(i,1)
|
|
jj = ii
|
|
kk = 1
|
|
j = i
|
|
j1 = min0(j-1,m1)
|
|
do k=1,j1
|
|
j = j-1
|
|
if(a(j,1)<=tol) cycle
|
|
kk = kk+1
|
|
jj = jj-1
|
|
aa(jj,kk) = a(j,k+1)
|
|
end do
|
|
end do
|
|
! form successively in h the columns of a with a zero diagonal element.
|
|
ii = 0
|
|
do i=1,n
|
|
ii = ii+1
|
|
if(a(i,1)>tol) cycle
|
|
ii = ii-1
|
|
if(ii==0) cycle
|
|
jj = 1
|
|
j = i
|
|
j1 = min0(j-1,m1)
|
|
do k=1,j1
|
|
j = j-1
|
|
if(a(j,1)<=tol) cycle
|
|
h(jj) = a(j,k+1)
|
|
jj = jj+1
|
|
end do
|
|
do kk=jj,m
|
|
h(kk) = 0.
|
|
end do
|
|
! rotate this column into aa by givens transformations.
|
|
jj = ii
|
|
do i1=1,ii
|
|
j1 = min0(jj-1,m1)
|
|
piv = h(1)
|
|
if(piv==0.) then
|
|
if(j1==0) exit
|
|
do j2=1,j1
|
|
j3 = j2+1
|
|
h(j2) = h(j3)
|
|
end do
|
|
else
|
|
call fpgivs(piv,aa(jj,1),cs,sn)
|
|
if(j1==0) exit
|
|
kk = jj
|
|
do j2=1,j1
|
|
j3 = j2+1
|
|
kk = kk-1
|
|
call fprota(cs,sn,h(j3),aa(kk,j3))
|
|
h(j2) = h(j3)
|
|
end do
|
|
end if
|
|
jj = jj-1
|
|
h(j3) = 0.
|
|
end do
|
|
end do
|
|
! solve the system (aa) (f1) = ff
|
|
ff(rank) = ff(rank)/aa(rank,1)
|
|
i = rank-1
|
|
do j=2,rank
|
|
store = ff(i)
|
|
i1 = min0(j-1,m1)
|
|
k = i
|
|
do ii=1,i1
|
|
k = k+1
|
|
stor1 = ff(k)
|
|
stor2 = aa(i,ii+1)
|
|
store = store-stor1*stor2
|
|
end do
|
|
stor1 = aa(i,1)
|
|
ff(i) = store/stor1
|
|
i = i-1
|
|
end do
|
|
! solve the system (aa)' (f2) = f1
|
|
ff(1) = ff(1)/aa(1,1)
|
|
do j=2,rank
|
|
store = ff(j)
|
|
i1 = min0(j-1,m1)
|
|
k = j
|
|
do ii=1,i1
|
|
k = k-1
|
|
stor1 = ff(k)
|
|
stor2 = aa(k,ii+1)
|
|
store = store-stor1*stor2
|
|
end do
|
|
stor1 = aa(j,1)
|
|
ff(j) = store/stor1
|
|
end do
|
|
! premultiply f2 by the transpoze of a.
|
|
k = 0
|
|
do i=1,n
|
|
store = 0.
|
|
if(a(i,1)>tol) k = k+1
|
|
j1 = min0(i,m)
|
|
kk = k
|
|
ij = i+1
|
|
do j=1,j1
|
|
ij = ij-1
|
|
if(a(ij,1)<=tol) cycle
|
|
stor1 = a(ij,j)
|
|
stor2 = ff(kk)
|
|
store = store+stor1*stor2
|
|
kk = kk-1
|
|
end do
|
|
c(i) = store
|
|
end do
|
|
! add to the sum of squared residuals the contribution of putting
|
|
! to zero the small diagonal elements of matrix (a).
|
|
stor3 = 0.
|
|
do i=1,n
|
|
if(a(i,1)>tol) cycle
|
|
store = f(i)
|
|
i1 = min0(n-i,m1)
|
|
do j=1,i1
|
|
ij = i+j
|
|
stor1 = c(ij)
|
|
stor2 = a(i,j+1)
|
|
store = store-stor1*stor2
|
|
end do
|
|
fac = a(i,1)*c(i)
|
|
stor1 = a(i,1)
|
|
stor2 = c(i)
|
|
stor1 = stor1*stor2
|
|
stor3 = stor3+stor1*(stor1-store-store)
|
|
end do
|
|
fac = stor3
|
|
sq = sq+fac
|
|
end subroutine fprank
|
|
|
|
subroutine fporde(x,y,m,kx,ky,tx,nx,ty,ny,nummer,idx,nreg)
|
|
! subroutine fporde sorts the data points (x(i),y(i)),i=1,2,...,m
|
|
! according to the panel tx(l)<=x<tx(l+1),ty(k)<=y<ty(k+1), they belong
|
|
! to. for each panel a stack is constructed containing the numbers
|
|
! of data points lying inside; index(j),j=1,2,...,nreg points to the
|
|
! first data point in the jth panel while nummer(i),i=1,2,...,m gives
|
|
! the number of the next data point in the panel.
|
|
! ..
|
|
implicit none
|
|
! ..scalar arguments..
|
|
integer m,kx,ky,nx,ny,nreg
|
|
! ..array arguments..
|
|
real(wp_) x(m),y(m),tx(nx),ty(ny)
|
|
integer nummer(m),idx(nreg)
|
|
! ..local scalars..
|
|
real(wp_) xi,yi
|
|
integer i,im,k,kx1,ky1,k1,l,l1,nk1x,nk1y,num,nyy
|
|
! ..
|
|
kx1 = kx+1
|
|
ky1 = ky+1
|
|
nk1x = nx-kx1
|
|
nk1y = ny-ky1
|
|
nyy = nk1y-ky
|
|
do i=1,nreg
|
|
idx(i) = 0
|
|
end do
|
|
do im=1,m
|
|
xi = x(im)
|
|
yi = y(im)
|
|
l = kx1
|
|
l1 = l+1
|
|
do
|
|
if(xi<tx(l1) .or. l==nk1x) exit
|
|
l = l1
|
|
l1 = l+1
|
|
end do
|
|
k = ky1
|
|
k1 = k+1
|
|
do
|
|
if(yi<ty(k1) .or. k==nk1y) exit
|
|
k = k1
|
|
k1 = k+1
|
|
end do
|
|
num = (l-kx1)*nyy+k-ky
|
|
nummer(im) = idx(num)
|
|
idx(num) = im
|
|
end do
|
|
end subroutine fporde
|
|
|
|
subroutine fpgivs(piv,ww,cs,sn)
|
|
! subroutine fpgivs calculates the parameters of a givens
|
|
! transformation .
|
|
! ..
|
|
implicit none
|
|
! arguments
|
|
real(wp_), intent(in) :: piv
|
|
real(wp_), intent(out) :: cs, sn
|
|
real(wp_), intent(inout) :: ww
|
|
! local variables
|
|
real(wp_) :: dd, store
|
|
! parameters
|
|
real(wp_), parameter :: one = 0.1e+01_wp_
|
|
! ..
|
|
store = abs(piv)
|
|
if(store>=ww) dd = store*sqrt(one+(ww/piv)**2)
|
|
if(store<ww) dd = ww*sqrt(one+(piv/ww)**2)
|
|
cs = ww/dd
|
|
sn = piv/dd
|
|
ww = dd
|
|
end subroutine fpgivs
|
|
|
|
subroutine fpgrre(ifsx,ifsy,ifbx,ifby,x,mx,y,my,z,mz,kx,ky,tx,nx, &
|
|
ty,ny,p,c,nc,fp,fpx,fpy,mm,mynx,kx1,kx2,ky1,ky2,spx,spy,right,q, &
|
|
ax,ay,bx,by,nrx,nry)
|
|
! ..
|
|
implicit none
|
|
! arguments
|
|
integer, intent(in) :: kx, kx1, kx2, ky, ky1, ky2, mm, mx, my, mz, &
|
|
mynx, nx, ny, nc
|
|
integer, intent(inout) :: ifsx, ifsy, ifbx, ifby, nrx(mx), nry(my)
|
|
real(wp_), intent(in) :: p, x(mx), y(my), z(mz), tx(nx), ty(ny)
|
|
real(wp_), intent(out) :: fp, c(nc), q(mynx), ax(nx,kx2), ay(ny,ky2), fpx(nx), fpy(ny)
|
|
real(wp_), intent(inout) :: spx(mx,kx1), spy(my,ky1), bx(nx,kx2), by(ny,ky2)
|
|
real(wp_), intent(out) :: right(mm)
|
|
! local variables
|
|
real(wp_) :: arg, cs, fac, pinv, piv, sn, term
|
|
integer :: i, ibandx, ibandy, ic, iq, irot, it, iz,i1, i2, i3, &
|
|
j, k, k1, k2, l, l1, l2, ncof, nk1x, nk1y, nrold, nroldx, nroldy, &
|
|
number, numx, numx1, numy, numy1, n1
|
|
real(wp_), dimension(7) :: h
|
|
! parameters
|
|
real(wp_), parameter :: one = 1.0_wp_, half = 0.5_wp_
|
|
! ..subroutine references..
|
|
! fpback,fpbspl,fpgivs,fpdisc,fprota
|
|
! ..
|
|
! the b-spline coefficients of the smoothing spline are calculated as
|
|
! the least-squares solution of the over-determined linear system of
|
|
! equations (ay) c (ax)' = q where
|
|
!
|
|
! | (spx) | | (spy) |
|
|
! (ax) = | ---------- | (ay) = | ---------- |
|
|
! | (1/p) (bx) | | (1/p) (by) |
|
|
!
|
|
! | z ' 0 |
|
|
! q = | ------ |
|
|
! | 0 ' 0 |
|
|
!
|
|
! with c : the (ny-ky-1) x (nx-kx-1) matrix which contains the
|
|
! b-spline coefficients.
|
|
! z : the my x mx matrix which contains the function values.
|
|
! spx,spy: the mx x (nx-kx-1) and my x (ny-ky-1) observation
|
|
! matrices according to the least-squares problems in
|
|
! the x- and y-direction.
|
|
! bx,by : the (nx-2*kx-1) x (nx-kx-1) and (ny-2*ky-1) x (ny-ky-1)
|
|
! matrices which contain the discontinuity jumps of the
|
|
! derivatives of the b-splines in the x- and y-direction.
|
|
nk1x = nx-kx1
|
|
nk1y = ny-ky1
|
|
if(p>0.0_wp_) pinv = one/p
|
|
! it depends on the value of the flags ifsx,ifsy,ifbx and ifby and on
|
|
! the value of p whether the matrices (spx),(spy),(bx) and (by) still
|
|
! must be determined.
|
|
if(ifsx==0) then
|
|
! calculate the non-zero elements of the matrix (spx) which is the
|
|
! observation matrix according to the least-squares spline approximat-
|
|
! ion problem in the x-direction.
|
|
l = kx1
|
|
l1 = kx2
|
|
number = 0
|
|
do it=1,mx
|
|
arg = x(it)
|
|
do
|
|
if(arg<tx(l1) .or. l==nk1x) exit
|
|
l = l1
|
|
l1 = l+1
|
|
number = number+1
|
|
end do
|
|
call fpbspl(tx,nx,kx,arg,l,h)
|
|
do i=1,kx1
|
|
spx(it,i) = h(i)
|
|
end do
|
|
nrx(it) = number
|
|
end do
|
|
ifsx = 1
|
|
end if
|
|
if(ifsy==0) then
|
|
! calculate the non-zero elements of the matrix (spy) which is the
|
|
! observation matrix according to the least-squares spline approximat-
|
|
! ion problem in the y-direction.
|
|
l = ky1
|
|
l1 = ky2
|
|
number = 0
|
|
do it=1,my
|
|
arg = y(it)
|
|
do
|
|
if(arg<ty(l1) .or. l==nk1y) exit
|
|
l = l1
|
|
l1 = l+1
|
|
number = number+1
|
|
end do
|
|
call fpbspl(ty,ny,ky,arg,l,h)
|
|
do i=1,ky1
|
|
spy(it,i) = h(i)
|
|
end do
|
|
nry(it) = number
|
|
end do
|
|
ifsy = 1
|
|
end if
|
|
if(p>0.0_wp_) then
|
|
! calculate the non-zero elements of the matrix (bx).
|
|
if(ifbx==0 .and. nx/=2*kx1) then
|
|
call fpdisc(tx,nx,kx2,bx,nx)
|
|
ifbx = 1
|
|
end if
|
|
! calculate the non-zero elements of the matrix (by).
|
|
if(ifby==0 .and. ny/=2*ky1) then
|
|
call fpdisc(ty,ny,ky2,by,ny)
|
|
ifby = 1
|
|
end if
|
|
end if
|
|
! reduce the matrix (ax) to upper triangular form (rx) using givens
|
|
! rotations. apply the same transformations to the rows of matrix q
|
|
! to obtain the my x (nx-kx-1) matrix g.
|
|
! store matrix (rx) into (ax) and g into q.
|
|
l = my*nk1x
|
|
! initialization.
|
|
do i=1,l
|
|
q(i) = 0.0_wp_
|
|
end do
|
|
do i=1,nk1x
|
|
do j=1,kx2
|
|
ax(i,j) = 0.0_wp_
|
|
end do
|
|
end do
|
|
l = 0
|
|
nrold = 0
|
|
! ibandx denotes the bandwidth of the matrices (ax) and (rx).
|
|
ibandx = kx1
|
|
do it=1,mx
|
|
number = nrx(it)
|
|
do
|
|
if(nrold/=number) then
|
|
if(p<=0.0_wp_) then
|
|
nrold = nrold+1
|
|
cycle
|
|
end if
|
|
ibandx = kx2
|
|
! fetch a new row of matrix (bx).
|
|
n1 = nrold+1
|
|
do j=1,kx2
|
|
h(j) = bx(n1,j)*pinv
|
|
end do
|
|
! find the appropriate column of q.
|
|
do j=1,my
|
|
right(j) = 0.0_wp_
|
|
end do
|
|
irot = nrold
|
|
else
|
|
! fetch a new row of matrix (spx).
|
|
h(ibandx) = 0.0_wp_
|
|
do j=1,kx1
|
|
h(j) = spx(it,j)
|
|
end do
|
|
! find the appropriate column of q.
|
|
do j=1,my
|
|
l = l+1
|
|
right(j) = z(l)
|
|
end do
|
|
irot = number
|
|
end if
|
|
! rotate the new row of matrix (ax) into triangle.
|
|
do i=1,ibandx
|
|
irot = irot+1
|
|
piv = h(i)
|
|
if(piv==0.0_wp_) cycle
|
|
! calculate the parameters of the givens transformation.
|
|
call fpgivs(piv,ax(irot,1),cs,sn)
|
|
! apply that transformation to the rows of matrix q.
|
|
iq = (irot-1)*my
|
|
do j=1,my
|
|
iq = iq+1
|
|
call fprota(cs,sn,right(j),q(iq))
|
|
end do
|
|
! apply that transformation to the columns of (ax).
|
|
if(i==ibandx) exit
|
|
i2 = 1
|
|
i3 = i+1
|
|
do j=i3,ibandx
|
|
i2 = i2+1
|
|
call fprota(cs,sn,h(j),ax(irot,i2))
|
|
end do
|
|
end do
|
|
if(nrold==number) exit
|
|
nrold = nrold+1
|
|
end do
|
|
end do
|
|
! reduce the matrix (ay) to upper triangular form (ry) using givens
|
|
! rotations. apply the same transformations to the columns of matrix g
|
|
! to obtain the (ny-ky-1) x (nx-kx-1) matrix h.
|
|
! store matrix (ry) into (ay) and h into c.
|
|
ncof = nk1x*nk1y
|
|
! initialization.
|
|
do i=1,ncof
|
|
c(i) = 0.0_wp_
|
|
end do
|
|
do i=1,nk1y
|
|
do j=1,ky2
|
|
ay(i,j) = 0.0_wp_
|
|
end do
|
|
end do
|
|
nrold = 0
|
|
! ibandy denotes the bandwidth of the matrices (ay) and (ry).
|
|
ibandy = ky1
|
|
do it=1,my
|
|
number = nry(it)
|
|
do
|
|
if(nrold/=number) then
|
|
if(p<=0.0_wp_) then
|
|
nrold = nrold+1
|
|
cycle
|
|
end if
|
|
ibandy = ky2
|
|
! fetch a new row of matrix (by).
|
|
n1 = nrold+1
|
|
do j=1,ky2
|
|
h(j) = by(n1,j)*pinv
|
|
end do
|
|
! find the appropiate row of g.
|
|
do j=1,nk1x
|
|
right(j) = 0.0_wp_
|
|
end do
|
|
irot = nrold
|
|
! fetch a new row of matrix (spy)
|
|
else
|
|
h(ibandy) = 0.0_wp_
|
|
do j=1,ky1
|
|
h(j) = spy(it,j)
|
|
end do
|
|
! find the appropiate row of g.
|
|
l = it
|
|
do j=1,nk1x
|
|
right(j) = q(l)
|
|
l = l+my
|
|
end do
|
|
irot = number
|
|
end if
|
|
! rotate the new row of matrix (ay) into triangle.
|
|
do i=1,ibandy
|
|
irot = irot+1
|
|
piv = h(i)
|
|
if(piv==0.0_wp_) cycle
|
|
! calculate the parameters of the givens transformation.
|
|
call fpgivs(piv,ay(irot,1),cs,sn)
|
|
! apply that transformation to the colums of matrix g.
|
|
ic = irot
|
|
do j=1,nk1x
|
|
call fprota(cs,sn,right(j),c(ic))
|
|
ic = ic+nk1y
|
|
end do
|
|
! apply that transformation to the columns of matrix (ay).
|
|
if(i==ibandy) exit
|
|
i2 = 1
|
|
i3 = i+1
|
|
do j=i3,ibandy
|
|
i2 = i2+1
|
|
call fprota(cs,sn,h(j),ay(irot,i2))
|
|
end do
|
|
end do
|
|
if(nrold==number) exit
|
|
nrold = nrold+1
|
|
end do
|
|
end do
|
|
! backward substitution to obtain the b-spline coefficients as the
|
|
! solution of the linear system (ry) c (rx)' = h.
|
|
! first step: solve the system (ry) (c1) = h.
|
|
k = 1
|
|
do i=1,nk1x
|
|
call fpback(ay,c(k),nk1y,ibandy,c(k),ny)
|
|
k = k+nk1y
|
|
end do
|
|
! second step: solve the system c (rx)' = (c1).
|
|
k = 0
|
|
do j=1,nk1y
|
|
k = k+1
|
|
l = k
|
|
do i=1,nk1x
|
|
right(i) = c(l)
|
|
l = l+nk1y
|
|
end do
|
|
call fpback(ax,right,nk1x,ibandx,right,nx)
|
|
l = k
|
|
do i=1,nk1x
|
|
c(l) = right(i)
|
|
l = l+nk1y
|
|
end do
|
|
end do
|
|
! calculate the quantities
|
|
! res(i,j) = (z(i,j) - s(x(i),y(j)))**2 , i=1,2,..,mx;j=1,2,..,my
|
|
! fp = sumi=1,mx(sumj=1,my(res(i,j)))
|
|
! fpx(r) = sum''i(sumj=1,my(res(i,j))) , r=1,2,...,nx-2*kx-1
|
|
! tx(r+kx) <= x(i) <= tx(r+kx+1)
|
|
! fpy(r) = sumi=1,mx(sum''j(res(i,j))) , r=1,2,...,ny-2*ky-1
|
|
! ty(r+ky) <= y(j) <= ty(r+ky+1)
|
|
fp = 0.0_wp_
|
|
do i=1,nx
|
|
fpx(i) = 0.0_wp_
|
|
end do
|
|
do i=1,ny
|
|
fpy(i) = 0.0_wp_
|
|
end do
|
|
nk1y = ny-ky1
|
|
iz = 0
|
|
nroldx = 0
|
|
! main loop for the different grid points.
|
|
do i1=1,mx
|
|
numx = nrx(i1)
|
|
numx1 = numx+1
|
|
nroldy = 0
|
|
do i2=1,my
|
|
numy = nry(i2)
|
|
numy1 = numy+1
|
|
iz = iz+1
|
|
! evaluate s(x,y) at the current grid point by making the sum of the
|
|
! cross products of the non-zero b-splines at (x,y), multiplied with
|
|
! the appropiate b-spline coefficients.
|
|
term = 0.0_wp_
|
|
k1 = numx*nk1y+numy
|
|
do l1=1,kx1
|
|
k2 = k1
|
|
fac = spx(i1,l1)
|
|
do l2=1,ky1
|
|
k2 = k2+1
|
|
term = term+fac*spy(i2,l2)*c(k2)
|
|
end do
|
|
k1 = k1+nk1y
|
|
end do
|
|
! calculate the squared residual at the current grid point.
|
|
term = (z(iz)-term)**2
|
|
! adjust the different parameters.
|
|
fp = fp+term
|
|
fpx(numx1) = fpx(numx1)+term
|
|
fpy(numy1) = fpy(numy1)+term
|
|
fac = term*half
|
|
if(numy/=nroldy) then
|
|
fpy(numy1) = fpy(numy1)-fac
|
|
fpy(numy) = fpy(numy)+fac
|
|
end if
|
|
nroldy = numy
|
|
if(numx/=nroldx) then
|
|
fpx(numx1) = fpx(numx1)-fac
|
|
fpx(numx) = fpx(numx)+fac
|
|
end if
|
|
end do
|
|
nroldx = numx
|
|
end do
|
|
end subroutine fpgrre
|
|
|
|
subroutine fpknot(x,m,t,n,fpint,nrdata,nrint,nest,istart)
|
|
! subroutine fpknot locates an additional knot for a spline of degree
|
|
! k and adjusts the corresponding parameters,i.e.
|
|
! t : the position of the knots.
|
|
! n : the number of knots.
|
|
! nrint : the number of knotintervals.
|
|
! fpint : the sum of squares of residual right hand sides
|
|
! for each knot interval.
|
|
! nrdata: the number of data points inside each knot interval.
|
|
! istart indicates that the smallest data point at which the new knot
|
|
! may be added is x(istart+1)
|
|
! ..
|
|
implicit none
|
|
! arguments
|
|
integer, intent(in) :: m, nest, istart
|
|
integer, intent(inout) :: n, nrint, nrdata(nest)
|
|
real(wp_), intent(in) :: x(m)
|
|
real(wp_), intent(inout) :: t(nest), fpint(nest)
|
|
! local variables
|
|
real(wp_) :: an, am, fpmax
|
|
integer :: ihalf, j, jbegin, jj, jk, jpoint, k, maxbeg, maxpt, next, nrx, number
|
|
! ..
|
|
k = (n-nrint-1)/2
|
|
! search for knot interval t(number+k) <= x <= t(number+k+1) where
|
|
! fpint(number) is maximal on the condition that nrdata(number)
|
|
! not equals zero.
|
|
fpmax = 0.0_wp_
|
|
jbegin = istart
|
|
do j=1,nrint
|
|
jpoint = nrdata(j)
|
|
if(fpmax<fpint(j) .and. jpoint/=0) then
|
|
fpmax = fpint(j)
|
|
number = j
|
|
maxpt = jpoint
|
|
maxbeg = jbegin
|
|
end if
|
|
jbegin = jbegin+jpoint+1
|
|
end do
|
|
! let coincide the new knot t(number+k+1) with a data point x(nrx)
|
|
! inside the old knot interval t(number+k) <= x <= t(number+k+1).
|
|
ihalf = maxpt/2+1
|
|
nrx = maxbeg+ihalf
|
|
next = number+1
|
|
if(next<=nrint) then
|
|
! adjust the different parameters.
|
|
do j=next,nrint
|
|
jj = next+nrint-j
|
|
fpint(jj+1) = fpint(jj)
|
|
nrdata(jj+1) = nrdata(jj)
|
|
jk = jj+k
|
|
t(jk+1) = t(jk)
|
|
end do
|
|
end if
|
|
nrdata(number) = ihalf-1
|
|
nrdata(next) = maxpt-ihalf
|
|
am = maxpt
|
|
an = nrdata(number)
|
|
fpint(number) = fpmax*an/am
|
|
an = nrdata(next)
|
|
fpint(next) = fpmax*an/am
|
|
jk = next+k
|
|
t(jk) = x(nrx)
|
|
n = n+1
|
|
nrint = nrint+1
|
|
end subroutine fpknot
|
|
|
|
subroutine fpregr(iopt,x,mx,y,my,z,mz,xb,xe,yb,ye,kx,ky,s, &
|
|
nxest,nyest,tol,maxit,nc,nx,tx,ny,ty,c,fp,fp0,fpold,reducx, &
|
|
reducy,fpintx,fpinty,lastdi,nplusx,nplusy,nrx,nry,nrdatx,nrdaty, &
|
|
wrk,lwrk,ier)
|
|
! ..
|
|
implicit none
|
|
! arguments
|
|
integer, intent(in) :: kx,ky,mx,my,mz,nxest,nyest,maxit,nc,lwrk,iopt
|
|
integer, intent(inout) :: nx,ny,lastdi,nplusx,nplusy, &
|
|
nrdatx(nxest),nrdaty(nyest),nrx(mx),nry(my)
|
|
integer, intent(out) :: ier
|
|
real(wp_), intent(in) :: xb,xe,yb,ye,s,tol,x(mx),y(my),z(mz)
|
|
real(wp_), intent(out) :: fp,c(nc)
|
|
real(wp_), intent(inout) :: fp0,fpold,reducx,reducy,tx(nxest), &
|
|
ty(nyest),fpintx(nxest),fpinty(nyest),wrk(lwrk)
|
|
! local variables
|
|
real(wp_) :: acc,fpms,f1,f2,f3,p,p1,p2,p3,rn
|
|
integer :: i,ich1,ich3,ifbx,ifby,ifsx,ifsy,iter,j,kx1,kx2,ky1,ky2,k3, &
|
|
l,lax,lay,lbx,lby,lq,lri,lsx,lsy,mk1,mm,mpm,mynx,nmaxx,nmaxy,nminx,nminy, &
|
|
nplx,nply,npl1,nrintx,nrinty,nxe,nxk,nye
|
|
! parameters
|
|
real(wp_), parameter :: one=1.0_wp_,half=0.5_wp_,con1=0.1_wp_,con9=0.9_wp_,con4=0.4e-01_wp_
|
|
|
|
! ..subroutine references..
|
|
! fpgrre,fpknot
|
|
! ..
|
|
! we partition the working space.
|
|
kx1 = kx+1
|
|
ky1 = ky+1
|
|
kx2 = kx1+1
|
|
ky2 = ky1+1
|
|
lsx = 1
|
|
lsy = lsx+mx*kx1
|
|
lri = lsy+my*ky1
|
|
mm = max(nxest,my)
|
|
lq = lri+mm
|
|
mynx = nxest*my
|
|
lax = lq+mynx
|
|
nxk = nxest*kx2
|
|
lbx = lax+nxk
|
|
lay = lbx+nxk
|
|
lby = lay+nyest*ky2
|
|
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
! part 1: determination of the number of knots and their position. c
|
|
! **************************************************************** c
|
|
! given a set of knots we compute the least-squares spline sinf(x,y), c
|
|
! and the corresponding sum of squared residuals fp=f(p=inf). c
|
|
! if iopt=-1 sinf(x,y) is the requested approximation. c
|
|
! if iopt=0 or iopt=1 we check whether we can accept the knots: c
|
|
! if fp <=s we will continue with the current set of knots. c
|
|
! if fp > s we will increase the number of knots and compute the c
|
|
! corresponding least-squares spline until finally fp<=s. c
|
|
! the initial choice of knots depends on the value of s and iopt. c
|
|
! if s=0 we have spline interpolation; in that case the number of c
|
|
! knots equals nmaxx = mx+kx+1 and nmaxy = my+ky+1. c
|
|
! if s>0 and c
|
|
! *iopt=0 we first compute the least-squares polynomial of degree c
|
|
! kx in x and ky in y; nx=nminx=2*kx+2 and ny=nymin=2*ky+2. c
|
|
! *iopt=1 we start with the knots found at the last call of the c
|
|
! routine, except for the case that s > fp0; then we can compute c
|
|
! the least-squares polynomial directly. c
|
|
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
! determine the number of knots for polynomial approximation.
|
|
nminx = 2*kx1
|
|
nminy = 2*ky1
|
|
if(iopt>=0) then
|
|
! acc denotes the absolute tolerance for the root of f(p)=s.
|
|
acc = tol*s
|
|
! find nmaxx and nmaxy which denote the number of knots in x- and y-
|
|
! direction in case of spline interpolation.
|
|
nmaxx = mx+kx1
|
|
nmaxy = my+ky1
|
|
! find nxe and nye which denote the maximum number of knots
|
|
! allowed in each direction
|
|
nxe = min(nmaxx,nxest)
|
|
nye = min(nmaxy,nyest)
|
|
if(s<=0.0_wp_) then
|
|
! if s = 0, s(x,y) is an interpolating spline.
|
|
nx = nmaxx
|
|
ny = nmaxy
|
|
! test whether the required storage space exceeds the available one.
|
|
if(ny>nyest .or. nx>nxest) then
|
|
ier = 1
|
|
return
|
|
end if
|
|
! find the position of the interior knots in case of interpolation.
|
|
! the knots in the x-direction.
|
|
mk1 = mx-kx1
|
|
if(mk1/=0) then
|
|
k3 = kx/2
|
|
i = kx1+1
|
|
j = k3+2
|
|
if(k3*2/=kx) then
|
|
do l=1,mk1
|
|
tx(i) = x(j)
|
|
i = i+1
|
|
j = j+1
|
|
end do
|
|
else
|
|
do l=1,mk1
|
|
tx(i) = (x(j)+x(j-1))*half
|
|
i = i+1
|
|
j = j+1
|
|
end do
|
|
end if
|
|
end if
|
|
! the knots in the y-direction.
|
|
mk1 = my-ky1
|
|
if(mk1/=0) then
|
|
k3 = ky/2
|
|
i = ky1+1
|
|
j = k3+2
|
|
if(k3*2/=ky) then
|
|
do l=1,mk1
|
|
ty(i) = y(j)
|
|
i = i+1
|
|
j = j+1
|
|
end do
|
|
else
|
|
do l=1,mk1
|
|
ty(i) = (y(j)+y(j-1))*half
|
|
i = i+1
|
|
j = j+1
|
|
end do
|
|
end if
|
|
end if
|
|
else
|
|
! if s > 0 our initial choice of knots depends on the value of iopt.
|
|
if(iopt/=0 .and. fp0>s) then
|
|
! if iopt=1 and fp0 > s we start computing the least- squares spline
|
|
! according to the set of knots found at the last call of the routine.
|
|
! we determine the number of grid coordinates x(i) inside each knot
|
|
! interval (tx(l),tx(l+1)).
|
|
l = kx2
|
|
j = 1
|
|
nrdatx(1) = 0
|
|
mpm = mx-1
|
|
do i=2,mpm
|
|
nrdatx(j) = nrdatx(j)+1
|
|
if(x(i)<tx(l)) cycle
|
|
nrdatx(j) = nrdatx(j)-1
|
|
l = l+1
|
|
j = j+1
|
|
nrdatx(j) = 0
|
|
end do
|
|
! we determine the number of grid coordinates y(i) inside each knot
|
|
! interval (ty(l),ty(l+1)).
|
|
l = ky2
|
|
j = 1
|
|
nrdaty(1) = 0
|
|
mpm = my-1
|
|
do i=2,mpm
|
|
nrdaty(j) = nrdaty(j)+1
|
|
if(y(i)<ty(l)) cycle
|
|
nrdaty(j) = nrdaty(j)-1
|
|
l = l+1
|
|
j = j+1
|
|
nrdaty(j) = 0
|
|
end do
|
|
else
|
|
! if iopt=0 or iopt=1 and s>=fp0, we start computing the least-squares
|
|
! polynomial of degree kx in x and ky in y (which is a spline without
|
|
! interior knots).
|
|
nx = nminx
|
|
ny = nminy
|
|
nrdatx(1) = mx-2
|
|
nrdaty(1) = my-2
|
|
lastdi = 0
|
|
nplusx = 0
|
|
nplusy = 0
|
|
fp0 = 0.0_wp_
|
|
fpold = 0.0_wp_
|
|
reducx = 0.0_wp_
|
|
reducy = 0.0_wp_
|
|
end if
|
|
end if
|
|
end if
|
|
mpm = mx+my
|
|
ifsx = 0
|
|
ifsy = 0
|
|
ifbx = 0
|
|
ifby = 0
|
|
p = -one
|
|
! main loop for the different sets of knots.mpm=mx+my is a save upper
|
|
! bound for the number of trials.
|
|
do iter=1,mpm
|
|
if(nx==nminx .and. ny==nminy) ier = -2
|
|
! find nrintx (nrinty) which is the number of knot intervals in the
|
|
! x-direction (y-direction).
|
|
nrintx = nx-nminx+1
|
|
nrinty = ny-nminy+1
|
|
! find ncof, the number of b-spline coefficients for the current set
|
|
! of knots.
|
|
! nk1x = nx-kx1
|
|
! nk1y = ny-ky1
|
|
! ncof = nk1x*nk1y
|
|
! find the position of the additional knots which are needed for the
|
|
! b-spline representation of s(x,y).
|
|
i = nx
|
|
do j=1,kx1
|
|
tx(j) = xb
|
|
tx(i) = xe
|
|
i = i-1
|
|
end do
|
|
i = ny
|
|
do j=1,ky1
|
|
ty(j) = yb
|
|
ty(i) = ye
|
|
i = i-1
|
|
end do
|
|
! find the least-squares spline sinf(x,y) and calculate for each knot
|
|
! interval tx(j+kx)<=x<=tx(j+kx+1) (ty(j+ky)<=y<=ty(j+ky+1)) the sum
|
|
! of squared residuals fpintx(j),j=1,2,...,nx-2*kx-1 (fpinty(j),j=1,2,
|
|
! ...,ny-2*ky-1) for the data points having their absciss (ordinate)-
|
|
! value belonging to that interval.
|
|
! fp gives the total sum of squared residuals.
|
|
call fpgrre(ifsx,ifsy,ifbx,ifby,x,mx,y,my,z,mz,kx,ky,tx,nx,ty, &
|
|
ny,p,c,nc,fp,fpintx,fpinty,mm,mynx,kx1,kx2,ky1,ky2,wrk(lsx), &
|
|
wrk(lsy),wrk(lri),wrk(lq),wrk(lax),wrk(lay),wrk(lbx),wrk(lby), &
|
|
nrx,nry)
|
|
if(ier==(-2)) fp0 = fp
|
|
! test whether the least-squares spline is an acceptable solution.
|
|
if(iopt<0) return
|
|
fpms = fp-s
|
|
if(abs(fpms) < acc) return
|
|
! if f(p=inf) < s, we accept the choice of knots.
|
|
if(fpms<0.0_wp_) exit
|
|
! if nx=nmaxx and ny=nmaxy, sinf(x,y) is an interpolating spline.
|
|
if(nx==nmaxx .and. ny==nmaxy) then
|
|
ier = -1
|
|
fp = 0.0_wp_
|
|
return
|
|
end if
|
|
! increase the number of knots.
|
|
! if nx=nxe and ny=nye we cannot further increase the number of knots
|
|
! because of the storage capacity limitation.
|
|
if(nx==nxe .and. ny==nye) then
|
|
ier = 1
|
|
return
|
|
end if
|
|
ier = 0
|
|
! adjust the parameter reducx or reducy according to the direction
|
|
! in which the last added knots were located.
|
|
if(lastdi<0) then
|
|
reducx = fpold-fp
|
|
else if(lastdi>0) then
|
|
reducy = fpold-fp
|
|
end if
|
|
! store the sum of squared residuals for the current set of knots.
|
|
fpold = fp
|
|
! find nplx, the number of knots we should add in the x-direction.
|
|
nplx = 1
|
|
if(nx/=nminx) then
|
|
npl1 = nplusx*2
|
|
rn = nplusx
|
|
if(reducx>acc) npl1 = int(rn*fpms/reducx)
|
|
nplx = min(nplusx*2,max(npl1,nplusx/2,1))
|
|
end if
|
|
! find nply, the number of knots we should add in the y-direction.
|
|
nply = 1
|
|
if(ny/=nminy) then
|
|
npl1 = nplusy*2
|
|
rn = nplusy
|
|
if(reducy>acc) npl1 = int(rn*fpms/reducy)
|
|
nply = min0(nplusy*2,max0(npl1,nplusy/2,1))
|
|
end if
|
|
if (ny==nye .or. (nx/=nxe .and. ny/=nye .and.&
|
|
(nplx<nply .or. (nplx==nply .and. lastdi>=0)))) then
|
|
! addition in the x-direction.
|
|
lastdi = -1
|
|
nplusx = nplx
|
|
ifsx = 0
|
|
do l=1,nplusx
|
|
! add a new knot in the x-direction
|
|
call fpknot(x,mx,tx,nx,fpintx,nrdatx,nrintx,nxest,1)
|
|
! test whether we cannot further increase the number of knots in the
|
|
! x-direction.
|
|
if (nx==nxe) exit
|
|
end do
|
|
else
|
|
! addition in the y-direction.
|
|
lastdi = 1
|
|
nplusy = nply
|
|
ifsy = 0
|
|
do l=1,nplusy
|
|
! add a new knot in the y-direction.
|
|
call fpknot(y,my,ty,ny,fpinty,nrdaty,nrinty,nyest,1)
|
|
! test whether we cannot further increase the number of knots in the
|
|
! y-direction.
|
|
if (ny==nye) exit
|
|
end do
|
|
end if
|
|
! restart the computations with the new set of knots.
|
|
end do
|
|
! test whether the least-squares polynomial is a solution of our
|
|
! approximation problem.
|
|
if(ier==(-2)) return
|
|
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
! part 2: determination of the smoothing spline sp(x,y) c
|
|
! ***************************************************** c
|
|
! we have determined the number of knots and their position. we now c
|
|
! compute the b-spline coefficients of the smoothing spline sp(x,y). c
|
|
! this smoothing spline varies with the parameter p in such a way thatc
|
|
! f(p) = sumi=1,mx(sumj=1,my((z(i,j)-sp(x(i),y(j)))**2) c
|
|
! is a continuous, strictly decreasing function of p. moreover the c
|
|
! least-squares polynomial corresponds to p=0 and the least-squares c
|
|
! spline to p=infinity. iteratively we then have to determine the c
|
|
! positive value of p such that f(p)=s. the process which is proposed c
|
|
! here makes use of rational interpolation. f(p) is approximated by a c
|
|
! rational function r(p)=(u*p+v)/(p+w); three values of p (p1,p2,p3) c
|
|
! with corresponding values of f(p) (f1=f(p1)-s,f2=f(p2)-s,f3=f(p3)-s)c
|
|
! are used to calculate the new value of p such that r(p)=s. c
|
|
! convergence is guaranteed by taking f1 > 0 and f3 < 0. c
|
|
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
! initial value for p.
|
|
p1 = 0.0_wp_
|
|
f1 = fp0-s
|
|
p3 = -one
|
|
f3 = fpms
|
|
p = one
|
|
ich1 = 0
|
|
ich3 = 0
|
|
! iteration process to find the root of f(p)=s.
|
|
do iter = 1,maxit
|
|
! find the smoothing spline sp(x,y) and the corresponding sum of
|
|
! squared residuals fp.
|
|
call fpgrre(ifsx,ifsy,ifbx,ifby,x,mx,y,my,z,mz,kx,ky,tx,nx,ty, &
|
|
ny,p,c,nc,fp,fpintx,fpinty,mm,mynx,kx1,kx2,ky1,ky2,wrk(lsx), &
|
|
wrk(lsy),wrk(lri),wrk(lq),wrk(lax),wrk(lay),wrk(lbx),wrk(lby), &
|
|
nrx,nry)
|
|
! test whether the approximation sp(x,y) is an acceptable solution.
|
|
fpms = fp-s
|
|
if(abs(fpms)<acc) return
|
|
! test whether the maximum allowable number of iterations has been
|
|
! reached.
|
|
if(iter==maxit) then
|
|
ier = 3
|
|
return
|
|
end if
|
|
! carry out one more step of the iteration process.
|
|
p2 = p
|
|
f2 = fpms
|
|
if(ich3==0) then
|
|
if((f2-f3)<=acc) then
|
|
! our initial choice of p is too large.
|
|
p3 = p2
|
|
f3 = f2
|
|
p = p*con4
|
|
if(p<=p1) p = p1*con9 + p2*con1
|
|
cycle
|
|
end if
|
|
if(f2<0.0_wp_) ich3 = 1
|
|
end if
|
|
if(ich1==0) then
|
|
if((f1-f2)<=acc) then
|
|
! our initial choice of p is too small
|
|
p1 = p2
|
|
f1 = f2
|
|
p = p/con4
|
|
if(p3>=0.0_wp_ .and. p>=p3) p = p2*con1 + p3*con9
|
|
cycle
|
|
end if
|
|
! test whether the iteration process proceeds as theoretically
|
|
! expected.
|
|
if(f2>0.0_wp_) ich1 = 1
|
|
end if
|
|
if(f2>=f1 .or. f2<=f3) then
|
|
ier = 2
|
|
return
|
|
end if
|
|
! find the new value of p.
|
|
p = fprati(p1,f1,p2,f2,p3,f3)
|
|
end do
|
|
end subroutine fpregr
|
|
|
|
subroutine fprota(cs,sn,a,b)
|
|
! subroutine fprota applies a givens rotation to a and b.
|
|
! ..
|
|
implicit none
|
|
! arguments
|
|
real(wp_), intent(in) :: cs, sn
|
|
real(wp_), intent(inout) :: a, b
|
|
! local variables
|
|
real(wp_) :: stor1,stor2
|
|
! ..
|
|
stor1 = a
|
|
stor2 = b
|
|
b = cs*stor2+sn*stor1
|
|
a = cs*stor1-sn*stor2
|
|
end subroutine fprota
|
|
|
|
function fprati(p1,f1,p2,f2,p3,f3)
|
|
! given three points (p1,f1),(p2,f2) and (p3,f3), function fprati
|
|
! gives the value of p such that the rational interpolating function
|
|
! of the form r(p) = (u*p+v)/(p+w) equals zero at p.
|
|
! ..
|
|
implicit none
|
|
real(wp_) :: fprati
|
|
! arguments
|
|
real(wp_), intent(in) :: p2, f2
|
|
real(wp_), intent(inout) :: p1, f1, p3, f3
|
|
! local variables
|
|
real(wp_) :: h1, h2, h3, p
|
|
! ..
|
|
if(p3<=0.0_wp_) then
|
|
! value of p in case p3 = infinity.
|
|
p = (p1*(f1-f3)*f2-p2*(f2-f3)*f1)/((f1-f2)*f3)
|
|
else
|
|
! value of p in case p3 ^= infinity.
|
|
h1 = f1*(f2-f3)
|
|
h2 = f2*(f3-f1)
|
|
h3 = f3*(f1-f2)
|
|
p = -(p1*p2*h3+p2*p3*h1+p3*p1*h2)/(p1*h1+p2*h2+p3*h3)
|
|
! adjust the value of p1,f1,p3 and f3 such that f1 > 0 and f3 < 0.
|
|
end if
|
|
if(f2>=0.0_wp_) then
|
|
p1 = p2
|
|
f1 = f2
|
|
else
|
|
p3 = p2
|
|
f3 = f2
|
|
end if
|
|
fprati = p
|
|
end function fprati
|
|
|
|
subroutine regrid(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)
|
|
! given the set of values z(i,j) on the rectangular grid (x(i),y(j)),
|
|
! i=1,...,mx;j=1,...,my, subroutine regrid determines a smooth bivar-
|
|
! iate spline approximation s(x,y) of degrees kx and ky on the rect-
|
|
! angle xb <= x <= xe, yb <= y <= ye.
|
|
! if iopt = -1 regrid calculates the least-squares spline according
|
|
! to a given set of knots.
|
|
! if iopt >= 0 the total numbers nx and ny of these knots and their
|
|
! position tx(j),j=1,...,nx and ty(j),j=1,...,ny are chosen automatic-
|
|
! ally by the routine. the smoothness of s(x,y) is then achieved by
|
|
! minimalizing the discontinuity jumps in the derivatives of s(x,y)
|
|
! across the boundaries of the subpanels (tx(i),tx(i+1))*(ty(j),ty(j+1).
|
|
! the amounth of smoothness is determined by the condition that f(p) =
|
|
! sum ((z(i,j)-s(x(i),y(j))))**2) be <= s, with s a given non-negative
|
|
! constant, called the smoothing factor.
|
|
! the fit is given in the b-spline representation (b-spline coefficients
|
|
! c((ny-ky-1)*(i-1)+j),i=1,...,nx-kx-1;j=1,...,ny-ky-1) and can be eval-
|
|
! uated by means of subroutine bispev.
|
|
!
|
|
! calling sequence:
|
|
! call regrid(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)
|
|
!
|
|
! parameters:
|
|
! iopt : integer flag. on entry iopt must specify whether a least-
|
|
! squares spline (iopt=-1) or a smoothing spline (iopt=0 or 1)
|
|
! must be determined.
|
|
! if iopt=0 the routine will start with an initial set of knots
|
|
! tx(i)=xb,tx(i+kx+1)=xe,i=1,...,kx+1;ty(i)=yb,ty(i+ky+1)=ye,i=
|
|
! 1,...,ky+1. if iopt=1 the routine will continue with the set
|
|
! of knots found at the last call of the routine.
|
|
! attention: a call with iopt=1 must always be immediately pre-
|
|
! ceded by another call with iopt=1 or iopt=0 and
|
|
! s/=0.
|
|
! unchanged on exit.
|
|
! mx : integer. on entry mx must specify the number of grid points
|
|
! along the x-axis. mx > kx . unchanged on exit.
|
|
! x : real array of dimension at least (mx). before entry, x(i)
|
|
! must be set to the x-co-ordinate of the i-th grid point
|
|
! along the x-axis, for i=1,2,...,mx. these values must be
|
|
! supplied in strictly ascending order. unchanged on exit.
|
|
! my : integer. on entry my must specify the number of grid points
|
|
! along the y-axis. my > ky . unchanged on exit.
|
|
! y : real array of dimension at least (my). before entry, y(j)
|
|
! must be set to the y-co-ordinate of the j-th grid point
|
|
! along the y-axis, for j=1,2,...,my. these values must be
|
|
! supplied in strictly ascending order. unchanged on exit.
|
|
! z : real array of dimension at least (mx*my).
|
|
! before entry, z(my*(i-1)+j) must be set to the data value at
|
|
! the grid point (x(i),y(j)) for i=1,...,mx and j=1,...,my.
|
|
! unchanged on exit.
|
|
! xb,xe : real values. on entry xb,xe,yb and ye must specify the bound-
|
|
! yb,ye aries of the rectangular approximation domain.
|
|
! xb<=x(i)<=xe,i=1,...,mx; yb<=y(j)<=ye,j=1,...,my.
|
|
! unchanged on exit.
|
|
! kx,ky : integer values. on entry kx and ky must specify the degrees
|
|
! of the spline. 1<=kx,ky<=5. it is recommended to use bicubic
|
|
! (kx=ky=3) splines. unchanged on exit.
|
|
! s : real. on entry (in case iopt>=0) s must specify the smoothing
|
|
! factor. s >=0. unchanged on exit.
|
|
! for advice on the choice of s see further comments
|
|
! nxest : integer. unchanged on exit.
|
|
! nyest : integer. unchanged on exit.
|
|
! on entry, nxest and nyest must specify an upper bound for the
|
|
! number of knots required in the x- and y-directions respect.
|
|
! these numbers will also determine the storage space needed by
|
|
! the routine. nxest >= 2*(kx+1), nyest >= 2*(ky+1).
|
|
! in most practical situation nxest = mx/2, nyest=my/2, will
|
|
! be sufficient. always large enough are nxest=mx+kx+1, nyest=
|
|
! my+ky+1, the number of knots needed for interpolation (s=0).
|
|
! see also further comments.
|
|
! nx : integer.
|
|
! unless ier=10 (in case iopt >=0), nx will contain the total
|
|
! number of knots with respect to the x-variable, of the spline
|
|
! approximation returned. if the computation mode iopt=1 is
|
|
! used, the value of nx should be left unchanged between sub-
|
|
! sequent calls.
|
|
! in case iopt=-1, the value of nx should be specified on entry
|
|
! tx : real array of dimension nmax.
|
|
! on succesful exit, this array will contain the knots of the
|
|
! spline with respect to the x-variable, i.e. the position of
|
|
! the interior knots tx(kx+2),...,tx(nx-kx-1) as well as the
|
|
! position of the additional knots tx(1)=...=tx(kx+1)=xb and
|
|
! tx(nx-kx)=...=tx(nx)=xe needed for the b-spline representat.
|
|
! if the computation mode iopt=1 is used, the values of tx(1),
|
|
! ...,tx(nx) should be left unchanged between subsequent calls.
|
|
! if the computation mode iopt=-1 is used, the values tx(kx+2),
|
|
! ...tx(nx-kx-1) must be supplied by the user, before entry.
|
|
! see also the restrictions (ier=10).
|
|
! ny : integer.
|
|
! unless ier=10 (in case iopt >=0), ny will contain the total
|
|
! number of knots with respect to the y-variable, of the spline
|
|
! approximation returned. if the computation mode iopt=1 is
|
|
! used, the value of ny should be left unchanged between sub-
|
|
! sequent calls.
|
|
! in case iopt=-1, the value of ny should be specified on entry
|
|
! ty : real array of dimension nmax.
|
|
! on succesful exit, this array will contain the knots of the
|
|
! spline with respect to the y-variable, i.e. the position of
|
|
! the interior knots ty(ky+2),...,ty(ny-ky-1) as well as the
|
|
! position of the additional knots ty(1)=...=ty(ky+1)=yb and
|
|
! ty(ny-ky)=...=ty(ny)=ye needed for the b-spline representat.
|
|
! if the computation mode iopt=1 is used, the values of ty(1),
|
|
! ...,ty(ny) should be left unchanged between subsequent calls.
|
|
! if the computation mode iopt=-1 is used, the values ty(ky+2),
|
|
! ...ty(ny-ky-1) must be supplied by the user, before entry.
|
|
! see also the restrictions (ier=10).
|
|
! c : real array of dimension at least (nxest-kx-1)*(nyest-ky-1).
|
|
! on succesful exit, c contains the coefficients of the spline
|
|
! approximation s(x,y)
|
|
! fp : real. unless ier=10, fp contains the sum of squared
|
|
! residuals of the spline approximation returned.
|
|
! wrk : real array of dimension (lwrk). used as workspace.
|
|
! if the computation mode iopt=1 is used the values of wrk(1),
|
|
! ...,wrk(4) should be left unchanged between subsequent calls.
|
|
! lwrk : integer. on entry lwrk must specify the actual dimension of
|
|
! the array wrk as declared in the calling (sub)program.
|
|
! lwrk must not be too small.
|
|
! lwrk >= 4+nxest*(my+2*kx+5)+nyest*(2*ky+5)+mx*(kx+1)+
|
|
! my*(ky+1) +u
|
|
! where u is the larger of my and nxest.
|
|
! iwrk : integer array of dimension (kwrk). used as workspace.
|
|
! if the computation mode iopt=1 is used the values of iwrk(1),
|
|
! ...,iwrk(3) should be left unchanged between subsequent calls
|
|
! kwrk : integer. on entry kwrk must specify the actual dimension of
|
|
! the array iwrk as declared in the calling (sub)program.
|
|
! kwrk >= 3+mx+my+nxest+nyest.
|
|
! ier : integer. unless the routine detects an error, ier contains a
|
|
! non-positive value on exit, i.e.
|
|
! ier=0 : normal return. the spline returned has a residual sum of
|
|
! squares fp such that abs(fp-s)/s <= tol with tol a relat-
|
|
! ive tolerance set to 0.001 by the program.
|
|
! ier=-1 : normal return. the spline returned is an interpolating
|
|
! spline (fp=0).
|
|
! ier=-2 : normal return. the spline returned is the least-squares
|
|
! polynomial of degrees kx and ky. in this extreme case fp
|
|
! gives the upper bound for the smoothing factor s.
|
|
! ier=1 : error. the required storage space exceeds the available
|
|
! storage space, as specified by the parameters nxest and
|
|
! nyest.
|
|
! probably causes : nxest or nyest too small. if these param-
|
|
! eters are already large, it may also indicate that s is
|
|
! too small
|
|
! the approximation returned is the least-squares spline
|
|
! according to the current set of knots. the parameter fp
|
|
! gives the corresponding sum of squared residuals (fp>s).
|
|
! ier=2 : error. a theoretically impossible result was found during
|
|
! the iteration proces for finding a smoothing spline with
|
|
! fp = s. probably causes : s too small.
|
|
! there is an approximation returned but the corresponding
|
|
! sum of squared residuals does not satisfy the condition
|
|
! abs(fp-s)/s < tol.
|
|
! ier=3 : error. the maximal number of iterations maxit (set to 20
|
|
! by the program) allowed for finding a smoothing spline
|
|
! with fp=s has been reached. probably causes : s too small
|
|
! there is an approximation returned but the corresponding
|
|
! sum of squared residuals does not satisfy the condition
|
|
! abs(fp-s)/s < tol.
|
|
! ier=10 : error. on entry, the input data are controlled on validity
|
|
! the following restrictions must be satisfied.
|
|
! -1<=iopt<=1, 1<=kx,ky<=5, mx>kx, my>ky, nxest>=2*kx+2,
|
|
! nyest>=2*ky+2, kwrk>=3+mx+my+nxest+nyest,
|
|
! lwrk >= 4+nxest*(my+2*kx+5)+nyest*(2*ky+5)+mx*(kx+1)+
|
|
! my*(ky+1) +max(my,nxest),
|
|
! xb<=x(i-1)<x(i)<=xe,i=2,..,mx,yb<=y(j-1)<y(j)<=ye,j=2,..,my
|
|
! if iopt=-1: 2*kx+2<=nx<=min(nxest,mx+kx+1)
|
|
! xb<tx(kx+2)<tx(kx+3)<...<tx(nx-kx-1)<xe
|
|
! 2*ky+2<=ny<=min(nyest,my+ky+1)
|
|
! yb<ty(ky+2)<ty(ky+3)<...<ty(ny-ky-1)<ye
|
|
! the schoenberg-whitney conditions, i.e. there must
|
|
! be subset of grid co-ordinates xx(p) and yy(q) such
|
|
! that tx(p) < xx(p) < tx(p+kx+1) ,p=1,...,nx-kx-1
|
|
! ty(q) < yy(q) < ty(q+ky+1) ,q=1,...,ny-ky-1
|
|
! if iopt>=0: s>=0
|
|
! if s=0 : nxest>=mx+kx+1, nyest>=my+ky+1
|
|
! if one of these conditions is found to be violated,control
|
|
! is immediately repassed to the calling program. in that
|
|
! case there is no approximation returned.
|
|
!
|
|
! further comments:
|
|
! regrid does not allow individual weighting of the data-values.
|
|
! so, if these were determined to widely different accuracies, then
|
|
! perhaps the general data set routine surfit should rather be used
|
|
! in spite of efficiency.
|
|
! by means of the parameter s, the user can control the tradeoff
|
|
! between closeness of fit and smoothness of fit of the approximation.
|
|
! if s is too large, the spline will be too smooth and signal will be
|
|
! lost ; if s is too small the spline will pick up too much noise. in
|
|
! the extreme cases the program will return an interpolating spline if
|
|
! s=0 and the least-squares polynomial (degrees kx,ky) if s is
|
|
! very large. between these extremes, a properly chosen s will result
|
|
! in a good compromise between closeness of fit and smoothness of fit.
|
|
! to decide whether an approximation, corresponding to a certain s is
|
|
! satisfactory the user is highly recommended to inspect the fits
|
|
! graphically.
|
|
! recommended values for s depend on the accuracy of the data values.
|
|
! if the user has an idea of the statistical errors on the data, he
|
|
! can also find a proper estimate for s. for, by assuming that, if he
|
|
! specifies the right s, regrid will return a spline s(x,y) which
|
|
! exactly reproduces the function underlying the data he can evaluate
|
|
! the sum((z(i,j)-s(x(i),y(j)))**2) to find a good estimate for this s
|
|
! for example, if he knows that the statistical errors on his z(i,j)-
|
|
! values is not greater than 0.1, he may expect that a good s should
|
|
! have a value not larger than mx*my*(0.1)**2.
|
|
! if nothing is known about the statistical error in z(i,j), s must
|
|
! be determined by trial and error, taking account of the comments
|
|
! above. the best is then to start with a very large value of s (to
|
|
! determine the least-squares polynomial and the corresponding upper
|
|
! bound fp0 for s) and then to progressively decrease the value of s
|
|
! ( say by a factor 10 in the beginning, i.e. s=fp0/10,fp0/100,...
|
|
! and more carefully as the approximation shows more detail) to
|
|
! obtain closer fits.
|
|
! to economize the search for a good s-value the program provides with
|
|
! different modes of computation. at the first call of the routine, or
|
|
! whenever he wants to restart with the initial set of knots the user
|
|
! must set iopt=0.
|
|
! if iopt=1 the program will continue with the set of knots found at
|
|
! the last call of the routine. this will save a lot of computation
|
|
! time if regrid is called repeatedly for different values of s.
|
|
! the number of knots of the spline returned and their location will
|
|
! depend on the value of s and on the complexity of the shape of the
|
|
! function underlying the data. if the computation mode iopt=1
|
|
! is used, the knots returned may also depend on the s-values at
|
|
! previous calls (if these were smaller). therefore, if after a number
|
|
! of trials with different s-values and iopt=1, the user can finally
|
|
! accept a fit as satisfactory, it may be worthwhile for him to call
|
|
! regrid once more with the selected value for s but now with iopt=0.
|
|
! indeed, regrid may then return an approximation of the same quality
|
|
! of fit but with fewer knots and therefore better if data reduction
|
|
! is also an important objective for the user.
|
|
! the number of knots may also depend on the upper bounds nxest and
|
|
! nyest. indeed, if at a certain stage in regrid the number of knots
|
|
! in one direction (say nx) has reached the value of its upper bound
|
|
! (nxest), then from that moment on all subsequent knots are added
|
|
! in the other (y) direction. this may indicate that the value of
|
|
! nxest is too small. on the other hand, it gives the user the option
|
|
! of limiting the number of knots the routine locates in any direction
|
|
! for example, by setting nxest=2*kx+2 (the lowest allowable value for
|
|
! nxest), the user can indicate that he wants an approximation which
|
|
! is a simple polynomial of degree kx in the variable x.
|
|
!
|
|
! other subroutines required:
|
|
! fpback,fpbspl,fpregr,fpdisc,fpgivs,fpgrre,fprati,fprota,fpchec,
|
|
! fpknot
|
|
!
|
|
! references:
|
|
! dierckx p. : a fast algorithm for smoothing data on a rectangular
|
|
! grid while using spline functions, siam j.numer.anal.
|
|
! 19 (1982) 1286-1304.
|
|
! dierckx p. : a fast algorithm for smoothing data on a rectangular
|
|
! grid while using spline functions, report tw53, dept.
|
|
! computer science,k.u.leuven, 1980.
|
|
! dierckx p. : curve and surface fitting with splines, monographs on
|
|
! numerical analysis, oxford university press, 1993.
|
|
!
|
|
! author:
|
|
! p.dierckx
|
|
! dept. computer science, k.u. leuven
|
|
! celestijnenlaan 200a, b-3001 heverlee, belgium.
|
|
! e-mail : Paul.Dierckx@cs.kuleuven.ac.be
|
|
!
|
|
! creation date : may 1979
|
|
! latest update : march 1989
|
|
!
|
|
! ..
|
|
implicit none
|
|
! arguments
|
|
integer, intent(in) :: iopt, mx, my, kx, ky, nxest, nyest, lwrk, kwrk
|
|
integer, intent(out) :: ier
|
|
integer, intent(inout) :: nx, ny, iwrk(kwrk)
|
|
real(wp_), intent(in) :: xb, xe, yb, ye, x(mx), y(my), z(mx*my), s
|
|
real(wp_), intent(out) :: fp, c((nxest-kx-1)*(nyest-ky-1))
|
|
real(wp_), intent(inout) :: tx(nxest), ty(nyest), wrk(lwrk)
|
|
! local variables
|
|
integer :: i, j, jwrk, kndx, kndy, knrx, knry, kwest, kx1, kx2, &
|
|
ky1, ky2, lfpx, lfpy, lwest, lww, nc, nminx, nminy, mz
|
|
! parameters
|
|
integer, parameter :: maxit = 20
|
|
real(wp_), parameter :: tol = 0.1e-02_wp_
|
|
! ..subroutine references..
|
|
! fpregr,fpchec
|
|
! ..
|
|
! before starting computations a data check is made. if the input data
|
|
! are invalid, control is immediately repassed to the calling program.
|
|
ier = 10
|
|
if(kx<=0 .or. kx>5) return
|
|
kx1 = kx+1
|
|
kx2 = kx1+1
|
|
if(ky<=0 .or. ky>5) return
|
|
ky1 = ky+1
|
|
ky2 = ky1+1
|
|
if(iopt<(-1) .or. iopt>1) return
|
|
nminx = 2*kx1
|
|
if(mx<kx1 .or. nxest<nminx) return
|
|
nminy = 2*ky1
|
|
if(my<ky1 .or. nyest<nminy) return
|
|
mz = mx*my
|
|
nc = (nxest-kx1)*(nyest-ky1)
|
|
lwest = 4+nxest*(my+2*kx2+1)+nyest*(2*ky2+1)+mx*kx1+ &
|
|
my*ky1+max(nxest,my)
|
|
kwest = 3+mx+my+nxest+nyest
|
|
if(lwrk<lwest .or. kwrk<kwest) return
|
|
if(xb>x(1) .or. xe<x(mx)) return
|
|
do i=2,mx
|
|
if(x(i-1)>=x(i)) return
|
|
end do
|
|
if(yb>y(1) .or. ye<y(my)) return
|
|
do i=2,my
|
|
if(y(i-1)>=y(i)) return
|
|
end do
|
|
if(iopt<0) then
|
|
if(nx<nminx .or. nx>nxest) return
|
|
j = nx
|
|
do i=1,kx1
|
|
tx(i) = xb
|
|
tx(j) = xe
|
|
j = j-1
|
|
end do
|
|
call fpchec(x,mx,tx,nx,kx,ier)
|
|
if(ier/=0) return
|
|
if(ny<nminy .or. ny>nyest) return
|
|
j = ny
|
|
do i=1,ky1
|
|
ty(i) = yb
|
|
ty(j) = ye
|
|
j = j-1
|
|
end do
|
|
call fpchec(y,my,ty,ny,ky,ier)
|
|
if(ier/=0) return
|
|
else
|
|
if(s<0.0_wp_) return
|
|
if(s==0.0_wp_ .and. (nxest<(mx+kx1) .or. nyest<(my+ky1)) ) &
|
|
return
|
|
ier = 0
|
|
end if
|
|
! we partition the working space and determine the spline approximation
|
|
lfpx = 5
|
|
lfpy = lfpx+nxest
|
|
lww = lfpy+nyest
|
|
jwrk = lwrk-4-nxest-nyest
|
|
knrx = 4
|
|
knry = knrx+mx
|
|
kndx = knry+my
|
|
kndy = kndx+nxest
|
|
call fpregr(iopt,x,mx,y,my,z,mz,xb,xe,yb,ye,kx,ky,s,nxest,nyest, &
|
|
tol,maxit,nc,nx,tx,ny,ty,c,fp,wrk(1),wrk(2),wrk(3),wrk(4), &
|
|
wrk(lfpx),wrk(lfpy),iwrk(1),iwrk(2),iwrk(3),iwrk(knrx), &
|
|
iwrk(knry),iwrk(kndx),iwrk(kndy),wrk(lww),jwrk,ier)
|
|
end subroutine regrid
|
|
|
|
subroutine parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,x,mx,y,my,z, &
|
|
wrk,lwrk,iwrk,kwrk,ier)
|
|
! subroutine parder evaluates on a grid (x(i),y(j)),i=1,...,mx; j=1,...
|
|
! ,my the partial derivative ( order nux,nuy) of a bivariate spline
|
|
! s(x,y) of degrees kx and ky, given in the b-spline representation.
|
|
!
|
|
! calling sequence:
|
|
! call parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,x,mx,y,my,z,wrk,lwrk,
|
|
! * iwrk,kwrk,ier)
|
|
!
|
|
! input parameters:
|
|
! tx : real array, length nx, which contains the position of the
|
|
! knots in the x-direction.
|
|
! nx : integer, giving the total number of knots in the x-direction
|
|
! ty : real array, length ny, which contains the position of the
|
|
! knots in the y-direction.
|
|
! ny : integer, giving the total number of knots in the y-direction
|
|
! c : real array, length (nx-kx-1)*(ny-ky-1), which contains the
|
|
! b-spline coefficients.
|
|
! kx,ky : integer values, giving the degrees of the spline.
|
|
! nux : integer values, specifying the order of the partial
|
|
! nuy derivative. 0<=nux<kx, 0<=nuy<ky.
|
|
! x : real array of dimension (mx).
|
|
! before entry x(i) must be set to the x co-ordinate of the
|
|
! i-th grid point along the x-axis.
|
|
! tx(kx+1)<=x(i-1)<=x(i)<=tx(nx-kx), i=2,...,mx.
|
|
! mx : on entry mx must specify the number of grid points along
|
|
! the x-axis. mx >=1.
|
|
! y : real array of dimension (my).
|
|
! before entry y(j) must be set to the y co-ordinate of the
|
|
! j-th grid point along the y-axis.
|
|
! ty(ky+1)<=y(j-1)<=y(j)<=ty(ny-ky), j=2,...,my.
|
|
! my : on entry my must specify the number of grid points along
|
|
! the y-axis. my >=1.
|
|
! wrk : real array of dimension lwrk. used as workspace.
|
|
! lwrk : integer, specifying the dimension of wrk.
|
|
! lwrk >= mx*(kx+1-nux)+my*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1)
|
|
! iwrk : integer array of dimension kwrk. used as workspace.
|
|
! kwrk : integer, specifying the dimension of iwrk. kwrk >= mx+my.
|
|
!
|
|
! output parameters:
|
|
! z : real array of dimension (mx*my).
|
|
! on succesful exit z(my*(i-1)+j) contains the value of the
|
|
! specified partial derivative of s(x,y) at the point
|
|
! (x(i),y(j)),i=1,...,mx;j=1,...,my.
|
|
! ier : integer error flag
|
|
! ier=0 : normal return
|
|
! ier=10: invalid input data (see restrictions)
|
|
!
|
|
! restrictions:
|
|
! mx >=1, my >=1, 0 <= nux < kx, 0 <= nuy < ky, kwrk>=mx+my
|
|
! lwrk>=mx*(kx+1-nux)+my*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1),
|
|
! tx(kx+1) <= x(i-1) <= x(i) <= tx(nx-kx), i=2,...,mx
|
|
! ty(ky+1) <= y(j-1) <= y(j) <= ty(ny-ky), j=2,...,my
|
|
!
|
|
! other subroutines required:
|
|
! fpbisp,fpbspl
|
|
!
|
|
! references :
|
|
! de boor c : on calculating with b-splines, j. approximation theory
|
|
! 6 (1972) 50-62.
|
|
! dierckx p. : curve and surface fitting with splines, monographs on
|
|
! numerical analysis, oxford university press, 1993.
|
|
!
|
|
! author :
|
|
! p.dierckx
|
|
! dept. computer science, k.u.leuven
|
|
! celestijnenlaan 200a, b-3001 heverlee, belgium.
|
|
! e-mail : Paul.Dierckx@cs.kuleuven.ac.be
|
|
!
|
|
! latest update : march 1989
|
|
!
|
|
implicit none
|
|
! arguments
|
|
integer, intent(in) :: nx, ny, kx, ky, nux, nuy, mx, my, lwrk, kwrk
|
|
integer, intent(out) :: ier
|
|
integer, intent(inout) :: iwrk(kwrk)
|
|
real(wp_), intent(in) :: tx(nx), ty(ny), c((nx-kx-1)*(ny-ky-1)), &
|
|
x(mx), y(my)
|
|
real(wp_), intent(out) :: z(mx*my)
|
|
real(wp_), intent(inout) :: wrk(lwrk)
|
|
! local variables
|
|
integer :: i, iwx, iwy, j, kkx, kky, kx1, ky1, lx, ly, lwest, &
|
|
l1, l2, m, m0, m1, nc, nkx1, nky1, nxx, nyy
|
|
real(wp_) :: ak, fac
|
|
! ..
|
|
! before starting computations a data check is made. if the input data
|
|
! are invalid control is immediately repassed to the calling program.
|
|
ier = 10
|
|
kx1 = kx+1
|
|
ky1 = ky+1
|
|
nkx1 = nx-kx1
|
|
nky1 = ny-ky1
|
|
nc = nkx1*nky1
|
|
if(nux<0 .or. nux>=kx) return
|
|
if(nuy<0 .or. nuy>=ky) return
|
|
lwest = nc +(kx1-nux)*mx+(ky1-nuy)*my
|
|
if(lwrk<lwest) return
|
|
if(kwrk<(mx+my)) return
|
|
if(mx<1) return
|
|
do i=2,mx
|
|
if(x(i)<x(i-1)) return
|
|
end do
|
|
if(my<1) return
|
|
do i=2,my
|
|
if(y(i)<y(i-1)) return
|
|
end do
|
|
ier = 0
|
|
nxx = nkx1
|
|
nyy = nky1
|
|
kkx = kx
|
|
kky = ky
|
|
! the partial derivative of order (nux,nuy) of a bivariate spline of
|
|
! degrees kx,ky is a bivariate spline of degrees kx-nux,ky-nuy.
|
|
! we calculate the b-spline coefficients of this spline
|
|
do i=1,nc
|
|
wrk(i) = c(i)
|
|
end do
|
|
if(nux/=0) then
|
|
lx = 1
|
|
do j=1,nux
|
|
ak = kkx
|
|
nxx = nxx-1
|
|
l1 = lx
|
|
m0 = 1
|
|
do i=1,nxx
|
|
l1 = l1+1
|
|
l2 = l1+kkx
|
|
fac = tx(l2)-tx(l1)
|
|
if(fac>0.0_wp_) then
|
|
do m=1,nyy
|
|
m1 = m0+nyy
|
|
wrk(m0) = (wrk(m1)-wrk(m0))*ak/fac
|
|
m0 = m0+1
|
|
end do
|
|
end if
|
|
end do
|
|
lx = lx+1
|
|
kkx = kkx-1
|
|
end do
|
|
end if
|
|
if(nuy/=0) then
|
|
ly = 1
|
|
do j=1,nuy
|
|
ak = kky
|
|
nyy = nyy-1
|
|
l1 = ly
|
|
do i=1,nyy
|
|
l1 = l1+1
|
|
l2 = l1+kky
|
|
fac = ty(l2)-ty(l1)
|
|
if(fac>0.0_wp_) then
|
|
m0 = i
|
|
do m=1,nxx
|
|
m1 = m0+1
|
|
wrk(m0) = (wrk(m1)-wrk(m0))*ak/fac
|
|
m0 = m0+nky1
|
|
end do
|
|
end if
|
|
end do
|
|
ly = ly+1
|
|
kky = kky-1
|
|
end do
|
|
m0 = nyy
|
|
m1 = nky1
|
|
do m=2,nxx
|
|
do i=1,nyy
|
|
m0 = m0+1
|
|
m1 = m1+1
|
|
wrk(m0) = wrk(m1)
|
|
end do
|
|
m1 = m1+nuy
|
|
end do
|
|
end if
|
|
! we partition the working space and evaluate the partial derivative
|
|
iwx = 1+nxx*nyy
|
|
iwy = iwx+mx*(kx1-nux)
|
|
call fpbisp(tx(nux+1),nx-2*nux,ty(nuy+1),ny-2*nuy,wrk,kkx,kky, &
|
|
x,mx,y,my,z,wrk(iwx),wrk(iwy),iwrk(1),iwrk(mx+1))
|
|
end subroutine parder
|
|
|
|
subroutine coeff_parder(tx,nx,ty,ny,c,kx,ky,nux,nuy, &
|
|
wrk,lwrk,ier)
|
|
! subroutine parder evaluates on a grid (x(i),y(j)),i=1,...,mx; j=1,...
|
|
! ,my the partial derivative ( order nux,nuy) of a bivariate spline
|
|
! s(x,y) of degrees kx and ky, given in the b-spline representation.
|
|
!
|
|
! calling sequence:
|
|
! call parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,x,mx,y,my,z,wrk,lwrk,
|
|
! * iwrk,kwrk,ier)
|
|
!
|
|
! input parameters:
|
|
! tx : real array, length nx, which contains the position of the
|
|
! knots in the x-direction.
|
|
! nx : integer, giving the total number of knots in the x-direction
|
|
! ty : real array, length ny, which contains the position of the
|
|
! knots in the y-direction.
|
|
! ny : integer, giving the total number of knots in the y-direction
|
|
! c : real array, length (nx-kx-1)*(ny-ky-1), which contains the
|
|
! b-spline coefficients.
|
|
! kx,ky : integer values, giving the degrees of the spline.
|
|
! nux : integer values, specifying the order of the partial
|
|
! nuy derivative. 0<=nux<kx, 0<=nuy<ky.
|
|
! wrk : real array of dimension lwrk. used as workspace.
|
|
! lwrk : integer, specifying the dimension of wrk.
|
|
! lwrk >= mx*(kx+1-nux)+my*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1)
|
|
!
|
|
! output parameters:
|
|
! ier : integer error flag
|
|
! ier=0 : normal return
|
|
! ier=10: invalid input data (see restrictions)
|
|
!
|
|
! restrictions:
|
|
! mx >=1, my >=1, 0 <= nux < kx, 0 <= nuy < ky
|
|
! lwrk>=mx*(kx+1-nux)+my*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1),
|
|
!
|
|
! other subroutines required:
|
|
! fpbisp,fpbspl
|
|
!
|
|
! references :
|
|
! de boor c : on calculating with b-splines, j. approximation theory
|
|
! 6 (1972) 50-62.
|
|
! dierckx p. : curve and surface fitting with splines, monographs on
|
|
! numerical analysis, oxford university press, 1993.
|
|
!
|
|
! author :
|
|
! p.dierckx
|
|
! dept. computer science, k.u.leuven
|
|
! celestijnenlaan 200a, b-3001 heverlee, belgium.
|
|
! e-mail : Paul.Dierckx@cs.kuleuven.ac.be
|
|
!
|
|
! latest update : march 1989
|
|
!
|
|
implicit none
|
|
! arguments
|
|
integer, intent(in) :: nx, ny, kx, ky, nux, nuy, lwrk
|
|
integer, intent(out) :: ier
|
|
real(wp_), intent(in) :: tx(nx), ty(ny), c((nx-kx-1)*(ny-ky-1))
|
|
real(wp_), intent(inout) :: wrk(lwrk)
|
|
! local variables
|
|
integer :: mx, my, i, j, kkx, kky, kx1, ky1, lx, ly, lwest, &
|
|
l1, l2, m, m0, m1, nc, nkx1, nky1, nxx, nyy
|
|
real(wp_) :: ak, fac
|
|
! ..
|
|
! before starting computations a data check is made. if the input data
|
|
! are invalid control is immediately repassed to the calling program.
|
|
mx=1
|
|
my=1
|
|
ier = 10
|
|
kx1 = kx+1
|
|
ky1 = ky+1
|
|
nkx1 = nx-kx1
|
|
nky1 = ny-ky1
|
|
nc = nkx1*nky1
|
|
if(nux<0 .or. nux>=kx) return
|
|
if(nuy<0 .or. nuy>=ky) return
|
|
lwest = nc +(kx1-nux)*mx+(ky1-nuy)*my
|
|
if(lwrk<lwest) return
|
|
ier = 0
|
|
nxx = nkx1
|
|
nyy = nky1
|
|
kkx = kx
|
|
kky = ky
|
|
! the partial derivative of order (nux,nuy) of a bivariate spline of
|
|
! degrees kx,ky is a bivariate spline of degrees kx-nux,ky-nuy.
|
|
! we calculate the b-spline coefficients of this spline
|
|
do i=1,nc
|
|
wrk(i) = c(i)
|
|
end do
|
|
if(nux/=0) then
|
|
lx = 1
|
|
do j=1,nux
|
|
ak = kkx
|
|
nxx = nxx-1
|
|
l1 = lx
|
|
m0 = 1
|
|
do i=1,nxx
|
|
l1 = l1+1
|
|
l2 = l1+kkx
|
|
fac = tx(l2)-tx(l1)
|
|
if(fac>0.0_wp_) then
|
|
do m=1,nyy
|
|
m1 = m0+nyy
|
|
wrk(m0) = (wrk(m1)-wrk(m0))*ak/fac
|
|
m0 = m0+1
|
|
end do
|
|
end if
|
|
end do
|
|
lx = lx+1
|
|
kkx = kkx-1
|
|
end do
|
|
end if
|
|
if(nuy/=0) then
|
|
ly = 1
|
|
do j=1,nuy
|
|
ak = kky
|
|
nyy = nyy-1
|
|
l1 = ly
|
|
do i=1,nyy
|
|
l1 = l1+1
|
|
l2 = l1+kky
|
|
fac = ty(l2)-ty(l1)
|
|
if(fac>0.0_wp_) then
|
|
m0 = i
|
|
do m=1,nxx
|
|
m1 = m0+1
|
|
wrk(m0) = (wrk(m1)-wrk(m0))*ak/fac
|
|
m0 = m0+nky1
|
|
end do
|
|
end if
|
|
end do
|
|
ly = ly+1
|
|
kky = kky-1
|
|
end do
|
|
m0 = nyy
|
|
m1 = nky1
|
|
do m=2,nxx
|
|
do i=1,nyy
|
|
m0 = m0+1
|
|
m1 = m1+1
|
|
wrk(m0) = wrk(m1)
|
|
end do
|
|
m1 = m1+nuy
|
|
end do
|
|
end if
|
|
end subroutine coeff_parder
|
|
|
|
subroutine curfit(iopt,m,x,y,w,xb,xe,k,s,nest,n,t,c,fp, &
|
|
wrk,lwrk,iwrk,ier)
|
|
! given the set of data points (x(i),y(i)) and the set of positive
|
|
! numbers w(i),i=1,2,...,m,subroutine curfit determines a smooth spline
|
|
! approximation of degree k on the interval xb <= x <= xe.
|
|
! if iopt=-1 curfit calculates the weighted least-squares spline
|
|
! according to a given set of knots.
|
|
! if iopt>=0 the number of knots of the spline s(x) and the position
|
|
! t(j),j=1,2,...,n is chosen automatically by the routine. the smooth-
|
|
! ness of s(x) is then achieved by minimalizing the discontinuity
|
|
! jumps of the k-th derivative of s(x) at the knots t(j),j=k+2,k+3,...,
|
|
! n-k-1. the amount of smoothness is determined by the condition that
|
|
! f(p)=sum((w(i)*(y(i)-s(x(i))))**2) be <= s, with s a given non-
|
|
! negative constant, called the smoothing factor.
|
|
! the fit s(x) is given in the b-spline representation (b-spline coef-
|
|
! ficients c(j),j=1,2,...,n-k-1) and can be evaluated by means of
|
|
! subroutine splev.
|
|
!
|
|
! calling sequence:
|
|
! call curfit(iopt,m,x,y,w,xb,xe,k,s,nest,n,t,c,fp,wrk,
|
|
! * lwrk,iwrk,ier)
|
|
!
|
|
! parameters:
|
|
! iopt : integer flag. on entry iopt must specify whether a weighted
|
|
! least-squares spline (iopt=-1) or a smoothing spline (iopt=
|
|
! 0 or 1) must be determined. if iopt=0 the routine will start
|
|
! with an initial set of knots t(i)=xb, t(i+k+1)=xe, i=1,2,...
|
|
! k+1. if iopt=1 the routine will continue with the knots
|
|
! found at the last call of the routine.
|
|
! attention: a call with iopt=1 must always be immediately
|
|
! preceded by another call with iopt=1 or iopt=0.
|
|
! unchanged on exit.
|
|
! m : integer. on entry m must specify the number of data points.
|
|
! m > k. unchanged on exit.
|
|
! x : real array of dimension at least (m). before entry, x(i)
|
|
! must be set to the i-th value of the independent variable x,
|
|
! for i=1,2,...,m. these values must be supplied in strictly
|
|
! ascending order. unchanged on exit.
|
|
! y : real array of dimension at least (m). before entry, y(i)
|
|
! must be set to the i-th value of the dependent variable y,
|
|
! for i=1,2,...,m. unchanged on exit.
|
|
! w : real array of dimension at least (m). before entry, w(i)
|
|
! must be set to the i-th value in the set of weights. the
|
|
! w(i) must be strictly positive. unchanged on exit.
|
|
! see also further comments.
|
|
! xb,xe : real values. on entry xb and xe must specify the boundaries
|
|
! of the approximation interval. xb<=x(1), xe>=x(m).
|
|
! unchanged on exit.
|
|
! k : integer. on entry k must specify the degree of the spline.
|
|
! 1<=k<=5. it is recommended to use cubic splines (k=3).
|
|
! the user is strongly dissuaded from choosing k even,together
|
|
! with a small s-value. unchanged on exit.
|
|
! s : real.on entry (in case iopt>=0) s must specify the smoothing
|
|
! factor. s >=0. unchanged on exit.
|
|
! for advice on the choice of s see further comments.
|
|
! nest : integer. on entry nest must contain an over-estimate of the
|
|
! total number of knots of the spline returned, to indicate
|
|
! the storage space available to the routine. nest >=2*k+2.
|
|
! in most practical situation nest=m/2 will be sufficient.
|
|
! always large enough is nest=m+k+1, the number of knots
|
|
! needed for interpolation (s=0). unchanged on exit.
|
|
! n : integer.
|
|
! unless ier =10 (in case iopt >=0), n will contain the
|
|
! total number of knots of the spline approximation returned.
|
|
! if the computation mode iopt=1 is used this value of n
|
|
! should be left unchanged between subsequent calls.
|
|
! in case iopt=-1, the value of n must be specified on entry.
|
|
! t : real array of dimension at least (nest).
|
|
! on succesful exit, this array will contain the knots of the
|
|
! spline,i.e. the position of the interior knots t(k+2),t(k+3)
|
|
! ...,t(n-k-1) as well as the position of the additional knots
|
|
! t(1)=t(2)=...=t(k+1)=xb and t(n-k)=...=t(n)=xe needed for
|
|
! the b-spline representation.
|
|
! if the computation mode iopt=1 is used, the values of t(1),
|
|
! t(2),...,t(n) should be left unchanged between subsequent
|
|
! calls. if the computation mode iopt=-1 is used, the values
|
|
! t(k+2),...,t(n-k-1) must be supplied by the user, before
|
|
! entry. see also the restrictions (ier=10).
|
|
! c : real array of dimension at least (nest).
|
|
! on succesful exit, this array will contain the coefficients
|
|
! c(1),c(2),..,c(n-k-1) in the b-spline representation of s(x)
|
|
! fp : real. unless ier=10, fp contains the weighted sum of
|
|
! squared residuals of the spline approximation returned.
|
|
! wrk : real array of dimension at least (m*(k+1)+nest*(7+3*k)).
|
|
! used as working space. if the computation mode iopt=1 is
|
|
! used, the values wrk(1),...,wrk(n) should be left unchanged
|
|
! between subsequent calls.
|
|
! lwrk : integer. on entry,lwrk must specify the actual dimension of
|
|
! the array wrk as declared in the calling (sub)program.lwrk
|
|
! must not be too small (see wrk). unchanged on exit.
|
|
! iwrk : integer array of dimension at least (nest).
|
|
! used as working space. if the computation mode iopt=1 is
|
|
! used,the values iwrk(1),...,iwrk(n) should be left unchanged
|
|
! between subsequent calls.
|
|
! ier : integer. unless the routine detects an error, ier contains a
|
|
! non-positive value on exit, i.e.
|
|
! ier=0 : normal return. the spline returned has a residual sum of
|
|
! squares fp such that abs(fp-s)/s <= tol with tol a relat-
|
|
! ive tolerance set to 0.001 by the program.
|
|
! ier=-1 : normal return. the spline returned is an interpolating
|
|
! spline (fp=0).
|
|
! ier=-2 : normal return. the spline returned is the weighted least-
|
|
! squares polynomial of degree k. in this extreme case fp
|
|
! gives the upper bound fp0 for the smoothing factor s.
|
|
! ier=1 : error. the required storage space exceeds the available
|
|
! storage space, as specified by the parameter nest.
|
|
! probably causes : nest too small. if nest is already
|
|
! large (say nest > m/2), it may also indicate that s is
|
|
! too small
|
|
! the approximation returned is the weighted least-squares
|
|
! spline according to the knots t(1),t(2),...,t(n). (n=nest)
|
|
! the parameter fp gives the corresponding weighted sum of
|
|
! squared residuals (fp>s).
|
|
! ier=2 : error. a theoretically impossible result was found during
|
|
! the iteration proces for finding a smoothing spline with
|
|
! fp = s. probably causes : s too small.
|
|
! there is an approximation returned but the corresponding
|
|
! weighted sum of squared residuals does not satisfy the
|
|
! condition abs(fp-s)/s < tol.
|
|
! ier=3 : error. the maximal number of iterations maxit (set to 20
|
|
! by the program) allowed for finding a smoothing spline
|
|
! with fp=s has been reached. probably causes : s too small
|
|
! there is an approximation returned but the corresponding
|
|
! weighted sum of squared residuals does not satisfy the
|
|
! condition abs(fp-s)/s < tol.
|
|
! ier=10 : error. on entry, the input data are controlled on validity
|
|
! the following restrictions must be satisfied.
|
|
! -1<=iopt<=1, 1<=k<=5, m>k, nest>2*k+2, w(i)>0,i=1,2,...,m
|
|
! xb<=x(1)<x(2)<...<x(m)<=xe, lwrk>=(k+1)*m+nest*(7+3*k)
|
|
! if iopt=-1: 2*k+2<=n<=min(nest,m+k+1)
|
|
! xb<t(k+2)<t(k+3)<...<t(n-k-1)<xe
|
|
! the schoenberg-whitney conditions, i.e. there
|
|
! must be a subset of data points xx(j) such that
|
|
! t(j) < xx(j) < t(j+k+1), j=1,2,...,n-k-1
|
|
! if iopt>=0: s>=0
|
|
! if s=0 : nest >= m+k+1
|
|
! if one of these conditions is found to be violated,control
|
|
! is immediately repassed to the calling program. in that
|
|
! case there is no approximation returned.
|
|
!
|
|
! further comments:
|
|
! by means of the parameter s, the user can control the tradeoff
|
|
! between closeness of fit and smoothness of fit of the approximation.
|
|
! if s is too large, the spline will be too smooth and signal will be
|
|
! lost ; if s is too small the spline will pick up too much noise. in
|
|
! the extreme cases the program will return an interpolating spline if
|
|
! s=0 and the weighted least-squares polynomial of degree k if s is
|
|
! very large. between these extremes, a properly chosen s will result
|
|
! in a good compromise between closeness of fit and smoothness of fit.
|
|
! to decide whether an approximation, corresponding to a certain s is
|
|
! satisfactory the user is highly recommended to inspect the fits
|
|
! graphically.
|
|
! recommended values for s depend on the weights w(i). if these are
|
|
! taken as 1/d(i) with d(i) an estimate of the standard deviation of
|
|
! y(i), a good s-value should be found in the range (m-sqrt(2*m),m+
|
|
! sqrt(2*m)). if nothing is known about the statistical error in y(i)
|
|
! each w(i) can be set equal to one and s determined by trial and
|
|
! error, taking account of the comments above. the best is then to
|
|
! start with a very large value of s ( to determine the least-squares
|
|
! polynomial and the corresponding upper bound fp0 for s) and then to
|
|
! progressively decrease the value of s ( say by a factor 10 in the
|
|
! beginning, i.e. s=fp0/10, fp0/100,...and more carefully as the
|
|
! approximation shows more detail) to obtain closer fits.
|
|
! to economize the search for a good s-value the program provides with
|
|
! different modes of computation. at the first call of the routine, or
|
|
! whenever he wants to restart with the initial set of knots the user
|
|
! must set iopt=0.
|
|
! if iopt=1 the program will continue with the set of knots found at
|
|
! the last call of the routine. this will save a lot of computation
|
|
! time if curfit is called repeatedly for different values of s.
|
|
! the number of knots of the spline returned and their location will
|
|
! depend on the value of s and on the complexity of the shape of the
|
|
! function underlying the data. but, if the computation mode iopt=1
|
|
! is used, the knots returned may also depend on the s-values at
|
|
! previous calls (if these were smaller). therefore, if after a number
|
|
! of trials with different s-values and iopt=1, the user can finally
|
|
! accept a fit as satisfactory, it may be worthwhile for him to call
|
|
! curfit once more with the selected value for s but now with iopt=0.
|
|
! indeed, curfit may then return an approximation of the same quality
|
|
! of fit but with fewer knots and therefore better if data reduction
|
|
! is also an important objective for the user.
|
|
!
|
|
! other subroutines required:
|
|
! fpback,fpbspl,fpchec,fpcurf,fpdisc,fpgivs,fpknot,fprati,fprota
|
|
!
|
|
! references:
|
|
! dierckx p. : an algorithm for smoothing, differentiation and integ-
|
|
! ration of experimental data using spline functions,
|
|
! j.comp.appl.maths 1 (1975) 165-184.
|
|
! dierckx p. : a fast algorithm for smoothing data on a rectangular
|
|
! grid while using spline functions, siam j.numer.anal.
|
|
! 19 (1982) 1286-1304.
|
|
! dierckx p. : an improved algorithm for curve fitting with spline
|
|
! functions, report tw54, dept. computer science,k.u.
|
|
! leuven, 1981.
|
|
! dierckx p. : curve and surface fitting with splines, monographs on
|
|
! numerical analysis, oxford university press, 1993.
|
|
!
|
|
! author:
|
|
! p.dierckx
|
|
! dept. computer science, k.u. leuven
|
|
! celestijnenlaan 200a, b-3001 heverlee, belgium.
|
|
! e-mail : Paul.Dierckx@cs.kuleuven.ac.be
|
|
!
|
|
! creation date : may 1979
|
|
! latest update : march 1987
|
|
!
|
|
! ..
|
|
implicit none
|
|
! arguments
|
|
integer, intent(in) :: iopt, m, k, nest, lwrk
|
|
integer, intent(out) :: ier
|
|
integer, intent(inout) :: n, iwrk(nest)
|
|
real(wp_), intent(in) :: xb, xe, s, x(m), y(m), w(m)
|
|
real(wp_), intent(out) :: fp, c(nest)
|
|
real(wp_), intent(inout) :: t(nest), wrk(lwrk)
|
|
! local variables
|
|
integer :: i, ia, ib, ifp, ig, iq, iz, j, k1, k2, lwest, nmin
|
|
! parameters
|
|
integer, parameter :: maxit = 20
|
|
real(wp_), parameter :: tol = 0.1e-02_wp_
|
|
! ..
|
|
! before starting computations a data check is made. if the input data
|
|
! are invalid, control is immediately repassed to the calling program.
|
|
ier = 10
|
|
if(k<=0 .or. k>5) return
|
|
k1 = k+1
|
|
k2 = k1+1
|
|
if(iopt<(-1) .or. iopt>1) return
|
|
nmin = 2*k1
|
|
if(m<k1 .or. nest<nmin) return
|
|
lwest = m*k1+nest*(7+3*k)
|
|
if(lwrk<lwest) return
|
|
if(xb>x(1) .or. xe<x(m) .or. w(1)<=0.0_wp_) return
|
|
do i=2,m
|
|
if(x(i-1)>=x(i) .or. w(i)<=0.0_wp_) return
|
|
end do
|
|
if(iopt<0) then
|
|
if(n<nmin .or. n>nest) return
|
|
j = n
|
|
do i=1,k1
|
|
t(i) = xb
|
|
t(j) = xe
|
|
j = j-1
|
|
end do
|
|
call fpchec(x,m,t,n,k,ier)
|
|
if(ier/=0) return
|
|
else
|
|
if(s<0.0_wp_) return
|
|
if(s==0.0_wp_ .and. nest<(m+k1)) return
|
|
ier = 0
|
|
end if
|
|
! we partition the working space and determine the spline approximation.
|
|
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)
|
|
end subroutine curfit
|
|
|
|
subroutine fpcurf(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2, &
|
|
n,t,c,fp,fpint,z,a,b,g,q,nrdata,ier)
|
|
! ..
|
|
implicit none
|
|
! arguments
|
|
integer, intent(in) :: iopt, m, k, nest, maxit, k1, k2
|
|
integer, intent(out) :: ier
|
|
integer, intent(inout) :: n, nrdata(nest)
|
|
real(wp_), intent(in) :: xb, xe, s, tol, x(m), y(m), w(m)
|
|
real(wp_), intent(out) :: fp, c(nest)
|
|
real(wp_), intent(inout) :: t(nest), fpint(nest), z(nest), &
|
|
a(nest,k1), b(nest,k2), g(nest,k2), q(m,k1)
|
|
! local variables
|
|
real(wp_) :: acc, cs, fpart, fpms, fpold, fp0, f1, f2, f3, p, pinv, &
|
|
piv, p1, p2, p3, rn, sn, store, term, wi, xi, yi, h(7)
|
|
integer :: i, ich1, ich3, it, iter, i1, i2, i3, j, k3, l, l0, mk1, new, &
|
|
nk1, nmax, nmin, nplus, npl1, nrint, n8
|
|
logical :: rstart
|
|
! parameters
|
|
real(wp_), parameter :: one=1.0_wp_,con1=0.1_wp_,con9=0.9_wp_,con4=0.4e-01_wp_,half=0.5_wp_
|
|
! ..function references
|
|
! real(8) abs
|
|
! integer max0,min0
|
|
! ..subroutine references..
|
|
! fpback,fpbspl,fpgivs,fpdisc,fpknot,fprota
|
|
! ..
|
|
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
! part 1: determination of the number of knots and their position c
|
|
! ************************************************************** c
|
|
! given a set of knots we compute the least-squares spline sinf(x), c
|
|
! and the corresponding sum of squared residuals fp=f(p=inf). c
|
|
! if iopt=-1 sinf(x) is the requested approximation. c
|
|
! if iopt=0 or iopt=1 we check whether we can accept the knots: c
|
|
! if fp <=s we will continue with the current set of knots. c
|
|
! if fp > s we will increase the number of knots and compute the c
|
|
! corresponding least-squares spline until finally fp<=s. c
|
|
! the initial choice of knots depends on the value of s and iopt. c
|
|
! if s=0 we have spline interpolation; in that case the number of c
|
|
! knots equals nmax = m+k+1. c
|
|
! if s > 0 and c
|
|
! iopt=0 we first compute the least-squares polynomial of c
|
|
! degree k; n = nmin = 2*k+2 c
|
|
! iopt=1 we start with the set of knots found at the last c
|
|
! call of the routine, except for the case that s > fp0; then c
|
|
! we compute directly the least-squares polynomial of degree k. c
|
|
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
! determine nmin, the number of knots for polynomial approximation.
|
|
nmin = 2*k1
|
|
if(iopt>=0) then
|
|
! calculation of acc, the absolute tolerance for the root of f(p)=s.
|
|
acc = tol*s
|
|
! determine nmax, the number of knots for spline interpolation.
|
|
nmax = m+k1
|
|
if(s<=0.0_wp_) then
|
|
! if s=0, s(x) is an interpolating spline.
|
|
! test whether the required storage space exceeds the available one.
|
|
n = nmax
|
|
if(nmax>nest) then
|
|
ier = 1
|
|
return
|
|
end if
|
|
! find the position of the interior knots in case of interpolation.
|
|
mk1 = m-k1
|
|
if(mk1/=0) then
|
|
k3 = k/2
|
|
i = k2
|
|
j = k3+2
|
|
if(k3*2/=k) then
|
|
do l=1,mk1
|
|
t(i) = x(j)
|
|
i = i+1
|
|
j = j+1
|
|
end do
|
|
else
|
|
do l=1,mk1
|
|
t(i) = (x(j)+x(j-1))*half
|
|
i = i+1
|
|
j = j+1
|
|
end do
|
|
end if
|
|
end if
|
|
else
|
|
! if s>0 our initial choice of knots depends on the value of iopt.
|
|
! if iopt=0 or iopt=1 and s>=fp0, we start computing the least-squares
|
|
! polynomial of degree k which is a spline without interior knots.
|
|
! if iopt=1 and fp0>s we start computing the least squares spline
|
|
! according to the set of knots found at the last call of the routine.
|
|
if (iopt==0 .or. n==nmin) then
|
|
n = nmin
|
|
fpold = 0.0_wp_
|
|
nplus = 0
|
|
nrdata(1) = m-2
|
|
else
|
|
fp0 = fpint(n)
|
|
if (s>=fp0) then
|
|
n = nmin
|
|
fpold = 0.0_wp_
|
|
nplus = 0
|
|
nrdata(1) = m-2
|
|
else
|
|
fpold = fpint(n-1)
|
|
nplus = nrdata(n)
|
|
end if
|
|
end if
|
|
end if
|
|
end if
|
|
do
|
|
rstart=.false.
|
|
! main loop for the different sets of knots. m is a save upper bound
|
|
! for the number of trials.
|
|
do iter = 1,m
|
|
if(n==nmin) ier = -2
|
|
! find nrint, tne number of knot intervals.
|
|
nrint = n-nmin+1
|
|
! find the position of the additional knots which are needed for
|
|
! the b-spline representation of s(x).
|
|
nk1 = n-k1
|
|
i = n
|
|
do j=1,k1
|
|
t(j) = xb
|
|
t(i) = xe
|
|
i = i-1
|
|
end do
|
|
! compute the b-spline coefficients of the least-squares spline
|
|
! sinf(x). the observation matrix a is built up row by row and
|
|
! reduced to upper triangular form by givens transformations.
|
|
! at the same time fp=f(p=inf) is computed.
|
|
fp = 0.0_wp_
|
|
! initialize the observation matrix a.
|
|
do i=1,nk1
|
|
z(i) = 0.0_wp_
|
|
do j=1,k1
|
|
a(i,j) = 0.0_wp_
|
|
end do
|
|
end do
|
|
l = k1
|
|
do it=1,m
|
|
! fetch the current data point x(it),y(it).
|
|
xi = x(it)
|
|
wi = w(it)
|
|
yi = y(it)*wi
|
|
! search for knot interval t(l) <= xi < t(l+1).
|
|
do
|
|
if(xi<t(l+1) .or. l==nk1) exit
|
|
l = l+1
|
|
end do
|
|
! evaluate the (k+1) non-zero b-splines at xi and store them in q.
|
|
call fpbspl(t,n,k,xi,l,h)
|
|
do i=1,k1
|
|
q(it,i) = h(i)
|
|
h(i) = h(i)*wi
|
|
end do
|
|
! rotate the new row of the observation matrix into triangle.
|
|
j = l-k1
|
|
do i=1,k1
|
|
j = j+1
|
|
piv = h(i)
|
|
if(piv==0.0_wp_) cycle
|
|
! calculate the parameters of the givens transformation.
|
|
call fpgivs(piv,a(j,1),cs,sn)
|
|
! transformations to right hand side.
|
|
call fprota(cs,sn,yi,z(j))
|
|
if(i==k1) exit
|
|
i2 = 1
|
|
i3 = i+1
|
|
do i1 = i3,k1
|
|
i2 = i2+1
|
|
! transformations to left hand side.
|
|
call fprota(cs,sn,h(i1),a(j,i2))
|
|
end do
|
|
end do
|
|
! add contribution of this row to the sum of squares of residual
|
|
! right hand sides.
|
|
fp = fp+yi**2
|
|
end do
|
|
if(ier==(-2)) fp0 = fp
|
|
fpint(n) = fp0
|
|
fpint(n-1) = fpold
|
|
nrdata(n) = nplus
|
|
! backward substitution to obtain the b-spline coefficients.
|
|
call fpback(a,z,nk1,k1,c,nest)
|
|
! test whether the approximation sinf(x) is an acceptable solution.
|
|
if(iopt<0) return
|
|
fpms = fp-s
|
|
if(abs(fpms)<acc) return
|
|
! if f(p=inf) < s accept the choice of knots.
|
|
if(fpms<0.0_wp_) exit
|
|
! if n = nmax, sinf(x) is an interpolating spline.
|
|
if(n==nmax) then
|
|
ier = -1
|
|
return
|
|
end if
|
|
! increase the number of knots.
|
|
! if n=nest we cannot increase the number of knots because of
|
|
! the storage capacity limitation.
|
|
if(n==nest) then
|
|
ier = 1
|
|
return
|
|
end if
|
|
! determine the number of knots nplus we are going to add.
|
|
if(ier/=0) then
|
|
nplus = 1
|
|
ier = 0
|
|
else
|
|
npl1 = nplus*2
|
|
rn = nplus
|
|
if(fpold-fp>acc) npl1 = int(rn*fpms/(fpold-fp))
|
|
nplus = min0(nplus*2,max0(npl1,nplus/2,1))
|
|
end if
|
|
fpold = fp
|
|
! compute the sum((w(i)*(y(i)-s(x(i))))**2) for each knot interval
|
|
! t(j+k) <= x(i) <= t(j+k+1) and store it in fpint(j),j=1,2,...nrint.
|
|
fpart = 0.0_wp_
|
|
i = 1
|
|
l = k2
|
|
new = 0
|
|
do it=1,m
|
|
if(x(it)>=t(l) .and. l<=nk1) then
|
|
new = 1
|
|
l = l+1
|
|
end if
|
|
term = 0.0_wp_
|
|
l0 = l-k2
|
|
do j=1,k1
|
|
l0 = l0+1
|
|
term = term+c(l0)*q(it,j)
|
|
end do
|
|
term = (w(it)*(term-y(it)))**2
|
|
fpart = fpart+term
|
|
if(new==0) cycle
|
|
store = term*half
|
|
fpint(i) = fpart-store
|
|
i = i+1
|
|
fpart = store
|
|
new = 0
|
|
end do
|
|
fpint(nrint) = fpart
|
|
do l=1,nplus
|
|
! add a new knot.
|
|
call fpknot(x,m,t,n,fpint,nrdata,nrint,nest,1)
|
|
! if n=nmax we locate the knots as for interpolation.
|
|
if(n==nmax) then
|
|
rstart=.true.
|
|
exit
|
|
end if
|
|
! test whether we cannot further increase the number of knots.
|
|
if(n==nest) exit
|
|
end do
|
|
if(rstart) exit
|
|
! restart the computations with the new set of knots.
|
|
end do
|
|
if(rstart) then
|
|
! find the position of the interior knots in case of interpolation.
|
|
mk1 = m-k1
|
|
if (mk1/=0) then
|
|
k3 = k/2
|
|
i = k2
|
|
j = k3+2
|
|
if(k3*2/=k) then
|
|
do l=1,mk1
|
|
t(i) = x(j)
|
|
i = i+1
|
|
j = j+1
|
|
end do
|
|
else
|
|
do l=1,mk1
|
|
t(i) = (x(j)+x(j-1))*half
|
|
i = i+1
|
|
j = j+1
|
|
end do
|
|
end if
|
|
end if
|
|
else
|
|
exit
|
|
end if
|
|
end do
|
|
! test whether the least-squares kth degree polynomial is a solution
|
|
! of our approximation problem.
|
|
if(ier==(-2)) return
|
|
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
! part 2: determination of the smoothing spline sp(x). c
|
|
! *************************************************** c
|
|
! we have determined the number of knots and their position. c
|
|
! we now compute the b-spline coefficients of the smoothing spline c
|
|
! sp(x). the observation matrix a is extended by the rows of matrix c
|
|
! b expressing that the kth derivative discontinuities of sp(x) at c
|
|
! the interior knots t(k+2),...t(n-k-1) must be zero. the corres- c
|
|
! ponding weights of these additional rows are set to 1/p. c
|
|
! iteratively we then have to determine the value of p such that c
|
|
! f(p)=sum((w(i)*(y(i)-sp(x(i))))**2) be = s. we already know that c
|
|
! the least-squares kth degree polynomial corresponds to p=0, and c
|
|
! that the least-squares spline corresponds to p=infinity. the c
|
|
! iteration process which is proposed here, makes use of rational c
|
|
! interpolation. since f(p) is a convex and strictly decreasing c
|
|
! function of p, it can be approximated by a rational function c
|
|
! r(p) = (u*p+v)/(p+w). three values of p(p1,p2,p3) with correspond- c
|
|
! ing values of f(p) (f1=f(p1)-s,f2=f(p2)-s,f3=f(p3)-s) are used c
|
|
! to calculate the new value of p such that r(p)=s. convergence is c
|
|
! guaranteed by taking f1>0 and f3<0. c
|
|
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
! evaluate the discontinuity jump of the kth derivative of the
|
|
! b-splines at the knots t(l),l=k+2,...n-k-1 and store in b.
|
|
call fpdisc(t,n,k2,b,nest)
|
|
! initial value for p.
|
|
p1 = 0.0_wp_
|
|
f1 = fp0-s
|
|
p3 = -one
|
|
f3 = fpms
|
|
p = 0.0_wp_
|
|
do i=1,nk1
|
|
p = p+a(i,1)
|
|
end do
|
|
rn = nk1
|
|
p = rn/p
|
|
ich1 = 0
|
|
ich3 = 0
|
|
n8 = n-nmin
|
|
! iteration process to find the root of f(p) = s.
|
|
do iter=1,maxit
|
|
! the rows of matrix b with weight 1/p are rotated into the
|
|
! triangularised observation matrix a which is stored in g.
|
|
pinv = one/p
|
|
do i=1,nk1
|
|
c(i) = z(i)
|
|
g(i,k2) = 0.0_wp_
|
|
do j=1,k1
|
|
g(i,j) = a(i,j)
|
|
end do
|
|
end do
|
|
do it=1,n8
|
|
! the row of matrix b is rotated into triangle by givens transformation
|
|
do i=1,k2
|
|
h(i) = b(it,i)*pinv
|
|
end do
|
|
yi = 0.0_wp_
|
|
do j=it,nk1
|
|
piv = h(1)
|
|
! calculate the parameters of the givens transformation.
|
|
call fpgivs(piv,g(j,1),cs,sn)
|
|
! transformations to right hand side.
|
|
call fprota(cs,sn,yi,c(j))
|
|
if(j==nk1) exit
|
|
i2 = k1
|
|
if(j>n8) i2 = nk1-j
|
|
do i=1,i2
|
|
! transformations to left hand side.
|
|
i1 = i+1
|
|
call fprota(cs,sn,h(i1),g(j,i1))
|
|
h(i) = h(i1)
|
|
end do
|
|
h(i2+1) = 0.0_wp_
|
|
end do
|
|
end do
|
|
! backward substitution to obtain the b-spline coefficients.
|
|
call fpback(g,c,nk1,k2,c,nest)
|
|
! computation of f(p).
|
|
fp = 0.0_wp_
|
|
l = k2
|
|
do it=1,m
|
|
if(x(it)>=t(l) .and. l<=nk1) l = l+1
|
|
l0 = l-k2
|
|
term = 0.0_wp_
|
|
do j=1,k1
|
|
l0 = l0+1
|
|
term = term+c(l0)*q(it,j)
|
|
end do
|
|
fp = fp+(w(it)*(term-y(it)))**2
|
|
end do
|
|
! test whether the approximation sp(x) is an acceptable solution.
|
|
fpms = fp-s
|
|
if(abs(fpms)<acc) return
|
|
! test whether the maximal number of iterations is reached.
|
|
if(iter==maxit) then
|
|
ier = 3
|
|
return
|
|
end if
|
|
! carry out one more step of the iteration process.
|
|
p2 = p
|
|
f2 = fpms
|
|
if(ich3==0) then
|
|
if((f2-f3)<=acc) then
|
|
! our initial choice of p is too large.
|
|
p3 = p2
|
|
f3 = f2
|
|
p = p*con4
|
|
if(p<=p1) p=p1*con9 + p2*con1
|
|
cycle
|
|
end if
|
|
if(f2<0.0_wp_) ich3=1
|
|
end if
|
|
if(ich1==0) then
|
|
if((f1-f2)<=acc) then
|
|
! our initial choice of p is too small
|
|
p1 = p2
|
|
f1 = f2
|
|
p = p/con4
|
|
if(p3>=0.0_wp_ .and. p>=p3) p = p2*con1 + p3*con9
|
|
cycle
|
|
end if
|
|
if(f2>0.0_wp_) ich1=1
|
|
end if
|
|
! test whether the iteration process proceeds as theoretically
|
|
! expected.
|
|
if(f2>=f1 .or. f2<=f3) then
|
|
ier = 2
|
|
return
|
|
end if
|
|
! find the new value for p.
|
|
p = fprati(p1,f1,p2,f2,p3,f3)
|
|
end do
|
|
end subroutine fpcurf
|
|
|
|
subroutine splder(t,n,c,k,nu,x,y,m,wrk,ier)
|
|
! subroutine splder evaluates in a number of points x(i),i=1,2,...,m
|
|
! the derivative of order nu of a spline s(x) of degree k,given in
|
|
! its b-spline representation.
|
|
!
|
|
! calling sequence:
|
|
! call splder(t,n,c,k,nu,x,y,m,wrk,ier)
|
|
!
|
|
! input parameters:
|
|
! t : array,length n, which contains the position of the knots.
|
|
! n : integer, giving the total number of knots of s(x).
|
|
! c : array,length n, which contains the b-spline coefficients.
|
|
! k : integer, giving the degree of s(x).
|
|
! nu : integer, specifying the order of the derivative. 0<=nu<=k
|
|
! x : array,length m, which contains the points where the deriv-
|
|
! ative of s(x) must be evaluated.
|
|
! m : integer, giving the number of points where the derivative
|
|
! of s(x) must be evaluated
|
|
! wrk : real array of dimension n. used as working space.
|
|
!
|
|
! output parameters:
|
|
! y : array,length m, giving the value of the derivative of s(x)
|
|
! at the different points.
|
|
! ier : error flag
|
|
! ier = 0 : normal return
|
|
! ier =10 : invalid input data (see restrictions)
|
|
!
|
|
! restrictions:
|
|
! 0 <= nu <= k
|
|
! m >= 1
|
|
! t(k+1) <= x(i) <= x(i+1) <= t(n-k) , i=1,2,...,m-1.
|
|
!
|
|
! other subroutines required: fpbspl
|
|
!
|
|
! references :
|
|
! de boor c : on calculating with b-splines, j. approximation theory
|
|
! 6 (1972) 50-62.
|
|
! cox m.g. : the numerical evaluation of b-splines, j. inst. maths
|
|
! applics 10 (1972) 134-149.
|
|
! dierckx p. : curve and surface fitting with splines, monographs on
|
|
! numerical analysis, oxford university press, 1993.
|
|
!
|
|
! author :
|
|
! p.dierckx
|
|
! dept. computer science, k.u.leuven
|
|
! celestijnenlaan 200a, b-3001 heverlee, belgium.
|
|
! e-mail : Paul.Dierckx@cs.kuleuven.ac.be
|
|
!
|
|
! latest update : march 1987
|
|
!
|
|
implicit none
|
|
! arguments
|
|
integer, intent(in) :: n, k, nu, m
|
|
integer, intent(out) :: ier
|
|
real(wp_), intent(in) :: t(n), c(n), x(m)
|
|
real(wp_), intent(out) :: y(m), wrk(n)
|
|
! local variables
|
|
integer :: i, j, kk, k1, k2, l, ll, l1, l2, nk1, nk2
|
|
real(wp_) :: ak, arg, fac, sp, tb, te, h(6)
|
|
! before starting computations a data check is made. if the input data
|
|
! are invalid control is immediately repassed to the calling program.
|
|
ier = 10
|
|
if(nu<0 .or. nu>k) return
|
|
if(m<1) return
|
|
do i=2,m
|
|
if(x(i)<x(i-1)) return
|
|
end do
|
|
ier = 0
|
|
! fetch tb and te, the boundaries of the approximation interval.
|
|
k1 = k+1
|
|
nk1 = n-k1
|
|
tb = t(k1)
|
|
te = t(nk1+1)
|
|
! the derivative of order nu of a spline of degree k is a spline of
|
|
! degree k-nu,the b-spline coefficients wrk(i) of which can be found
|
|
! using the recurrence scheme of de boor.
|
|
l = 1
|
|
kk = k
|
|
! nn = n
|
|
do i=1,nk1
|
|
wrk(i) = c(i)
|
|
end do
|
|
if(nu/=0) then
|
|
nk2 = nk1
|
|
do j=1,nu
|
|
ak = kk
|
|
nk2 = nk2-1
|
|
l1 = l
|
|
do i=1,nk2
|
|
l1 = l1+1
|
|
l2 = l1+kk
|
|
fac = t(l2)-t(l1)
|
|
if(fac>0.0_wp_) wrk(i) = ak*(wrk(i+1)-wrk(i))/fac
|
|
end do
|
|
l = l+1
|
|
kk = kk-1
|
|
end do
|
|
if(kk==0) then
|
|
! if nu=k the derivative is a piecewise constant function
|
|
j = 1
|
|
do i=1,m
|
|
arg = x(i)
|
|
do
|
|
if(arg<t(l+1) .or. l==nk1) exit
|
|
l = l+1
|
|
j = j+1
|
|
end do
|
|
y(i) = wrk(j)
|
|
end do
|
|
return
|
|
end if
|
|
end if
|
|
l = k1
|
|
l1 = l+1
|
|
k2 = k1-nu
|
|
! main loop for the different points.
|
|
do i=1,m
|
|
! fetch a new x-value arg.
|
|
arg = x(i)
|
|
if(arg<tb) arg = tb
|
|
if(arg>te) arg = te
|
|
! search for knot interval t(l) <= arg < t(l+1)
|
|
do
|
|
if(arg<t(l1) .or. l==nk1) exit
|
|
l = l1
|
|
l1 = l+1
|
|
end do
|
|
! evaluate the non-zero b-splines of degree k-nu at arg.
|
|
call fpbspl(t,n,kk,arg,l,h)
|
|
! find the value of the derivative at x=arg.
|
|
sp = 0.0_wp_
|
|
ll = l-k1
|
|
do j=1,k2
|
|
ll = ll+1
|
|
sp = sp+wrk(ll)*h(j)
|
|
end do
|
|
y(i) = sp
|
|
end do
|
|
end subroutine splder
|
|
|
|
subroutine splev(t,n,c,k,x,y,m,ier)
|
|
! subroutine splev evaluates in a number of points x(i),i=1,2,...,m
|
|
! a spline s(x) of degree k, given in its b-spline representation.
|
|
!
|
|
! calling sequence:
|
|
! call splev(t,n,c,k,x,y,m,ier)
|
|
!
|
|
! input parameters:
|
|
! t : array,length n, which contains the position of the knots.
|
|
! n : integer, giving the total number of knots of s(x).
|
|
! c : array,length n, which contains the b-spline coefficients.
|
|
! k : integer, giving the degree of s(x).
|
|
! x : array,length m, which contains the points where s(x) must
|
|
! be evaluated.
|
|
! m : integer, giving the number of points where s(x) must be
|
|
! evaluated.
|
|
!
|
|
! output parameter:
|
|
! y : array,length m, giving the value of s(x) at the different
|
|
! points.
|
|
! ier : error flag
|
|
! ier = 0 : normal return
|
|
! ier =10 : invalid input data (see restrictions)
|
|
!
|
|
! restrictions:
|
|
! m >= 1
|
|
! t(k+1) <= x(i) <= x(i+1) <= t(n-k) , i=1,2,...,m-1.
|
|
!
|
|
! other subroutines required: fpbspl.
|
|
!
|
|
! references :
|
|
! de boor c : on calculating with b-splines, j. approximation theory
|
|
! 6 (1972) 50-62.
|
|
! cox m.g. : the numerical evaluation of b-splines, j. inst. maths
|
|
! applics 10 (1972) 134-149.
|
|
! dierckx p. : curve and surface fitting with splines, monographs on
|
|
! numerical analysis, oxford university press, 1993.
|
|
!
|
|
! author :
|
|
! p.dierckx
|
|
! dept. computer science, k.u.leuven
|
|
! celestijnenlaan 200a, b-3001 heverlee, belgium.
|
|
! e-mail : Paul.Dierckx@cs.kuleuven.ac.be
|
|
!
|
|
! latest update : march 1987
|
|
!
|
|
implicit none
|
|
! arguments
|
|
integer, intent(in) :: n, k, m
|
|
integer, intent(out) :: ier
|
|
real(wp_), intent(in) :: t(n), c(n), x(m)
|
|
real(wp_), intent(out) :: y(m)
|
|
! local variables
|
|
integer :: i, j, k1, l, ll, l1, nk1
|
|
real(wp_) :: arg, sp, tb, te, h(6)
|
|
! ..
|
|
! before starting computations a data check is made. if the input data
|
|
! are invalid control is immediately repassed to the calling program.
|
|
ier = 10
|
|
if(m<1) return
|
|
do i=2,m
|
|
if(x(i)<x(i-1)) return
|
|
end do
|
|
ier = 0
|
|
! 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
|
|
! main loop for the different points.
|
|
do i=1,m
|
|
! fetch a new x-value arg.
|
|
arg = x(i)
|
|
if(arg<tb) arg = tb
|
|
if(arg>te) arg = te
|
|
! search for knot interval t(l) <= arg < t(l+1)
|
|
do
|
|
if(arg<t(l1) .or. l==nk1) exit
|
|
l = l1
|
|
l1 = l+1
|
|
end do
|
|
! evaluate the non-zero b-splines at arg.
|
|
call fpbspl(t,n,k,arg,l,h)
|
|
! find the value of s(x) at x=arg.
|
|
sp = 0.0_wp_
|
|
ll = l-k1
|
|
do j=1,k1
|
|
ll = ll+1
|
|
sp = sp+c(ll)*h(j)
|
|
end do
|
|
y(i) = sp
|
|
end do
|
|
end subroutine splev
|
|
|
|
subroutine sproota(val,t,n,c,zero,mest,m,ier)
|
|
! subroutine sproot finds the zeros of a cubic spline s(x),which is
|
|
! given in its normalized b-spline representation.
|
|
!
|
|
! calling sequence:
|
|
! call sproot(t,n,c,zero,mest,m,ier)
|
|
!
|
|
! input parameters:
|
|
! t : real array,length n, containing the knots of s(x).
|
|
! n : integer, containing the number of knots. n>=8
|
|
! c : real array,length n, containing the b-spline coefficients.
|
|
! mest : integer, specifying the dimension of array zero.
|
|
!
|
|
! output parameters:
|
|
! zero : real array,lenth mest, containing the zeros of s(x).
|
|
! m : integer,giving the number of zeros.
|
|
! ier : error flag:
|
|
! ier = 0: normal return.
|
|
! ier = 1: the number of zeros exceeds mest.
|
|
! ier =10: invalid input data (see restrictions).
|
|
!
|
|
! other subroutines required: fpcuro
|
|
!
|
|
! restrictions:
|
|
! 1) n>= 8.
|
|
! 2) t(4) < t(5) < ... < t(n-4) < t(n-3).
|
|
! t(1) <= t(2) <= t(3) <= t(4)
|
|
! t(n-3) <= t(n-2) <= t(n-1) <= t(n)
|
|
!
|
|
! author :
|
|
! p.dierckx
|
|
! dept. computer science, k.u.leuven
|
|
! celestijnenlaan 200a, b-3001 heverlee, belgium.
|
|
! e-mail : Paul.Dierckx@cs.kuleuven.ac.be
|
|
!
|
|
! latest update : march 1987
|
|
!
|
|
! ..
|
|
implicit none
|
|
! ..scalar arguments..
|
|
integer, intent(in) :: n,mest
|
|
integer, intent(out) :: m,ier
|
|
real(wp_), intent(in) :: val
|
|
! ..array arguments..
|
|
real(wp_), intent(in) :: t(n),c(n)
|
|
real(wp_), intent(out) :: zero(mest)
|
|
! ..local scalars..
|
|
integer :: i,j,j1,l,n4
|
|
real(wp_) :: ah,a0,a1,a2,a3,bh,b0,b1,c1,c2,c3,c4,c5,d4,d5,h1,h2, &
|
|
t1,t2,t3,t4,t5,zz
|
|
logical :: z0,z1,z2,z3,z4,nz0,nz1,nz2,nz3,nz4
|
|
! ..local array..
|
|
real(wp_) :: y(3)
|
|
! ..
|
|
! set some constants
|
|
real(wp_), parameter :: two = 0.2e+01_wp_, three = 0.3e+01_wp_
|
|
! before starting computations a data check is made. if the input data
|
|
! are invalid, control is immediately repassed to the calling program.
|
|
n4 = n-4
|
|
ier = 10
|
|
if(n<8) return
|
|
j = n
|
|
do i=1,3
|
|
if(t(i)>t(i+1)) return
|
|
if(t(j)<t(j-1)) return
|
|
j = j-1
|
|
end do
|
|
do i=4,n4
|
|
if(t(i)>=t(i+1)) return
|
|
end do
|
|
! the problem considered reduces to finding the zeros of the cubic
|
|
! polynomials pl(x) which define the cubic spline in each knot
|
|
! interval t(l)<=x<=t(l+1). a zero of pl(x) is also a zero of s(x) on
|
|
! the condition that it belongs to the knot interval.
|
|
! the cubic polynomial pl(x) is determined by computing s(t(l)),
|
|
! s'(t(l)),s(t(l+1)) and s'(t(l+1)). in fact we only have to compute
|
|
! s(t(l+1)) and s'(t(l+1)); because of the continuity conditions of
|
|
! splines and their derivatives, the value of s(t(l)) and s'(t(l))
|
|
! is already known from the foregoing knot interval.
|
|
ier = 0
|
|
! evaluate some constants for the first knot interval
|
|
h1 = t(4)-t(3)
|
|
h2 = t(5)-t(4)
|
|
t1 = t(4)-t(2)
|
|
t2 = t(5)-t(3)
|
|
t3 = t(6)-t(4)
|
|
t4 = t(5)-t(2)
|
|
t5 = t(6)-t(3)
|
|
! calculate a0 = s(t(4)) and ah = s'(t(4)).
|
|
c1 = c(1)
|
|
c2 = c(2)
|
|
c3 = c(3)
|
|
c4 = (c2-c1)/t4
|
|
c5 = (c3-c2)/t5
|
|
d4 = (h2*c1+t1*c2)/t4
|
|
d5 = (t3*c2+h1*c3)/t5
|
|
a0 = (h2*d4+h1*d5)/t2 - val
|
|
ah = three*(h2*c4+h1*c5)/t2
|
|
z1 = .true.
|
|
if(ah<0.0_wp_) z1 = .false.
|
|
nz1 = .not.z1
|
|
m = 0
|
|
! main loop for the different knot intervals.
|
|
do l=4,n4
|
|
! evaluate some constants for the knot interval t(l) <= x <= t(l+1).
|
|
h1 = h2
|
|
h2 = t(l+2)-t(l+1)
|
|
t1 = t2
|
|
t2 = t3
|
|
t3 = t(l+3)-t(l+1)
|
|
t4 = t5
|
|
t5 = t(l+3)-t(l)
|
|
! find a0 = s(t(l)), ah = s'(t(l)), b0 = s(t(l+1)) and bh = s'(t(l+1)).
|
|
c1 = c2
|
|
c2 = c3
|
|
c3 = c(l)
|
|
c4 = c5
|
|
c5 = (c3-c2)/t5
|
|
d4 = (h2*c1+t1*c2)/t4
|
|
d5 = (h1*c3+t3*c2)/t5
|
|
b0 = (h2*d4+h1*d5)/t2 - val
|
|
bh = three*(h2*c4+h1*c5)/t2
|
|
! calculate the coefficients a0,a1,a2 and a3 of the cubic polynomial
|
|
! pl(x) = ql(y) = a0+a1*y+a2*y**2+a3*y**3 ; y = (x-t(l))/(t(l+1)-t(l)).
|
|
a1 = ah*h1
|
|
b1 = bh*h1
|
|
a2 = three*(b0-a0)-b1-two*a1
|
|
a3 = two*(a0-b0)+b1+a1
|
|
! test whether or not pl(x) could have a zero in the range
|
|
! t(l) <= x <= t(l+1).
|
|
z3 = .true.
|
|
if(b1<0.0_wp_) z3 = .false.
|
|
nz3 = .not.z3
|
|
if(a0*b0>0.0_wp_) then
|
|
z0 = .true.
|
|
if(a0<0.0_wp_) z0 = .false.
|
|
nz0 = .not.z0
|
|
z2 = .true.
|
|
if(a2<0.0_wp_) z2 = .false.
|
|
nz2 = .not.z2
|
|
z4 = .true.
|
|
if(3.0_wp_*a3+a2<0.0_wp_) z4 = .false.
|
|
nz4 = .not.z4
|
|
else
|
|
z0 = .true.
|
|
nz0 = .not.z0
|
|
z2 = .true.
|
|
nz2 = .not.z2
|
|
z4 = .true.
|
|
nz4 = .not.z4
|
|
end if
|
|
if(( z0.and.(nz1.and.( z3.or. z2.and.nz4).or.nz2.and. z3.and. z4) &
|
|
.or.nz0.and.( z1.and.(nz3.or.nz2.and. z4).or. z2.and.nz3.and.nz4) &
|
|
) .or. (a0*b0<=0.0_wp_) ) then
|
|
! find the zeros of ql(y).
|
|
call fpcuro(a3,a2,a1,a0,y,j)
|
|
! find which zeros of pl(x) are zeros of s(x).
|
|
do i=1,j
|
|
if(y(i)<0.0_wp_ .or. y(i)>1.0_wp_) cycle
|
|
! test whether the number of zeros of s(x) exceeds mest.
|
|
if(m>=mest) then
|
|
ier = 1
|
|
return
|
|
end if
|
|
m = m+1
|
|
zero(m) = t(l)+h1*y(i)
|
|
end do
|
|
end if
|
|
a0 = b0
|
|
ah = bh
|
|
z1 = z3
|
|
nz1 = nz3
|
|
end do
|
|
! the zeros of s(x) are arranged in increasing order.
|
|
do i=2,m
|
|
j = i
|
|
do
|
|
j1 = j-1
|
|
if(j1==0) exit
|
|
if(zero(j)>=zero(j1)) exit
|
|
zz = zero(j)
|
|
zero(j) = zero(j1)
|
|
zero(j1) = zz
|
|
j = j1
|
|
end do
|
|
end do
|
|
j = m
|
|
m = 1
|
|
do i=2,j
|
|
if(zero(i)==zero(m)) cycle
|
|
m = m+1
|
|
zero(m) = zero(i)
|
|
end do
|
|
end subroutine sproota
|
|
|
|
subroutine profil(iopt,tx,nx,ty,ny,c,kx,ky,u,nu,cu,ier)
|
|
! if iopt=0 subroutine profil calculates the b-spline coefficients of
|
|
! the univariate spline f(y) = s(u,y) with s(x,y) a bivariate spline of
|
|
! degrees kx and ky, given in the b-spline representation.
|
|
! if iopt = 1 it calculates the b-spline coefficients of the univariate
|
|
! spline g(x) = s(x,u)
|
|
!
|
|
! calling sequence:
|
|
! call profil(iopt,tx,nx,ty,ny,c,kx,ky,u,nu,cu,ier)
|
|
!
|
|
! input parameters:
|
|
! iopt : integer flag, specifying whether the profile f(y) (iopt=0)
|
|
! or the profile g(x) (iopt=1) must be determined.
|
|
! tx : real array, length nx, which contains the position of the
|
|
! knots in the x-direction.
|
|
! nx : integer, giving the total number of knots in the x-direction
|
|
! ty : real array, length ny, which contains the position of the
|
|
! knots in the y-direction.
|
|
! ny : integer, giving the total number of knots in the y-direction
|
|
! c : real array, length (nx-kx-1)*(ny-ky-1), which contains the
|
|
! b-spline coefficients.
|
|
! kx,ky : integer values, giving the degrees of the spline.
|
|
! u : real value, specifying the requested profile.
|
|
! tx(kx+1)<=u<=tx(nx-kx), if iopt=0.
|
|
! ty(ky+1)<=u<=ty(ny-ky), if iopt=1.
|
|
! nu : on entry nu must specify the dimension of the array cu.
|
|
! nu >= ny if iopt=0, nu >= nx if iopt=1.
|
|
!
|
|
! output parameters:
|
|
! cu : real array of dimension (nu).
|
|
! on succesful exit this array contains the b-spline
|
|
! ier : integer error flag
|
|
! ier=0 : normal return
|
|
! ier=10: invalid input data (see restrictions)
|
|
!
|
|
! restrictions:
|
|
! if iopt=0 : tx(kx+1) <= u <= tx(nx-kx), nu >=ny.
|
|
! if iopt=1 : ty(ky+1) <= u <= ty(ny-ky), nu >=nx.
|
|
!
|
|
! other subroutines required:
|
|
! fpbspl
|
|
!
|
|
! author :
|
|
! p.dierckx
|
|
! dept. computer science, k.u.leuven
|
|
! celestijnenlaan 200a, b-3001 heverlee, belgium.
|
|
! e-mail : Paul.Dierckx@cs.kuleuven.ac.be
|
|
!
|
|
! latest update : march 1987
|
|
!
|
|
implicit none
|
|
! ..scalar arguments..
|
|
integer,intent(in) :: iopt,nx,ny,kx,ky,nu
|
|
integer,intent(out) :: ier
|
|
real(wp_), intent(in) :: u
|
|
! ..array arguments..
|
|
real(wp_),intent(in) :: tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1))
|
|
real(wp_),intent(out) :: cu(nu)
|
|
! ..local scalars..
|
|
integer :: i,j,kx1,ky1,l,l1,m,m0,nkx1,nky1
|
|
real(wp_) :: summ
|
|
! ..local array
|
|
real(wp_) :: h(6)
|
|
! ..
|
|
! before starting computations a data check is made. if the input data
|
|
! are invalid control is immediately repassed to the calling program.
|
|
kx1 = kx+1
|
|
ky1 = ky+1
|
|
nkx1 = nx-kx1
|
|
nky1 = ny-ky1
|
|
ier = 10
|
|
if(iopt==0) then
|
|
if(nu<ny) return
|
|
if(u<tx(kx1) .or. u>tx(nkx1+1)) return
|
|
! the b-splinecoefficients of f(y) = s(u,y).
|
|
ier = 0
|
|
l = kx1
|
|
l1 = l+1
|
|
do
|
|
if(u<tx(l1) .or. l==nkx1) exit
|
|
l = l1
|
|
l1 = l+1
|
|
end do
|
|
call fpbspl(tx,nx,kx,u,l,h)
|
|
m0 = (l-kx1)*nky1+1
|
|
do i=1,nky1
|
|
m = m0
|
|
summ = 0.0_wp_
|
|
do j=1,kx1
|
|
summ = summ+h(j)*c(m)
|
|
m = m+nky1
|
|
end do
|
|
cu(i) = summ
|
|
m0 = m0+1
|
|
end do
|
|
else
|
|
if(nu<nx) return
|
|
if(u<ty(ky1) .or. u>ty(nky1+1)) return
|
|
! the b-splinecoefficients of g(x) = s(x,u).
|
|
ier = 0
|
|
l = ky1
|
|
l1 = l+1
|
|
do
|
|
if(u<ty(l1) .or. l==nky1) exit
|
|
l = l1
|
|
l1 = l+1
|
|
end do
|
|
call fpbspl(ty,ny,ky,u,l,h)
|
|
m0 = l-ky
|
|
do i=1,nkx1
|
|
m = m0
|
|
summ = 0.0_wp_
|
|
do j=1,ky1
|
|
summ = summ+h(j)*c(m)
|
|
m = m+1
|
|
end do
|
|
cu(i) = summ
|
|
m0 = m0+nky1
|
|
end do
|
|
end if
|
|
end subroutine profil
|
|
|
|
subroutine fpcuro(a,b,c,d,x,n)
|
|
! subroutine fpcuro finds the real zeros of a cubic polynomial
|
|
! p(x) = a*x**3+b*x**2+c*x+d.
|
|
!
|
|
! calling sequence:
|
|
! call fpcuro(a,b,c,d,x,n)
|
|
!
|
|
! input parameters:
|
|
! a,b,c,d: real values, containing the coefficients of p(x).
|
|
!
|
|
! output parameters:
|
|
! x : real array,length 3, which contains the real zeros of p(x)
|
|
! n : integer, giving the number of real zeros of p(x).
|
|
! ..
|
|
implicit none
|
|
! ..scalar arguments..
|
|
real(wp_), intent(in) :: a,b,c,d
|
|
integer, intent(out) :: n
|
|
! ..array argument..
|
|
real(wp_), intent(out) :: x(3)
|
|
! ..local scalars..
|
|
integer :: i
|
|
real(wp_) :: a1,b1,c1,df,disc,d1,f,pi3,p3,q,r,step,u,u1,u2,y
|
|
! ..function references..
|
|
! real(8) abs,max,atan,atan2,cos,sign,sqrt
|
|
! set constants
|
|
real(wp_), parameter :: two =0.2e+01_wp_,three=0.3e+01_wp_,four=0.4e+01_wp_
|
|
real(wp_), parameter :: ovfl=0.1e+05_wp_,half =0.5e+00_wp_,tent=0.1e+00_wp_
|
|
real(wp_), parameter :: e3 = tent/0.3_wp_
|
|
pi3 = atan(0.1e+01_wp_)/0.75_wp_
|
|
a1 = abs(a)
|
|
b1 = abs(b)
|
|
c1 = abs(c)
|
|
d1 = abs(d)
|
|
! test whether p(x) is a third degree polynomial.
|
|
if(max(b1,c1,d1)>=a1*ovfl) then
|
|
! test whether p(x) is a second degree polynomial.
|
|
if(max(c1,d1)>=b1*ovfl) then
|
|
! test whether p(x) is a first degree polynomial.
|
|
if(d1>=c1*ovfl) then
|
|
! p(x) is a constant function.
|
|
n = 0
|
|
return
|
|
end if
|
|
! p(x) is a first degree polynomial.
|
|
n = 1
|
|
x(1) = -d/c
|
|
else
|
|
! p(x) is a second degree polynomial.
|
|
disc = c*c-four*b*d
|
|
n = 0
|
|
if(disc<0.0_wp_) return
|
|
n = 2
|
|
u = sqrt(disc)
|
|
b1 = b+b
|
|
x(1) = (-c+u)/b1
|
|
x(2) = (-c-u)/b1
|
|
end if
|
|
else
|
|
! p(x) is a third degree polynomial.
|
|
b1 = b/a*e3
|
|
c1 = c/a
|
|
d1 = d/a
|
|
q = c1*e3-b1*b1
|
|
r = b1*b1*b1+(d1-b1*c1)*half
|
|
disc = q*q*q+r*r
|
|
if(disc<=0.0_wp_) then
|
|
u = sqrt(abs(q))
|
|
if(r<0.0_wp_) u = -u
|
|
p3 = atan2(sqrt(-disc),abs(r))*e3
|
|
u2 = u+u
|
|
n = 3
|
|
x(1) = -u2*cos(p3)-b1
|
|
x(2) = u2*cos(pi3-p3)-b1
|
|
x(3) = u2*cos(pi3+p3)-b1
|
|
else
|
|
u = sqrt(disc)
|
|
u1 = -r+u
|
|
u2 = -r-u
|
|
n = 1
|
|
x(1) = sign(abs(u1)**e3,u1)+sign(abs(u2)**e3,u2)-b1
|
|
end if
|
|
end if
|
|
! apply a newton iteration to improve the accuracy of the roots.
|
|
do i=1,n
|
|
y = x(i)
|
|
f = ((a*y+b)*y+c)*y+d
|
|
df = (three*a*y+two*b)*y+c
|
|
step = 0.0_wp_
|
|
if(abs(f)<abs(df)*tent) step = f/df
|
|
x(i) = y-step
|
|
end do
|
|
end subroutine fpcuro
|
|
|
|
end module dierckx
|