dbnfac Subroutine

private subroutine dbnfac(w, nroww, nrow, nbandl, nbandu, iflag)

Returns in w the LU-factorization (without pivoting) of the banded matrix a of order nrow with (nbandl + 1 + nbandu) bands or diagonals in the work array w .

gauss elimination without pivoting is used. the routine is intended for use with matrices a which do not require row inter- changes during factorization, especially for the totally positive matrices which occur in spline calculations. the routine should not be used for an arbitrary banded matrix.

Work array

Input

    w array of size (nroww,nrow) contains the interesting
    part of a banded matrix  a , with the diagonals or bands of  a
    stored in the rows of  w , while columns of  a  correspond to
    columns of  w . this is the storage mode used in  linpack  and
    results in efficient innermost loops.
       explicitly,  a  has  nbandl  bands below the diagonal
                        +     1     (main) diagonal
                        +   nbandu  bands above the diagonal
    and thus, with    middle = nbandu + 1,
      a(i+j,j)  is in  w(i+middle,j)  for i=-nbandu,...,nbandl
                                          j=1,...,nrow .
    for example, the interesting entries of a (1,2)-banded matrix
    of order  9  would appear in the first  1+1+2 = 4  rows of  w
    as follows.
                      13 24 35 46 57 68 79
                   12 23 34 45 56 67 78 89
                11 22 33 44 55 66 77 88 99
                21 32 43 54 65 76 87 98

    all other entries of  w  not identified in this way with an en-
    try of  a  are never referenced .

Output

  • if iflag = 1, then w contains the lu-factorization of a into a unit lower triangu- lar matrix l and an upper triangular matrix u (both banded) and stored in customary fashion over the corresponding entries of a . this makes it possible to solve any particular linear system a*x = b for x by a call dbnslv ( w, nroww, nrow, nbandl, nbandu, b ) with the solution x contained in b on return .
  • if iflag = 2, then one of nrow-1, nbandl,nbandu failed to be nonnegative, or else one of the potential pivots was found to be zero indicating that a does not have an lu-factorization. this implies that a is singular in case it is totally positive .

History

  • banfac written by carl de boor [5]
  • dbnfac from CMLIB [1]
  • Jacob Williams, 5/10/2015 : converted to free-form Fortran.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(inout), dimension(nroww,nrow) :: w

work array. See header for details.
integer, intent(in) :: nroww

row dimension of the work array w. must be >= nbandl + 1 + nbandu.
integer, intent(in) :: nrow

matrix order
integer, intent(in) :: nbandl

number of bands of a below the main diagonal
integer, intent(in) :: nbandu

number of bands of a above the main diagonal
integer, intent(out) :: iflag

indicating success(=1) or failure (=2)

Called by

proc~~dbnfac~~CalledByGraph proc~dbnfac dbnfac proc~dbintk dbintk proc~dbintk->proc~dbnfac proc~dbtpcf dbtpcf proc~dbtpcf->proc~dbintk proc~db1ink db1ink proc~db1ink->proc~dbtpcf proc~db2ink db2ink proc~db2ink->proc~dbtpcf proc~interp_surf interp_surf proc~interp_surf->proc~db2ink program~test_bspline test_bspline program~test_bspline->proc~interp_surf

Source Code

    subroutine dbnfac(w,nroww,nrow,nbandl,nbandu,iflag)

    integer,intent(in) :: nroww   !! row dimension of the work array w. must be >= nbandl + 1 + nbandu.
    integer,intent(in) :: nrow    !! matrix order
    integer,intent(in) :: nbandl  !! number of bands of a below the main diagonal
    integer,intent(in) :: nbandu  !! number of bands of a above the main diagonal
    integer,intent(out) :: iflag  !! indicating success(=1) or failure (=2)
    real(wp),dimension(nroww,nrow),intent(inout) :: w  !! work array. See header for details.

    integer :: i, ipk, j, jmax, k, kmax, middle, midmk, nrowm1
    real(wp) :: factor, pivot

    iflag = 1
    middle = nbandu + 1   ! w(middle,.) contains the main diagonal of a.
    nrowm1 = nrow - 1

    if (nrowm1 < 0) then
        iflag = 2
        return
    elseif (nrowm1 == 0) then
        if (w(middle,nrow)==0.0_wp) iflag = 2
        return
    endif

    if (nbandl<=0) then
        ! a is upper triangular. check that diagonal is nonzero .
        do i=1,nrowm1
            if (w(middle,i)==0.0_wp) then
                iflag = 2
                return
            endif
        enddo
        if (w(middle,nrow)==0.0_wp) iflag = 2
        return
    endif

    if (nbandu<=0) then
        ! a is lower triangular. check that diagonal is nonzero and
        ! divide each column by its diagonal.
        do i=1,nrowm1
            pivot = w(middle,i)
            if (pivot==0.0_wp) then
                iflag = 2
                return
            endif
            jmax = min(nbandl,nrow-i)
            do j=1,jmax
                w(middle+j,i) = w(middle+j,i)/pivot
            enddo
        enddo
        return
    endif

    ! a is not just a triangular matrix. construct lu factorization
    do i=1,nrowm1
        ! w(middle,i)  is pivot for i-th step .
        pivot = w(middle,i)
        if (pivot==0.0_wp) then
            iflag = 2
            return
        endif
        ! jmax is the number of (nonzero) entries in column i
        ! below the diagonal.
        jmax = min(nbandl,nrow-i)
        ! divide each entry in column i below diagonal by pivot.
        do j=1,jmax
            w(middle+j,i) = w(middle+j,i)/pivot
        enddo
        ! kmax is the number of (nonzero) entries in row i to
        ! the right of the diagonal.
        kmax = min(nbandu,nrow-i)
        ! subtract a(i,i+k)*(i-th column) from (i+k)-th column
        ! (below row i).
        do k=1,kmax
            ipk = i + k
            midmk = middle - k
            factor = w(midmk,ipk)
            do j=1,jmax
                w(midmk+j,ipk) = w(midmk+j,ipk) - w(middle+j,i)*factor
            enddo
        enddo
    enddo

    ! check the last diagonal entry.
    if (w(middle,nrow)==0.0_wp) iflag = 2

    endsubroutine dbnfac