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

!*************************************************************************
!
!! SUBBAK carries out backsubstitution for the current block.
!
!  Reference:
!
!    Carl DeBoor,
!    A Practical Guide to Splines,
!    Springer Verlag.
!
!  Parameters:
!
!    w, ipivot, nrow, ncol, last  are as on return from factrb.
!
!    x(1),...,x(ncol)  contains, on input, the right side for the
!            equations in this block after backsubstitution has been
!            carried up to but not including equation ipivot(last).
!            means that x(j) contains the right side of equation ipi-
!            vot(j) as modified during elimination, j=1,...,last, while
!            for j > last, x(j) is already a component of the solut-
!            ion vector.
!
!    x(1),...,x(ncol) contains, on output, the components of the solut-
!            ion corresponding to the present block.
!
  implicit none

  integer ncol
  integer nrow

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

  do k = last, 1, -1

    ip = ipivot(k)

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

    x(k) = ( x(k) - s ) / w(ip,k)

  end do

end