banfac.f

Go to the documentation of this file.

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

c*********************************************************************72
c
cc BANFAC factors a banded matrix without pivoting.
c
c  from  * a practical guide to splines *  by c. de boor    
c  returns in  w  the lu-factorization (without pivoting) of the banded
c  matrix  a  of order  nrow  with  (nbandl + 1 + nbandu) bands or diag-
c  onals in the work array  w .
c
c******  i n p u t  ******
c  w.....work array of size  (nroww,nrow)  containing the interesting
c        part of a banded matrix  a , with the diagonals or bands of  a
c        stored in the rows of  w , while columns of  a  correspond to
c        columns of  w . this is the storage mode used in  linpack  and
c        results in efficient innermost loops.
c           explicitly,  a  has  nbandl  bands below the diagonal
c                            +     1     (main) diagonal
c                            +   nbandu  bands above the diagonal
c        and thus, with    middle = nbandu + 1,
c          a(i+j,j)  is in  w(i+middle,j)  for i=-nbandu,...,nbandl
c                                              j=1,...,nrow .
c        for example, the interesting entries of a (1,2)-banded matrix
c        of order  9  would appear in the first  1+1+2 = 4  rows of  w
c        as follows.
c
c
c
c
c
c        all other entries of  w  not identified in this way with an en-
c        try of  a  are never referenced .
c  nroww.....row dimension of the work array  w .
c        must be  .ge.  nbandl + 1 + nbandu  .
c  nbandl.....number of bands of  a  below the main diagonal
c  nbandu.....number of bands of  a  above the main diagonal .
c
c******  o u t p u t  ******
c  iflag.....integer indicating success( = 1) or failure ( = 2) .
c     if  iflag = 1, then
c  w.....contains the lu-factorization of  a  into a unit lower triangu-
c        lar matrix  l  and an upper triangular matrix  u (both banded)
c        and stored in customary fashion over the corresponding entries
c        of  a . this makes it possible to solve any particular linear
c        system  a*x = b  for  x  by a
c              call banslv ( w, nroww, nrow, nbandl, nbandu, b )
c        with the solution x  contained in  b  on return .
c     if  iflag = 2, then
c        one of  nrow-1, nbandl,nbandu failed to be nonnegative, or else
c        one of the potential pivots was found to be zero indicating
c        that  a  does not have an lu-factorization. this implies that
c        a  is singular in case it is totally positive .
c
c******  m e t h o d  ******
c     gauss elimination  w i t h o u t  pivoting is used. the routine is
c  intended for use with matrices  a  which do not require row inter-
c  changes during factorization, especially for the  t o t a l l y
c  p o s i t i v e  matrices which occur in spline calculations.
c     the routine should not be used for an arbitrary banded matrix.
c
      implicit none

      integer iflag,nbandl,nbandu,nrow,nroww,   i,ipk,j,jmax,k,kmax
     &                                        ,middle,midmk,nrowm1
      double precision w(nroww,nrow),   factor,pivot
c
      iflag = 1
      middle = nbandu + 1
c                         w(middle,.) contains the main diagonal of  a .
      nrowm1 = nrow - 1
      if (nrowm1)                       999,900,1
    1 if (nbandl .gt. 0)                go to 10
c                a is upper triangular. check that diagonal is nonzero .
      do 5 i=1,nrowm1
         if (w(middle,i) .eq. 0.0d+00)       go to 999
    5    continue
                                        go to 900
   10 if (nbandu .gt. 0)                go to 20
c              a is lower triangular. check that diagonal is nonzero and
c                 divide each column by its diagonal .
      do 15 i=1,nrowm1
         pivot = w(middle,i)
         if(pivot .eq. 0.0d+00)              go to 999
         jmax = min0(nbandl, nrow - i)
         do 15 j=1,jmax
   15       w(middle+j,i) = w(middle+j,i)/pivot
                                        return
c
c        a  is not just a triangular matrix. construct lu factorization
   20 do 50 i=1,nrowm1
c                                  w(middle,i)  is pivot for i-th step .
         pivot = w(middle,i)
         if (pivot .eq. 0.0d+00)             go to 999
c                 jmax  is the number of (nonzero) entries in column  i
c                     below the diagonal .
         jmax = min0(nbandl,nrow - i)
c              divide each entry in column  i  below diagonal by pivot .
         do 32 j=1,jmax
   32       w(middle+j,i) = w(middle+j,i)/pivot
c                 kmax  is the number of (nonzero) entries in row  i  to
c                     the right of the diagonal .
         kmax = min0(nbandu,nrow - i)
c                  subtract  a(i,i+k)*(i-th column) from (i+k)-th column
c                  (below row  i ) .
         do 40 k=1,kmax
            ipk = i + k
            midmk = middle - k
            factor = w(midmk,ipk)
            do 40 j=1,jmax
   40          w(midmk+j,ipk) = w(midmk+j,ipk) - w(middle+j,i)*factor
   50    continue
c                                       check the last diagonal entry .
  900 if (w(middle,nrow) .ne. 0.0d+00)       return
  999 iflag = 2
                                        return
      end