Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/pppack/smooth.f90
03540012fedefeature/refactoring

Last Change
Raw File

smooth.f90
View Options

!>
!> @file smooth.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/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
subroutine smooth ( x, y, dy, npoint, s, v, a, sfp )
!*************************************************************************
!
!! SMOOTH constructs the cubic smoothing spline to given data.
!
! Discussion:
!
! The data is of the form
!
! (x(i),y(i)), i=1,...,npoint,
!
! The cubic smoothing spline has as small a second derivative as
! possible while
!
! s(f)=sum( ((y(i)-f(x(i)))/dy(i))**2 , i=1,...,npoint ) <= s .
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! input
!
! Input, real ( kind = 8 ) X(NPOINT), the abscissas, assumed to be strictly
! increasing .
!
! Input, real ( kind = 8 ) Y(NPOINT), the corresponding ordinates.
!
! dy(1),...,dy(npoint) estimate of uncertainty in data, assumed
! to be positive .
!
! npoint.....number of data points, assumed greater than 1
!
! s.....upper bound on the discrete weighted mean square distance of
! the approximation f from the data .
!
! work arrays:
!
! v of size (npoint,7)
! a of size (npoint,4)
!
! output
!
! a(.,1).....contains the sequence of smoothed ordinates .
! a(i,j)=d**(j-1)f(x(i)), j=2,3,4, i=1,...,npoint-1 , i.e., the
! first three derivatives of the smoothing spline f at the
! left end of each of the data intervals .
! Warning . . . a would have to be transposed before it
! could be used in ppvalu .
!
! Method:
!
! the matrices q-transp*d and q-transp*d**2*q are constructed in
! SETUPQ from x and dy , as is the vector qty=q-transp*y .
! then, for given p , the vector U is determined in CHOL1D as
! the solution of the linear system
! (6(1-p)q-transp*d**2*q+p*r)u =qty .
! from u , the smoothing spline f (for this choice of smoothing par-
! ameter p ) is obtained in the sense that
! f(x(.)) = y-6(1-p)d**2*q*u and
! (d**2)f(x(.)) = 6*p*u .
!
! the smoothing parameter p is found (if possible) so that
! sf(p) = s ,
! with sf(p)=s(f) , where f is the smoothing spline as it depends
! on p . if s=0, then p = 1 . if sf(0) <= s , then p = 0 .
! otherwise, the secant method is used to locate an appropriate p in
! the open interval (0,1) . specifically,
! p(0)=0, p(1) = (s-sf(0))/dsf
! with dsf=-24*u-transp*r*u a good approximation to d(sf(0)) = dsf
! +60*(d*q*u)-transp*(d*q*u) , and u as obtained for p=0 .
! after that, for n=1,2,... until sf(p(n)) <= 1.01*s, do....
! determine p(n+1) as the point at which the secant to sf at the
! points p(n) and p(n-1) takes on the value s .
! if p(n+1) >= 1 , choose instead p(n+1) as the point at which
! the parabola sf(p(n))*((1-.)/(1-p(n)))**2 takes on the value s.
!
! Note that, in exact arithmetic, always p(n+1) < p(n) , hence
! sf(p(n+1)) < sf(p(n)) . therefore, also stop the iteration,
! with final p=1 , in case sf(p(n+1)) >= sf(p(n)) .
!
implicit none
integer npoint
real ( kind = 8 ) a(npoint,4)
real ( kind = 8 ) change
real ( kind = 8 ) dy(npoint)
integer i
real ( kind = 8 ) p
real ( kind = 8 ) prevp
real ( kind = 8 ) prevsf
real ( kind = 8 ) s
real ( kind = 8 ) sfp
real ( kind = 8 ) utru
real ( kind = 8 ) v(npoint,7)
real ( kind = 8 ) x(npoint)
real ( kind = 8 ) y(npoint)
call setupq ( x, dy, y, npoint, v, a(1,4) )
if ( 0.0D+00 < s ) then
go to 20
end if
10 continue
p = 1.0D+00
call chol1d ( p, v, a(1,4), npoint, a(1,3), a(1,1) )
sfp = 0.0D+00
go to 70
20 continue
p = 0.0D+00
call chol1d ( p, v, a(1,4), npoint, a(1,3), a(1,1) )
sfp = 0.0D+00
do i = 1, npoint
sfp = sfp + ( a(i,1) * dy(i) )**2
end do
sfp = sfp * 36.0D+00
if ( sfp <= s ) then
go to 70
end if
prevp = 0.0D+00
prevsf = sfp
utru = 0.0D+00
do i = 2, npoint
utru = utru + v(i-1,4) * ( a(i-1,3) * ( a(i-1,3) + a(i,3) ) + a(i,3)**2 )
end do
p = ( sfp - s ) / ( 24.0D+00 * utru )
!
! Secant iteration for the determination of p starts here.
!
50 continue
call chol1d ( p, v, a(1,4), npoint, a(1,3), a(1,1) )
sfp = 0.0D+00
do i = 1, npoint
sfp = sfp + ( a(i,1) * dy(i) )**2
end do
sfp = sfp * 36.0D+00 * ( 1.0D+00 - p )**2
if ( sfp <= 1.01D+00 * s ) then
go to 70
end if
if ( prevsf <= sfp ) then
go to 10
end if
change = ( p - prevp ) / ( sfp - prevsf ) * ( sfp - s )
prevp = p
p = p - change
prevsf = sfp
if ( 1.0D+00 <= p ) then
p = 1.0D+00 - sqrt ( s / prevsf ) * ( 1.0D+00 - prevp )
end if
go to 50
!
! The correct value of p has been found.
! Compute polynomial coefficients from q*u (in a(.,1)).
!
70 continue
do i = 1, npoint
a(i,1) = y(i) - 6.0D+00 * ( 1.0D+00 - p ) * dy(i)**2 * a(i,1)
end do
do i = 1, npoint
a(i,3) = 6.0D+00 * a(i,3) * p
end do
do i = 1, npoint-1
a(i,4) = ( a(i+1,3) - a(i,3) ) / v(i,4)
a(i,2) = ( a(i+1,1) - a(i,1) ) / v(i,4) &
- ( a(i,3) + a(i,4) / 3.0D+00 * v(i,4) ) / 2.0D+00 * v(i,4)
end do
return
end
c4science · Help