!>
!> @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