Homec4science
Diffusion SPClibs (feature/fft-split) Browse

spclibs/pppack/splopt.f90
ff4fbc3f5f6cfeature/fft-split

Last Change
Raw File

splopt.f90
View Options

!>
!> @file splopt.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
subroutine splopt ( tau, n, k, scrtch, t, iflag )
!*************************************************************************
!
!! SPLOPT computes the knots for an optimal recovery scheme.
!
! Discussion:
!
! The optimal recovery scheme is of order K for data at TAU(1:N).
!
! The interior knots T(K+1:N) are determined by Newton's method in
! such a way that the signum function which changes sign at
! T(K+1), ..., T(N) and nowhere else in ( TAU(1), TAU(n) ) is
! orthogonal to the spline space SPLINE ( K, TAU ) on that interval.
!
! Let XI(J) be the current guess for T(K+J), j=1,...,n-k. Then
! the next Newton iterate is of the form
! xi(j) + (-)**(n-k-j)*x(j) , j=1,...,n-k,
! with X the solution of the linear system
! C * X = D.
!
! Here, c(i,j)=b(i)(xi(j)), all j, with b(i) the i-th b-spline of
! order K for the knot sequence TAU, all i, and D is the vector
! given by d(i)=sum( -a(j) , j=i,...,n )*(TAU(i+k)-TAU(i))/k, all i,
! with a(i)=sum ( (-)**(n-k-j)*b(i,k+1,tau)(xi(j)) , j=1,...,n-k )
! for i=1,...,n-1, and a(n)=-.5 .
!
! See chapter XIII of text and references there for a derivation.
!
! The first guess for t(k+j) is (TAU(j+1)+...+TAU(j+k-1))/(k-1) .
! iteration terminates if max(abs(x(j))) < t o l , with
! TOL = t o l r t e *(TAU(n)-TAU(1))/(n-k) ,
! or else after NEWTMX iterations , currently,
! newtmx, tolrte / 10, .000001
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! Input, real ( kind = 8 ) TAU(N), the interpolation points.
! assumed to be nondecreasing, with tau(i) < tau(i+k),all i.
!
! Input, integer N, the number of data points.
!
! Input, integer K, the order of the optimal recovery scheme to be used.
!
! Workspace, real ( kind = 8 ) SCRTCH((N-K)*(2*K+3)+5*K+3). The various
! contents are specified in the text below .
!
! Output, real ( kind = 8 ) T(N+K), the optimal knots ready for
! use in optimal recovery. specifically, t(1)=... = t(k) =
! tau(1) and t(n+1)=... = t(n+k) = tau(n) , while the n-k
! interior knots t(k+1), ..., t(n) are calculated.
!
! Output, integer IFLAG, error indicator.
! = 1, success. T contains the optimal knots.
! = 2, failure. K < 3 or N < K or the linear system was singular.
!
implicit none
integer k
integer n
real ( kind = 8 ) del
real ( kind = 8 ) delmax
real ( kind = 8 ) floatk
integer i
integer id
integer iflag
integer index
integer j
integer kp1
integer kpkm1
integer kpn
integer l
integer left
integer leftmk
integer lenw
integer ll
integer llmax
integer llmin
integer na
integer nb
integer nc
integer nd
integer, parameter :: newtmx = 10
integer newton
integer nmk
integer nx
real ( kind = 8 ) scrtch((n-k)*(2*k+3)+5*k+3)
real ( kind = 8 ) t(n+k)
real ( kind = 8 ) tau(n)
real ( kind = 8 ) sign
real ( kind = 8 ) signst
real ( kind = 8 ) sum1
real ( kind = 8 ) tol
real ( kind = 8 ), parameter :: tolrte = 0.000001D+00
real ( kind = 8 ) xij
nmk = n - k
if ( n < k ) then
write ( *, * ) ' '
write ( *, * ) 'SPLOPT - Fatal error!'
write ( *, * ) ' N < K, N = ',n,' K = ',k
iflag = 2
return
end if
if ( n == k ) then
do i = 1, k
t(i) = tau(1)
t(n+i) = tau(n)
end do
return
end if
if ( k <= 2 ) then
write(*,*)' '
write(*,*)'SPLOPT - Fatal error!'
write(*,*)' K < 2, K = ',k
iflag = 2
stop
end if
floatk = k
kp1 = k+1
kpkm1 = k+k-1
kpn = k+n
signst = -1.0D+00
if ( (nmk/2) * 2 < nmk ) then
signst = 1.0D+00
end if
!
! scrtch(i)=tau-extended(i), i=1,...,n+k+k
!
nx = n + k + k + 1
!
! scrtch(i+nx)=xi(i),i=0,...,n-k+1
!
na = nx + nmk + 1
!
! scrtch(i+na)=-a(i), i=1,...,n
!
nd = na + n
!
! scrtch(i+nd)=x(i) or d(i), i=1,...,n-k
!
nb = nd+nmk
!
! scrtch(i+nb)=biatx(i),i=1,...,k+1
!
nc = nb+kp1
!
! scrtch(i+(j-1)*(2k-1)+nc)=w(i,j) = c(i-k+j,j), i=j-k,...,j+k,
! j=1,...,n-k.
!
lenw = kpkm1*nmk
!
! Extend TAU to a knot sequence and store in scrtch.
!
do j = 1, k
scrtch(j) = tau(1)
scrtch(kpn+j) = tau(n)
end do
do j = 1, n
scrtch(k+j) = tau(j)
end do
!
! First guess for scrtch (.+nx) = xi .
!
scrtch(nx) = tau(1)
scrtch(nmk+1+nx) = tau(n)
do j = 1, nmk
sum1 = 0.0D+00
do l = 1, k-1
sum1 = sum1 + tau(j+l)
end do
scrtch(j+nx) = sum1 / real ( k - 1, kind = 8 )
end do
!
! last entry of scrtch (.+na) =-a is always ...
!
scrtch(n+na) = 0.5D+00
!
! Start the Newton iteration.
!
newton = 1
tol = tolrte * ( tau(n) - tau(1) ) / real ( nmk, kind = 8 )
!
! Start the Newton step.
! compute the 2k-1 bands of the matrix c and store in scrtch(.+nc),
! and compute the vector scrtch(.+na)=-a.
!
100 continue
do i = 1, lenw
scrtch(i+nc) = 0.0D+00
end do
do i = 2, n
scrtch(i-1+na) = 0.0D+00
end do
sign = signst
left = kp1
do j = 1, nmk
xij = scrtch(j+nx)
130 continue
if ( xij < scrtch(left+1) ) then
go to 140
end if
left = left+1
if ( left < kpn ) then
go to 130
end if
left = left-1
140 continue
call bsplvb(scrtch,k,1,xij,left,scrtch(1+nb))
!
! The TAU sequence in scrtch is preceded by k additional knots
! therefore, scrtch(ll+nb) now contains b(left-2k+ll)(xij)
! which is destined for c(left-2k+ll,j), and therefore for
! w(left-k-j+ll,j)= scrtch(left-k-j+ll+(j-1)*kpkm1 + nc)
! since we store the 2k-1 bands of c in the 2k-1 r o w s of
! the work array w, and w in turn is stored in s c r t c h ,
! with w(1,1)=scrtch(1+nc).
!
! also, c being of order n-k, we would want
! 1 <= left-2k+ll .le. n-k or
! llmin=2k-left <= ll .le. n-left+k = llmax .
!
leftmk = left-k
index = leftmk-j+(j-1)*kpkm1+nc
llmin = max(1,k-leftmk)
llmax = min(k,n-leftmk)
do ll = llmin, llmax
scrtch(ll+index)=scrtch(ll+nb)
end do
call bsplvb (scrtch,kp1,2,xij,left,scrtch(1+nb))
id=max(0,leftmk-kp1)
llmin=1-min(0,leftmk-kp1)
do ll=llmin, kp1
id=id+1
scrtch(id+na)=scrtch(id+na)-sign*scrtch(ll+nb)
end do
sign=-sign
end do
call banfac(scrtch(1+nc),kpkm1,nmk,k-1,k-1,iflag)
if ( iflag == 2 ) then
write ( *, * ) ' '
write ( *, * ) 'SPLOPT - Fatal error!'
write ( *, * ) ' Matrix C is not invertible.'
stop
end if
!
! compute scrtch (.+nd)= d from scrtch (.+na) =-a .
!
do i=n,2,-1
scrtch(i-1+na)=scrtch(i-1+na)+scrtch(i+na)
end do
do i=1,nmk
scrtch(i+nd)=scrtch(i+na)*(tau(i+k)-tau(i))/floatk
end do
!
! Compute scrtch (.+nd)= x .
!
call banslv(scrtch(1+nc),kpkm1,nmk,k-1,k-1,scrtch(1+nd))
!
! Compute scrtch (.+nd)=change in xi . modify, if necessary, to
! prevent new xi from moving more than 1/3 of the way to its
! neighbors. then add to xi to obtain new xi in scrtch(.+nx).
!
delmax = 0.0D+00
sign = signst
do i = 1, nmk
del = sign * scrtch(i+nd)
delmax = max ( delmax, abs ( del ) )
if ( 0.0D+00 < del ) then
go to 230
end if
del = max ( del, ( scrtch(i-1+nx) - scrtch(i+nx) ) / 3.0D+00 )
go to 240
230 del = min (del,(scrtch(i+1+nx)-scrtch(i+nx))/3.0D+00 )
240 sign = -sign
scrtch(i+nx) = scrtch(i+nx)+del
end do
!
! Call it a day in case change in xi was small enough or too many
! steps were taken.
!
if ( delmax < tol ) then
go to 270
end if
newton = newton + 1
if ( newton <= newtmx ) then
go to 100
end if
write ( *, * ) ' '
write ( *, * ) 'SPLOPT - Warning!'
write ( *, * ) ' No convergence. Number of Newton steps was ', newtmx
270 continue
do i = 1, nmk
t(k+i) = scrtch(i+nx)
end do
290 continue
do i=1,k
t(i)=tau(1)
t(n+i)=tau(n)
end do
return
! 310 iflag=2
!
! return
end
c4science · Help