pure subroutine dgbfa(abd, lda, n, ml, mu, ipvt, info)
!! DGBFA factors a real band matrix by elimination.
!
! Discussion:
!
! DGBFA is usually called by DGBCO, but it can be called
! directly with a saving in time if RCOND is not needed.
!
! 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/output, real ( kind = 8 ) ABD(LDA,N). On input, the matrix in band
! storage. The columns of the matrix are stored in the columns of ABD
! and the diagonals of the matrix are stored in rows ML+1 through
! 2*ML+MU+1 of ABD. On output, an upper triangular matrix in band storage
! and the multipliers which were used to obtain it. The factorization
! can be written A = L*U where L is a product of permutation and unit lower
! triangular matrices and U is upper triangular.
!
! Input, integer ( kind = 4 ) LDA, the leading dimension of the array ABD.
! 2*ML + MU + 1 <= LDA is required.
!
! 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.
!
! Output, integer ( kind = 4 ) IPVT(N), the pivot indices.
!
! Output, integer ( kind = 4 ) INFO, error flag.
! 0, normal value.
! K, if U(K,K) == 0.0. This is not an error condition for this
! subroutine, but it does indicate that DGBSL will divide by zero if
! called. Use RCOND in DGBCO for a reliable indication of singularity.
use blas_interfaces, only: idamax, daxpy, dscal
real(rk), intent(inout) :: abd(lda, n)
integer, intent(in) :: lda
integer, intent(in) :: n
integer, intent(in) :: ml
integer, intent(in) :: mu
integer, intent(out) :: ipvt(n)
integer, intent(out) :: info
integer :: i
integer :: i0
integer :: j
integer :: j0
integer :: j1
integer :: ju
integer :: jz
integer :: k
integer :: l
integer :: lm
integer :: m
integer :: mm
real(rk) :: t
m = ml + mu + 1
info = 0
!
! Zero initial fill-in columns.
!
j0 = mu + 2
j1 = min(n, m) - 1
do jz = j0, j1
i0 = m + 1 - jz
do i = i0, ml
abd(i, jz) = zero
end do
end do
jz = j1
ju = 0
!
! Gaussian elimination with partial pivoting.
!
do k = 1, n - 1
!
! Zero out the next fill-in column.
!
jz = jz + 1
if (jz <= n) then
abd(1:ml, jz) = zero
end if
!
! Find L = pivot index.
!
lm = min(ml, n - k)
l = idamax(lm + 1, abd(m, k), 1) + m - 1
ipvt(k) = l + k - m
!
! Zero pivot implies this column already triangularized.
!
if (abd(l, k) == zero) then
info = k
!
! Interchange if necessary.
!
else
if (l /= m) then
t = abd(l, k)
abd(l, k) = abd(m, k)
abd(m, k) = t
end if
!
! Compute multipliers.
!
t = -one/abd(m, k)
call dscal(lm, t, abd(m + 1, k), 1)
!
! Row elimination with column indexing.
!
ju = min(max(ju, mu + ipvt(k)), n)
mm = m
do j = k + 1, ju
l = l - 1
mm = mm - 1
t = abd(l, j)
if (l /= mm) then
abd(l, j) = abd(mm, j)
abd(mm, j) = t
end if
call daxpy(lm, t, abd(m + 1, k), 1, abd(mm + 1, j), 1)
end do
end if
end do
ipvt(n) = n
if (abd(m, n) == zero) then
info = n
end if
end subroutine dgbfa