pure subroutine dgbsl(abd, lda, n, ml, mu, ipvt, b, job)
!! [[dgbsl]] solves a real banded system factored by [[dgbco]] or [[dgbfa]].
!
! Discussion:
!
! DGBSL can solve either A * X = B or A' * X = B.
!
! A division by zero will occur if the input factor contains a
! zero on the diagonal. Technically this indicates singularity
! but it is often caused by improper arguments or improper
! setting of LDA. It will not occur if the subroutines are
! called correctly and if DGBCO has set 0.0 < RCOND
! or DGBFA has set INFO == 0.
!
! To compute inverse(A) * C where C is a matrix with P columns:
!
! call dgbco ( abd, lda, n, ml, mu, ipvt, rcond, z )
!
! if ( rcond is too small ) then
! exit
! end if
!
! do j = 1, p
! call dgbsl ( abd, lda, n, ml, mu, ipvt, c(1,j), 0 )
! end do
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 17 May 2005
!
! Author:
!
! FORTRAN90 version by John Burkardt.
!
! Reference:
!
! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
! LINPACK User's Guide,
! SIAM, 1979,
! ISBN13: 978-0-898711-72-1,
! LC: QA214.L56.
!
! Parameters:
!
! Input, real ( kind = 8 ) ABD(LDA,N), the output from DGBCO or DGBFA.
!
! Input, integer ( kind = 4 ) LDA, the leading dimension of the array ABD.
!
! Input, integer ( kind = 4 ) N, the order of the matrix.
!
! Input, integer ( kind = 4 ) ML, MU, the number of diagonals below and above
! the main diagonal. 0 <= ML < N, 0 <= MU < N.
!
! Input, integer ( kind = 4 ) IPVT(N), the pivot vector from DGBCO or DGBFA.
!
! Input/output, real ( kind = 8 ) B(N). On input, the right hand side.
! On output, the solution.
!
! Input, integer ( kind = 4 ) JOB, job choice.
! 0, solve A*X=B.
! nonzero, solve A'*X=B.
use blas_interfaces, only: daxpy, ddot
real(rk), intent(in) :: abd(lda, n)
integer, intent(in) :: lda
integer, intent(in) :: n
integer, intent(in) :: ml
integer, intent(in) :: mu
integer, intent(in) :: ipvt(n)
real(rk), intent(inout) :: b(n)
integer, intent(in) :: job
integer :: k
integer :: l
integer :: la
integer :: lb
integer :: lm
integer :: m
real(rk) :: t
m = mu + ml + 1
!
! JOB = 0, Solve A * x = b.
!
! First solve L * y = b.
!
if (job == 0) then
if (0 < ml) then
do k = 1, n - 1
lm = min(ml, n - k)
l = ipvt(k)
t = b(l)
if (l /= k) then
b(l) = b(k)
b(k) = t
end if
call daxpy(lm, t, abd(m + 1, k), 1, b(k + 1), 1)
end do
end if
!
! Now solve U * x = y.
!
do k = n, 1, -1
b(k) = b(k)/abd(m, k)
lm = min(k, m) - 1
la = m - lm
lb = k - lm
t = -b(k)
call daxpy(lm, t, abd(la, k), 1, b(lb), 1)
end do
!
! JOB nonzero, solve A' * x = b.
!
! First solve U' * y = b.
!
else
do k = 1, n
lm = min(k, m) - 1
la = m - lm
lb = k - lm
t = ddot(lm, abd(la, k), 1, b(lb), 1)
b(k) = (b(k) - t)/abd(m, k)
end do
!
! Now solve L' * x = y.
!
if (0 < ml) then
do k = n - 1, 1, -1
lm = min(ml, n - k)
b(k) = b(k) + ddot(lm, abd(m + 1, k), 1, b(k + 1), 1)
l = ipvt(k)
if (l /= k) then
t = b(l)
b(l) = b(k)
b(k) = t
end if
end do
end if
end if
end subroutine dgbsl