dbnfac Subroutine

private pure 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(kind=ip), intent(in) :: nroww

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

matrix order
integer(kind=ip), intent(in) :: nbandl

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

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

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

Called by

proc~~dbnfac~~CalledByGraph proc~dbnfac bspline_sub_module::dbnfac proc~dbint4 bspline_sub_module::dbint4 proc~dbint4->proc~dbnfac proc~dbintk bspline_sub_module::dbintk proc~dbintk->proc~dbnfac proc~db1ink_alt bspline_sub_module::db1ink_alt proc~db1ink_alt->proc~dbint4 proc~db1ink_alt_2 bspline_sub_module::db1ink_alt_2 proc~db1ink_alt_2->proc~dbint4 proc~dbtpcf bspline_sub_module::dbtpcf proc~dbtpcf->proc~dbintk interface~db1ink bspline_sub_module::db1ink interface~db1ink->proc~db1ink_alt interface~db1ink->proc~db1ink_alt_2 proc~db1ink_default bspline_sub_module::db1ink_default interface~db1ink->proc~db1ink_default proc~db1ink_default->proc~dbtpcf proc~db2ink bspline_sub_module::db2ink proc~db2ink->proc~dbtpcf proc~db3ink bspline_sub_module::db3ink proc~db3ink->proc~dbtpcf proc~db4ink bspline_sub_module::db4ink proc~db4ink->proc~dbtpcf proc~db5ink bspline_sub_module::db5ink proc~db5ink->proc~dbtpcf proc~db6ink bspline_sub_module::db6ink proc~db6ink->proc~dbtpcf proc~initialize_1d_auto_knots bspline_oo_module::bspline_1d%initialize_1d_auto_knots proc~initialize_1d_auto_knots->interface~db1ink proc~initialize_1d_specify_knots bspline_oo_module::bspline_1d%initialize_1d_specify_knots proc~initialize_1d_specify_knots->interface~db1ink proc~initialize_2d_auto_knots bspline_oo_module::bspline_2d%initialize_2d_auto_knots proc~initialize_2d_auto_knots->proc~db2ink proc~initialize_2d_specify_knots bspline_oo_module::bspline_2d%initialize_2d_specify_knots proc~initialize_2d_specify_knots->proc~db2ink proc~initialize_3d_auto_knots bspline_oo_module::bspline_3d%initialize_3d_auto_knots proc~initialize_3d_auto_knots->proc~db3ink proc~initialize_3d_specify_knots bspline_oo_module::bspline_3d%initialize_3d_specify_knots proc~initialize_3d_specify_knots->proc~db3ink proc~initialize_4d_auto_knots bspline_oo_module::bspline_4d%initialize_4d_auto_knots proc~initialize_4d_auto_knots->proc~db4ink proc~initialize_4d_specify_knots bspline_oo_module::bspline_4d%initialize_4d_specify_knots proc~initialize_4d_specify_knots->proc~db4ink proc~initialize_5d_auto_knots bspline_oo_module::bspline_5d%initialize_5d_auto_knots proc~initialize_5d_auto_knots->proc~db5ink proc~initialize_5d_specify_knots bspline_oo_module::bspline_5d%initialize_5d_specify_knots proc~initialize_5d_specify_knots->proc~db5ink proc~initialize_6d_auto_knots bspline_oo_module::bspline_6d%initialize_6d_auto_knots proc~initialize_6d_auto_knots->proc~db6ink proc~initialize_6d_specify_knots bspline_oo_module::bspline_6d%initialize_6d_specify_knots proc~initialize_6d_specify_knots->proc~db6ink proc~bspline_1d_constructor_auto_knots bspline_oo_module::bspline_1d_constructor_auto_knots proc~bspline_1d_constructor_auto_knots->proc~initialize_1d_auto_knots proc~bspline_1d_constructor_specify_knots bspline_oo_module::bspline_1d_constructor_specify_knots proc~bspline_1d_constructor_specify_knots->proc~initialize_1d_specify_knots proc~bspline_2d_constructor_auto_knots bspline_oo_module::bspline_2d_constructor_auto_knots proc~bspline_2d_constructor_auto_knots->proc~initialize_2d_auto_knots proc~bspline_2d_constructor_specify_knots bspline_oo_module::bspline_2d_constructor_specify_knots proc~bspline_2d_constructor_specify_knots->proc~initialize_2d_specify_knots proc~bspline_3d_constructor_auto_knots bspline_oo_module::bspline_3d_constructor_auto_knots proc~bspline_3d_constructor_auto_knots->proc~initialize_3d_auto_knots proc~bspline_3d_constructor_specify_knots bspline_oo_module::bspline_3d_constructor_specify_knots proc~bspline_3d_constructor_specify_knots->proc~initialize_3d_specify_knots proc~bspline_4d_constructor_auto_knots bspline_oo_module::bspline_4d_constructor_auto_knots proc~bspline_4d_constructor_auto_knots->proc~initialize_4d_auto_knots proc~bspline_4d_constructor_specify_knots bspline_oo_module::bspline_4d_constructor_specify_knots proc~bspline_4d_constructor_specify_knots->proc~initialize_4d_specify_knots proc~bspline_5d_constructor_auto_knots bspline_oo_module::bspline_5d_constructor_auto_knots proc~bspline_5d_constructor_auto_knots->proc~initialize_5d_auto_knots proc~bspline_5d_constructor_specify_knots bspline_oo_module::bspline_5d_constructor_specify_knots proc~bspline_5d_constructor_specify_knots->proc~initialize_5d_specify_knots proc~bspline_6d_constructor_auto_knots bspline_oo_module::bspline_6d_constructor_auto_knots proc~bspline_6d_constructor_auto_knots->proc~initialize_6d_auto_knots proc~bspline_6d_constructor_specify_knots bspline_oo_module::bspline_6d_constructor_specify_knots proc~bspline_6d_constructor_specify_knots->proc~initialize_6d_specify_knots interface~bspline_1d bspline_oo_module::bspline_1d interface~bspline_1d->proc~bspline_1d_constructor_auto_knots interface~bspline_1d->proc~bspline_1d_constructor_specify_knots interface~bspline_2d bspline_oo_module::bspline_2d interface~bspline_2d->proc~bspline_2d_constructor_auto_knots interface~bspline_2d->proc~bspline_2d_constructor_specify_knots interface~bspline_3d bspline_oo_module::bspline_3d interface~bspline_3d->proc~bspline_3d_constructor_auto_knots interface~bspline_3d->proc~bspline_3d_constructor_specify_knots interface~bspline_4d bspline_oo_module::bspline_4d interface~bspline_4d->proc~bspline_4d_constructor_auto_knots interface~bspline_4d->proc~bspline_4d_constructor_specify_knots interface~bspline_5d bspline_oo_module::bspline_5d interface~bspline_5d->proc~bspline_5d_constructor_auto_knots interface~bspline_5d->proc~bspline_5d_constructor_specify_knots interface~bspline_6d bspline_oo_module::bspline_6d interface~bspline_6d->proc~bspline_6d_constructor_auto_knots interface~bspline_6d->proc~bspline_6d_constructor_specify_knots

Source Code

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

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

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

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

    if (nrowm1 < 0_ip) then
        iflag = 2_ip
        return
    else if (nrowm1 == 0_ip) then
        if (w(middle,nrow)==0.0_wp) iflag = 2_ip
        return
    end if

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

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

    ! a is not just a triangular matrix. construct lu factorization
    do i=1_ip,nrowm1
        ! w(middle,i)  is pivot for i-th step .
        pivot = w(middle,i)
        if (pivot==0.0_wp) then
            iflag = 2_ip
            return
        end if
        ! 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_ip,jmax
            w(middle+j,i) = w(middle+j,i)/pivot
        end do
        ! 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_ip,kmax
            ipk = i + k
            midmk = middle - k
            factor = w(midmk,ipk)
            do j=1_ip,jmax
                w(midmk+j,ipk) = w(midmk+j,ipk) - w(middle+j,i)*factor
            end do
        end do
    end do

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

    end subroutine dbnfac