!>
!> @file subfor.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 subfor ( w, ipivot, nrow, last, b, x )

!*************************************************************************
!
!! SUBFOR carries out the forward pass of substitution for the current block.
!
!  Discussion:
!
!    The forward pass is the action on the right side corresponding to the
!    elimination carried out in FACTRB for this block.
!
!    At the end, x(j) contains the right side of the transformed
!    ipivot(j)-th equation in this block, j=1,...,nrow. then, since
!    for i=1,...,nrow-last, b(nrow+i) is going to be used as the right
!    side of equation I in the next block (shifted over there from
!    this block during factorization), it is set equal to x(last+i) here.
!
!  Reference:
!
!    Carl DeBoor,
!    A Practical Guide to Splines,
!    Springer Verlag.
!
!  Parameters:
!
!    w, ipivot, nrow, last  are as on return from factrb.
!
!    b(j)   is expected to contain, on input, the right side of j-th
!           equation for this block, j=1,...,nrow.
!    b(nrow+j)   contains, on output, the appropriately modified right
!           side for equation j in next block, j=1,...,nrow-last.
!
!    x(j)   contains, on output, the appropriately modified right
!           side of equation ipivot(j) in this block, j=1,...,last (and
!           even for j=last+1,...,nrow).
!
  implicit none

  integer last
  integer nrow

  real ( kind = 8 ) b(nrow+nrow-last)
  integer ip
  integer ipivot(nrow)
  integer j
  integer k
  real ( kind = 8 ) s
  real ( kind = 8 ) w(nrow,last)
  real ( kind = 8 ) x(nrow)

  ip = ipivot(1)
  x(1) = b(ip)

  do k = 2, nrow

    ip = ipivot(k)

    s = 0.0D+00
    do j = 1, min ( k-1, last )
      s = s + w(ip,j) * x(j)
    end do

    x(k) = b(ip) - s

  end do
!
!  Transfer modified right sides of equations ipivot(last+1),...,
!  ipivot(nrow) to next block.
!
  do k = last+1, nrow
    b(nrow-last+k) = x(k)
  end do

  return
end