dgefa Subroutine

public subroutine dgefa(a, lda, n, ipvt, info)

factor a matrix using gaussian elimination.

dgefa is usually called by dgeco, but it can be called directly with a saving in time if rcond is not needed. (time for dgeco) = (1 + 9/n)*(time for dgefa) .

Author

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

Revision history

  • 780814 date written
  • 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) :: a(lda,*)
  • input: the matrix to be factored.
  • output: an upper triangular matrix 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 a.
integer, intent(in) :: n

the order of the matrix a.
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 dgesl or dgedi will divide by zero if called. use rcond in dgeco for a reliable indication of singularity.

Calls

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

Called by

proc~~dgefa~~CalledByGraph proc~dgefa dvode_linpack_module::dgefa proc~dvjac dvode_module::dvode_t%dvjac proc~dvjac->proc~dgefa 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 dgefa(a,lda,n,ipvt,info)

       integer,intent(in) :: lda !! the leading dimension of the array `a`.
       integer,intent(in) :: n !! the order of the matrix `a`.
       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 dgesl or dgedi will divide by zero
                                   !!      if called.  use  rcond  in dgeco for a reliable
                                   !!      indication of singularity.
       real(wp),intent(inout) :: a(lda,*) !! * **input:** the matrix to be factored.
                                          !! * **output:** an upper triangular matrix 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 :: j , k , kp1 , l , nm1

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

      ! gaussian elimination with partial pivoting

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

            ! find l = pivot index

             l = idamax(n-k+1,a(k,k),1) + k - 1
             ipvt(k) = l

            ! zero pivot implies this column already triangularized

             if ( a(l,k)==0.0_wp ) then
                info = k
             else

               ! interchange if necessary

                if ( l/=k ) then
                   t = a(l,k)
                   a(l,k) = a(k,k)
                   a(k,k) = t
                endif

               ! compute multipliers

                t = -1.0_wp/a(k,k)
                call dscal(n-k,t,a(k+1,k),1)

               ! row elimination with column indexing

                do j = kp1 , n
                   t = a(l,j)
                   if ( l/=k ) then
                      a(l,j) = a(k,j)
                      a(k,j) = t
                   endif
                   call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
                enddo
             endif
          enddo
       endif
       ipvt(n) = n
       if ( a(n,n)==0.0_wp ) info = n

   end subroutine dgefa