!>
!> @file bchfac.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 bchfac ( w, nbands, nrow, diag )

!*************************************************************************
!
!! BCHFAC constructs a Cholesky factorization of a matrix.
!
!  Discussion:
!
!    The factorization has the form
!
!      C = L * D * L'
!
!    with L unit lower triangular and D diagonal, for a given matrix C of
!    order NROW, where C is symmetric positive semidefinite and banded,
!    having NBANDS diagonals at and below the main diagonal.
!
!    Gauss elimination is used, adapted to the symmetry and bandedness of C.
!
!    Near-zero pivots are handled in a special way.  The diagonal
!    element C(N,N)=W(1,N) is saved initially in DIAG(N), all N.
!
!    At the N-th elimination step, the current pivot element, W(1,N),
!    is compared with its original value, DIAG(N).  If, as the result
!    of prior elimination steps, this element has been reduced by about
!    a word length, (i.e., if W(1,N)+DIAG(N) <= DIAG(N)), then the pivot
!    is declared to be zero, and the entire N-th row is declared to
!    be linearly dependent on the preceding rows.  This has the effect
!    of producing X(N) = 0 when solving C*X = B for X, regardless of B.
!
!    Justification for this is as follows.  In contemplated applications
!    of this program, the given equations are the normal equations for
!    some least-squares approximation problem, DIAG(N) = C(N,N) gives
!    the norm-square of the N-th basis function, and, at this point,
!    W(1,N) contains the norm-square of the error in the least-squares
!    approximation to the N-th basis function by linear combinations
!    of the first N-1.
!
!    Having W(1,N)+DIAG(N) <= DIAG(N) signifies that the N-th function
!    is linearly dependent to machine accuracy on the first N-1
!    functions, therefore can safely be left out from the basis of
!    approximating functions.
!
!    The solution of a linear system C*X=B is effected by the
!    succession of the following two calls:
!
!      CALL BCHFAC(W,NBANDS,NROW,DIAG)
!
!      CALL BCHSLV(W,NBANDS,NROW,B,X)
!
!  Reference:
!
!    Carl DeBoor,
!    A Practical Guide to Splines,
!    Springer Verlag.
!
!  Parameters:
!
!    Input/output, real ( kind = 8 ) W(NBANDS,NROW).
!
!    On input, W contains the NBANDS diagonals in its rows,
!    with the main diagonal in row 1.  Precisely, W(I,J)
!    contains C(I+J-1,J), I=1,...,NBANDS, J=1,...,NROW.
!
!    For example, the interesting entries of a seven diagonal
!    symmetric matrix C of order 9 would be stored in W as
!
!      11 22 33 44 55 66 77 88 99
!      21 32 43 54 65 76 87 98  *
!      31 42 53 64 75 86 97  *  *
!      41 52 63 74 85 96  *  *  *
!
!    Entries of the array not associated with an
!    entry of C are never referenced.
!
!    On output, W contains the Cholesky factorization
!    C = L*D*L-transp, with W(1,I) containing 1/D(I,I) and W(I,J)
!    containing L(I-1+J,J), I=2,...,NBANDS.
!
!    Input, integer NBANDS, indicates the bandwidth of the
!    matrix C, i.e., C(I,J) = 0 for NBANDS < ABS(I-J).
!
!    Input, integer NROW, is the order of the matrix C.
!
!    Work array, real ( kind = 8 ) DIAG(NROW).
!
  implicit none

  integer nbands
  integer nrow

  real ( kind = 8 ) diag(nrow)
  integer i
  integer imax
  integer j
  integer jmax
  integer n
  real ( kind = 8 ) ratio
  real ( kind = 8 ) w(nbands,nrow)

  if ( nrow <= 1 ) then
    if ( 0.0D+00 < w(1,1) ) then
      w(1,1) = 1.0D+00 / w(1,1)
    end if
    return
  end if
!
!  Store the diagonal.
!
  diag(1:nrow) = w(1,1:nrow)
!
!  Factorization.
!
  do n = 1, nrow

    if ( w(1,n) + diag(n) <= diag(n) ) then
      w(1:nbands,n) = 0.0D+00
    else

      w(1,n) = 1.0D+00 / w(1,n)

      imax = min ( nbands-1, nrow-n )

      jmax = imax

      do i = 1, imax

        ratio = w(i+1,n) * w(1,n)

        do j = 1, jmax
          w(j,n+i) = w(j,n+i) - w(j+i,n) * ratio
        end do

        jmax = jmax-1
        w(i+1,n) = ratio

      end do

    end if

  end do

  return
end