dgbfa Subroutine

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

factor a band matrix using gaussian elimination.

dgbfa is usually called by dgbco, but it can be called directly with a saving in time if rcond is not needed.

Band storage

if a is a band matrix, the following program segment will set up the input.

  ml = (band width below the diagonal)
  mu = (band width above the diagonal)
  m = ml + mu + 1
  do j = 1, n
     i1 = max(1, j-mu)
     i2 = min(n, j+ml)
     do i = i1, i2
        k = i - j + m
        abd(k,j) = a(i,j)
     end do
  end do

this uses rows ml+1 through 2*ml+mu+1 of abd . in addition, the first ml rows in abd are used for elements generated during the triangularization. the total number of rows needed in abd is 2*ml+mu+1. the ml+mu by ml+mu upper left triangle and the ml by ml lower right triangle are not referenced.

Author

  • moler, c. b., (u. of new mexico)

Revision history

  • 780814 date written
  • 890531 changed all specific intrinsics to generic. (wrb)
  • 890831 modified array declarations. (wrb)
  • 890831 revision date from version 3.2
  • 891214 prologue converted to version 4.0 format. (bab)
  • 900326 removed duplicate information from description section. (wrb)
  • 920501 reformatted the references section. (wrb)

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(inout) :: abd(lda,*)
  • input: contains 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 . see the comments below for details.
  • 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.
integer, intent(in) :: lda

the leading dimension of the array abd . lda must be >= 2*ml + mu + 1 .
integer, intent(in) :: n

the order of the original matrix.
integer, intent(in) :: ml

number of diagonals below the main diagonal. 0 <= ml < n .
integer, intent(in) :: mu

number of diagonals above the main diagonal. 0 <= mu < n . more efficient if ml <= mu .
integer, intent(out) :: ipvt(*)

an integer vector of pivot indices.
integer, intent(out) :: info
  • 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.

Calls

proc~~dgbfa~~CallsGraph proc~dgbfa dvode_linpack_module::dgbfa proc~daxpy dvode_blas_module::daxpy proc~dgbfa->proc~daxpy proc~dscal dvode_blas_module::dscal proc~dgbfa->proc~dscal proc~idamax dvode_blas_module::idamax proc~dgbfa->proc~idamax

Called by

proc~~dgbfa~~CalledByGraph proc~dgbfa dvode_linpack_module::dgbfa proc~dvjac dvode_module::dvode_t%dvjac proc~dvjac->proc~dgbfa proc~dvnlsd dvode_module::dvode_t%dvnlsd proc~dvnlsd->proc~dvjac proc~dvstep dvode_module::dvode_t%dvstep proc~dvstep->proc~dvnlsd proc~dvode dvode_module::dvode_t%dvode proc~dvode->proc~dvstep

Source Code

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

       integer,intent(in) :: lda !! the leading dimension of the array  `abd` .
                                 !! `lda` must be `>= 2*ml + mu + 1` .
       integer,intent(in) :: n !! the order of the original matrix.
       integer,intent(in) :: ml !! number of diagonals below the main diagonal.
                                !! `0 <= ml < n` .
       integer,intent(in) :: mu !! number of diagonals above the main diagonal.
                                !! `0 <= mu < n` .
                                !! more efficient if  `ml <= mu` .
       integer,intent(out) :: ipvt(*) !! an integer vector of pivot indices.
       integer,intent(out) :: info !! * 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.
       real(wp),intent(inout) :: abd(lda,*) !! * **input:** contains 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` .
                                            !!   see the comments below for details.
                                            !! * **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.


       real(wp) :: t
       integer :: i , i0 , j , ju , jz , j0 , j1 , k , kp1 , l , &
                  lm , m , mm , nm1

#ifdef HAS_BLAS
       ! user is linking against an external BLAS library
       integer,external :: idamax 
#endif

       m = ml + mu + 1
       info = 0

      ! zero initial fill-in columns

       j0 = mu + 2
       j1 = min(n,m) - 1
       if ( j1>=j0 ) then
          do jz = j0 , j1
             i0 = m + 1 - jz
             do i = i0 , ml
                abd(i,jz) = 0.0_wp
             enddo
          enddo
       endif
       jz = j1
       ju = 0

      ! gaussian elimination with partial pivoting

       nm1 = n - 1
       if ( nm1>=1 ) then
          do k = 1 , nm1
             kp1 = k + 1

            ! zero next fill-in column

             jz = jz + 1
             if ( jz<=n ) then
                if ( ml>=1 ) then
                   do i = 1 , ml
                      abd(i,jz) = 0.0_wp
                   enddo
                endif
             endif

            ! 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)==0.0_wp ) then
                info = k
             else

               ! interchange if necessary

                if ( l/=m ) then
                   t = abd(l,k)
                   abd(l,k) = abd(m,k)
                   abd(m,k) = t
                endif

               ! compute multipliers

                t = -1.0_wp/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
                if ( ju>=kp1 ) then
                   do j = kp1 , 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
                      endif
                      call daxpy(lm,t,abd(m+1,k),1,abd(mm+1,j),1)
                   enddo
                endif
             endif
          enddo
       endif
       ipvt(n) = n
       if ( abd(m,n)==0.0_wp ) info = n

   end subroutine dgbfa