dgbfa Subroutine

public pure subroutine dgbfa(abd, lda, n, ml, mu, ipvt, info)

Uses

  • proc~~dgbfa~~UsesGraph proc~dgbfa dgbfa module~blas_interfaces blas_interfaces proc~dgbfa->module~blas_interfaces module~daskr_kinds daskr_kinds module~blas_interfaces->module~daskr_kinds iso_fortran_env iso_fortran_env module~daskr_kinds->iso_fortran_env

DGBFA factors a real band matrix by elimination.

Arguments

Type IntentOptional Attributes Name
real(kind=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

Calls

proc~~dgbfa~~CallsGraph proc~dgbfa dgbfa interface~daxpy daxpy proc~dgbfa->interface~daxpy interface~dscal dscal proc~dgbfa->interface~dscal interface~idamax idamax proc~dgbfa->interface~idamax

Called by

proc~~dgbfa~~CalledByGraph proc~dgbfa dgbfa proc~dmatd dmatd proc~dmatd->proc~dgbfa proc~jac_banpre jac_banpre proc~jac_banpre->proc~dgbfa

Variables

Type Visibility Attributes Name Initial
integer, public :: i
integer, public :: i0
integer, public :: j
integer, public :: j0
integer, public :: j1
integer, public :: ju
integer, public :: jz
integer, public :: k
integer, public :: l
integer, public :: lm
integer, public :: m
integer, public :: mm
real(kind=rk), public :: t

Source Code

   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