Homec4science
Diffusion SPClibs (master) Browse

spclibs/pppack/bchslv.f90
93a1e9712d04master

Last Change
Raw File

bchslv.f90
View Options

!>
!> @file bchslv.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 bchslv ( w, nbands, nrow, b )
!*************************************************************************
!
!! BCHSLV solves a banded symmetric positive definite system.
!
! Discussion:
!
! The system is of the form:
!
! C * X = B
!
! and the Cholesky factorization of C has been constructed
! by BCHFAC.
!
! With the factorization
!
! C = L * D * L'
!
! available, where L is unit lower triangular and D is diagonal,
! the triangular system
!
! L * Y = B
!
! is solved for Y (forward substitution), Y is stored in B, the
! vector D**(-1)*Y is computed and stored in B, then the
! triangular system L'*X = D**(-1)*Y is solved for X
! (backsubstitution).
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! Input, real ( kind = 8 ) W(NBANDS,NROW), the Cholesky factorization for C,
! as computed by BCHFAC.
!
! Input, integer NBANDS, the bandwidth of C.
!
! Input, integer NROW, the order of the matrix C.
!
! Input/output, real ( kind = 8 ) B(NROW).
! On input, the right hand side.
! On output, the solution.
!
implicit none
integer nbands
integer nrow
real ( kind = 8 ) b(nrow)
integer j
integer n
real ( kind = 8 ) w(nbands,nrow)
if ( nrow <= 1 ) then
b(1) = b(1) * w(1,1)
return
end if
!
! Forward substitution.
! Solve L*Y=B.
!
do n = 1, nrow
do j = 1, min(nbands-1,nrow-n)
b(j+n) = b(j+n) - w(j+1,n) * b(n)
end do
end do
!
! Backsubstitution.
! Solve L'*X=D**(-1)*Y.
!
do n = nrow, 1, -1
b(n) = b(n)*w(1,n)
do j = 1, min(nbands-1,nrow-n)
b(n) = b(n) - w(j+1,n) * b(j+n)
end do
end do
return
end
c4science · Help