!>
!> @file splint.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 splint ( tau, gtau, t, n, k, q, bcoef, iflag )

!*************************************************************************
!
!! SPLINT produces the B-spline coefficients BCOEF of an interpolating spline.
!
!  Discussion:
!
!    The spline is of order K with knots T(1:N+K), and takes on the
!    value GTAU(I) at TAU(I), for I = 1 to N.
!
!    The I-th equation of the linear system
!      A * BCOEF = B
!    for the b-coefficients of the interpolant enforces interpolation
!    at TAU(1:N).
!
!    Hence,  b(i)=gtau(i), all i, and  a  is a band matrix with  2k-1
!    bands (if it is invertible).
!
!    The matrix A is generated row by row and stored, diagonal by di-
!    agonal, in the rows of the array  q , with the main diagonal go-
!    ing into row K.  see comments in the program below.
!
!    The banded system is then solved by a call to  banfac (which con-
!    structs the triangular factorization for  a  and stores it again in
!    q ), followed by a call to  banslv (which then obtains the solution
!    bcoef  by substitution).
!
!    BANFAC does no pivoting, since the total positivity of the matrix
!    A makes this unnecessary.
!
!  Reference:
!
!    Carl DeBoor,
!    A Practical Guide to Splines,
!    Springer Verlag.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) TAU(N), the data point abscissas.  The entries in
!    TAU should be strictly increasing.
!
!    Input, real ( kind = 8 ) GTAU(N), the data ordinates.
!
!    Input, real ( kind = 8 ) T(N+K), the knot sequence.
!
!    Input, integer N, the number of data points.
!
!    Input, integer K, the order of the spline.
!
!  output
!
!  q, array of size  (2*k-1)*n , containing the triangular factoriz-
!        ation of the coefficient matrix of the linear system for the b-
!        coefficients of the spline interpolant.
!        the b-coefficients for the interpolant of an additional data set can
!        be obtained without going through all the calculations in this
!        routine, simply by loading  htau  into  bcoef  and then execut-
!        ing the    call banslv ( q, 2*k-1, n, k-1, k-1, bcoef )
!
!  bcoef, the b-coefficients of the interpolant, of length  n.
!
!  iflag, an integer indicating success (= 1)  or failure (= 2)
!        the linear system to be solved is (theoretically) invertible if
!        and only if
!              t(i) < tau(i) < tau(i+k),    all i.
!        violation of this condition is certain to lead to  iflag=2 .
!
  implicit none

  integer n

  real ( kind = 8 ) bcoef(n)
  real ( kind = 8 ) gtau(n)
  integer i
  integer iflag
  integer ilp1mx
  integer j
  integer jj
  integer k
  integer kpkm2
  integer left
  real ( kind = 8 ) q((2*k-1)*n)
  real ( kind = 8 ) t(n+k)
  real ( kind = 8 ) tau(n)
  real ( kind = 8 ) taui

  kpkm2 = 2*(k-1)
  left = k

  do i = 1, (2*k-1)*n
    q(i) = 0.0D+00
  end do
!
!  loop over i to construct the  n  interpolation equations
!
  do i = 1, n

    taui = tau(i)
    ilp1mx = min(i+k,n+1)
!
!  find  left  in the closed interval (i,i+k-1) such that
!    t(left) <= tau(i)  < t(left+1)
!  matrix is singular if this is not possible
!
    left = max(left,i)

    if ( taui < t(left)) then
      go to 70
    end if

   20   continue

    if ( taui < t(left+1)) then
      go to 30
    end if

    left = left+1
    if ( left < ilp1mx) then
      go to 20
    end if

    left = left-1
    if ( t(left+1) < taui ) then
      go to 70
    end if
!
!  The i-th equation enforces interpolation at taui, hence
!  a(i,j)=b(j,k,t)(taui), all j. only the  k  entries with  j =
!  left-k+1,...,left actually might be nonzero. these  k  numbers
!  are returned, in  bcoef (used for temp.storage here), by the
!  following
!
   30   continue

    call bsplvb(t,k,1,taui,left,bcoef)
!
!  We therefore want  bcoef(j)=b(left-k+j)(taui) to go into
!  a(i,left-k+j), i.e., into  q(i-(left+j)+2*k,(left+j)-k) since
!  a(i+j,j)  is to go into  q(i+k,j), all i,j,  if we consider  q
!  as a two-dim. array , with  2*k-1  rows (see comments in
!  banfac). in the present program, we treat  q  as an equivalent
!  one-dimensional array (because of fortran restrictions on
!  dimension statements) . we therefore want  bcoef(j) to go into
!  entry
!    i -(left+j)+2*k + ((left+j)-k-1)*(2*k-1)
!   = i-left+1+(left -k)*(2*k-1) + (2*k-2)*j
!  of  q .
!
    jj = i-left+1+(left-k)*(k+k-1)
    do j = 1,k
      jj = jj+kpkm2
      q(jj) = bcoef(j)
    end do

  end do
!
!  Obtain factorization of A, stored again in Q.
!
  call banfac ( q, k+k-1, n, k-1, k-1, iflag )

  if ( iflag == 2 ) then
    write(*,*)' '
    write(*,*)'SPLINT - Fatal Error!'
    write(*,*)'  The linear system is not invertible!'
    return
  end if
!
!  Solve  a*bcoef=gtau  by backsubstitution
!
  bcoef(1:n) = gtau(1:n)

  call banslv ( q, k+k-1, n, k-1, k-1, bcoef )
  return

   70 iflag=2

  write ( *, * ) ' '
  write ( *, * ) 'SPLINT - Fatal Error!'
  write ( *, * ) '  The linear system is not invertible!'

  return
end