dgbfa Subroutine

private pure subroutine dgbfa(abd, n, ml, mu, ipvt, tol, info)

Factors a banded matrix by Gaussian elimination with partial pivoting. The matrix A is stored in abd with A(i,j) in abd(i-j+ml+mu+1,j). Rows 1:ml of abd are used for the fill-in generated by pivoting.

History

  • Based on the LINPACK routine DGBFA (Cleve Moler, 1978). Added the tol argument for detecting singular matrices.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(inout), dimension(2_ip*ml+mu+1_ip,n) :: abd

IN: the matrix in band storage. OUT: its LU factorization.
integer(kind=ip), intent(in) :: n

order of the matrix
integer(kind=ip), intent(in) :: ml

number of diagonals below the main diagonal
integer(kind=ip), intent(in) :: mu

number of diagonals above the main diagonal
integer(kind=ip), intent(out), dimension(n) :: ipvt

pivot indices
real(kind=wp), intent(in) :: tol

a pivot with abs(u(k,k)) <= tol is treated as zero (tol=0 gives the original LINPACK behavior)
integer(kind=ip), intent(out) :: info

0 if successful, otherwise k if u(k,k) is zero (the matrix is singular)

Called by

proc~~dgbfa~~CalledByGraph proc~dgbfa dgbfa proc~dbint4 dbint4 proc~dbint4->proc~dgbfa proc~db1ink_alt db1ink_alt proc~db1ink_alt->proc~dbint4 proc~db1ink_alt_2 db1ink_alt_2 proc~db1ink_alt_2->proc~dbint4 interface~db1ink db1ink interface~db1ink->proc~db1ink_alt interface~db1ink->proc~db1ink_alt_2 proc~initialize_1d_auto_knots bspline_1d%initialize_1d_auto_knots proc~initialize_1d_auto_knots->interface~db1ink proc~initialize_1d_auto_knots_multi bspline_1d%initialize_1d_auto_knots_multi proc~initialize_1d_auto_knots_multi->interface~db1ink proc~initialize_1d_specify_knots bspline_1d%initialize_1d_specify_knots proc~initialize_1d_specify_knots->interface~db1ink proc~initialize_1d_specify_knots_multi bspline_1d%initialize_1d_specify_knots_multi proc~initialize_1d_specify_knots_multi->interface~db1ink proc~bspline_1d_constructor_auto_knots bspline_1d_constructor_auto_knots proc~bspline_1d_constructor_auto_knots->proc~initialize_1d_auto_knots proc~bspline_1d_constructor_specify_knots bspline_1d_constructor_specify_knots proc~bspline_1d_constructor_specify_knots->proc~initialize_1d_specify_knots interface~bspline_1d bspline_1d interface~bspline_1d->proc~bspline_1d_constructor_auto_knots interface~bspline_1d->proc~bspline_1d_constructor_specify_knots

Source Code

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

    implicit none

    integer(ip),intent(in)                :: n     !! order of the matrix
    integer(ip),intent(in)                :: ml    !! number of diagonals below the main diagonal
    integer(ip),intent(in)                :: mu    !! number of diagonals above the main diagonal
    real(wp),dimension(2_ip*ml+mu+1_ip,n),intent(inout) :: abd  !! **IN:** the matrix in band storage.
                                                                !! **OUT:** its LU factorization.
    integer(ip),dimension(n),intent(out)  :: ipvt  !! pivot indices
    real(wp),intent(in)                   :: tol   !! a pivot with `abs(u(k,k)) <= tol` is
                                                   !! treated as zero (`tol=0` gives the
                                                   !! original LINPACK behavior)
    integer(ip),intent(out)               :: info  !! 0 if successful, otherwise `k` if `u(k,k)`
                                                   !! is zero (the matrix is singular)

    integer(ip) :: i, j, jk, ju, jz, l, lm, m, mm
    real(wp) :: tmp

    m = ml + mu + 1_ip
    info = 0_ip

    ! zero the initial fill-in columns
    do jz = mu+2_ip, min(n,m)-1_ip
        abd(m+1_ip-jz:ml,jz) = 0.0_wp
    end do
    jz = min(n,m) - 1_ip
    ju = 0_ip

    do jk = 1_ip, n-1_ip

        ! zero the next fill-in column
        jz = jz + 1_ip
        if (jz<=n .and. ml>=1_ip) abd(1:ml,jz) = 0.0_wp

        ! find the pivot
        lm = min(ml,n-jk)
        l = maxloc(abs(abd(m:m+lm,jk)),dim=1,kind=ip) + m - 1_ip
        ipvt(jk) = l + jk - m

        if (abs(abd(l,jk))<=tol) then
            ! this column is already triangularized
            info = jk
            cycle
        end if

        ! interchange if necessary
        if (l/=m) then
            tmp = abd(l,jk)
            abd(l,jk) = abd(m,jk)
            abd(m,jk) = tmp
        end if

        ! compute the multipliers
        abd(m+1_ip:m+lm,jk) = abd(m+1_ip:m+lm,jk) * (-1.0_wp/abd(m,jk))

        ! row elimination with column indexing
        ju = min(max(ju,mu+ipvt(jk)),n)
        mm = m
        do j = jk+1_ip, ju
            l = l - 1_ip
            mm = mm - 1_ip
            tmp = abd(l,j)
            if (l/=mm) then
                abd(l,j) = abd(mm,j)
                abd(mm,j) = tmp
            end if
            do i = 1_ip, lm
                abd(mm+i,j) = abd(mm+i,j) + tmp*abd(m+i,jk)
            end do
        end do

    end do

    ipvt(n) = n
    if (abs(abd(m,n))<=tol) info = n

    end subroutine dgbfa