!>
!> @file l2appr.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 l2appr ( t, n, k, q, diag, bcoef )

!*************************************************************************
!
!! L2APPR constructs a weighted L2 spline approximation to given data.
!
!  Discussion:
!
!    The routine constructs the weighted discrete L2-approximation by
!    splines of order K with knot sequence T(1), ..., T(n+k) to
!    given data points ( TAU(1:NTAU), GTAU(1:NTAU) ).
!
!    The B-spline coefficients BCOEF of the approximating spline are
!    determined from the normal equations using Cholesky's method.
!
!  Method:
!
!    The B-spline coefficients of the L2-approximation are determined as the
!    solution of the normal equations
!
!      sum ( (b(i), b(j) ) * bcoef(j) : j=1,...,n) =(b(i),g),
!                                               i=1, ..., n .
!    Here, b(i) denotes the i-th B-spline, G denotes the function to
!    be approximated, and the inner product of two functions F and G
!    is given by
!
!      (f,g)  := sum ( f(tau(i))*g(tau(i))*weight(i) : i=1,...,ntau) .
!
!    The arrays TAU and WEIGHT are given in common block
!    DATA, as is the array GTAU containing the sequence
!    g(tau(i)), i=1,..., NTAU.
!
!    The relevant function values of the B-splines  b(i), i=1,...,n, are
!    supplied by the subprogram BSPLVB.
!
!    The coefficient matrix C, with
!           c(i,j)  := (b(i), b(j)), i,j=1,...,n,
!    of the normal equations is symmetric and (2*k-1)-banded, therefore
!    can be specified by giving its K bands at or below the diagonal.
!    For i=1,...,n,  we store
!    (b(i),b(j)) = c(i,j)  in  q(i-j+1,j), j=i,...,min(i+k-1,n)
!    and the right side
!    (b(i), g )  in  bcoef(i).
!
!    Since B-spline values are most efficiently generated by finding sim-
!    ultaneously the value of every nonzero B-spline at one point,
!    the entries of C (i.e., of  Q ), are generated by computing, for
!    each ll, all the terms involving  tau(ll)  simultaneously and adding
!    them to all relevant entries.
!
!  Reference:
!
!    Carl DeBoor,
!    A Practical Guide to Splines,
!    Springer Verlag.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) T(N+K), the knot sequence.
!
!    Input, integer N, the dimension of the space of splines of order K
!    with knots t.
!
!    Input, integer K, the order.
!
!  work arrays
!
!    q....a work array of size (at least) k*n. its first  k  rows are used
!       for the  k  lower diagonals of the gramian matrix  c.
!
!    diag.....a work array of length  n  used in bchfac .
!
!  input  via  c o m m o n  /data/
!
!  ntau.....number of data points
!  (tau(i),gtau(i)), i=1,...,ntau     are the  ntau  data points to be
!        fitted .
!  weight(i), i=1,...,ntau    are the corresponding weights .
!
!  output
!  bcoef(1), ..., bcoef(n)  the b-spline coefficients of the l2-appr.
!
  implicit none

  integer k
  integer n
  integer, parameter :: ntmax = 200

  real ( kind = 8 ) bcoef(n)
  real ( kind = 8 ) biatx(k)
  real ( kind = 8 ) diag(n)
  real ( kind = 8 ) dw
  real ( kind = 8 ) gtau
  integer i
  integer j
  integer jj
  integer left
  integer leftmk
  integer ll
  integer mm
  integer ntau
  real ( kind = 8 ) q(k,n)
  real ( kind = 8 ) t(n+k)
  real ( kind = 8 ) tau
  real ( kind = 8 ) totalw
  real ( kind = 8 ) weight

  COMMON /DATA/ tau(ntmax),gtau(ntmax),weight(ntmax),totalw,ntau

  bcoef(1:n) = 0.0D+00
  q(1:k,1:n) = 0.0D+00

  left = k
  leftmk = 0

  do ll = 1, ntau
!
!  Locate LEFT such that tau(ll) in (t(left),t(left+1)).
!
    do

      if ( left == n ) then
        exit
      end if

      if ( tau(ll) < t(left+1) ) then
        exit
      end if

      left = left + 1
      leftmk = leftmk + 1

    end do

    call bsplvb ( t, k, 1, tau(ll), left, biatx )
!
!  biatx(mm) contains the value of b(left-k+mm) at tau(ll).
!  hence, with  dw := biatx(mm)*weight(ll), the number dw*gtau(ll)
!  is a summand in the inner product
!     (b(left-k+mm), g)  which goes into  bcoef(left-k+mm)
!  and the number biatx(jj)*dw is a summand in the inner product
!     (b(left-k+jj), b(left-k+mm)), into  q(jj-mm+1,left-k+mm)
!  since  (left-k+jj)-(left-k+mm)+1 = jj - mm + 1 .
!
    do mm = 1, k

      dw = biatx(mm)*weight(ll)
      j = leftmk+mm
      bcoef(j) = dw*gtau(ll)+bcoef(j)
      i = 1

      do jj = mm, k
        q(i,j) = biatx(jj)*dw+q(i,j)
        i = i+1
      end do

    end do

  end do
!
!  Construct the Cholesky factorization for C in  q , then
!  use it to solve the normal equations
!    c*x = bcoef
!  for X, and store X in BCOEF.
!
  call bchfac ( q, k, n, diag )

  call bchslv ( q, k, n, bcoef )

  return
end