SUBROUTINE sgefa (a, lda, n, ipvt, info)  
!     sgefa factors a real matrix by gaussian elimination.
!     sgefa is usually called by sgeco, but it can be called
!     directly with a saving in time if  rcond  is not needed.
!     (time for sgeco) = (1 + 9/n)*(time for sgefa) .
!
!     on entry
!        a       real(lda, n)
!                the matrix to be factored.
!        lda     integer
!                the leading dimension of the array  a .
!        n       integer
!                the order of the matrix  a .
!     on return
!        a       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.
!        ipvt    integer(n)
!                an integer vector of pivot indices.
!        info    integer
!                = 0  normal value.
!                = k  if  u(k,k) .eq. 0.0 .  this is not an error
!                     condition for this subroutine, but it does
!                     indicate that sgesl or sgedi will divide by zero
!                     if called.  use  rcond  in sgeco for a reliable
!                     indication of singularity.
!     linpack. this version dated 08/14/78 .
!     cleve moler, university of new mexico, argonne national lab.
!     subroutines and functions
!     blas saxpy,sscal,isamax
implicit none

REAL   , INTENT(IN OUT):: a(lda,1)
INTEGER, INTENT(IN)    :: n, lda
INTEGER, INTENT(OUT)   :: info, ipvt(1)
!     internal variables
real :: t
integer :: isamax, j, k, kp1, l, nm1
!    gaussian elimination with partial pivoting
info = 0
nm1 = n - 1
if (nm1.lt.1) goto 70
do 60 k = 1, nm1
   kp1 = k + 1
!        find l = pivot index
   l = isamax (n - k + 1, a (k, k), 1) + k - 1
   ipvt (k) = l
!        zero pivot implies this column already triangularized
   if (a (l, k) .eq.0.0e0) goto 40
!           interchange if necessary
   if (l.eq.k) goto 10
   t = a (l, k)
   a (l, k) = a (k, k)
   a (k, k) = t
   10    continue
!           compute multipliers
   t = - 1.0e0 / a (k, k)
   call sscal (n - k, t, a (k + 1, k), 1)
!           row elimination with column indexing
   do 30 j = kp1, n
      t = a (l, j)
      if (l.eq.k) goto 20
      a (l, j) = a (k, j)
      a (k, j) = t
   20       continue
      call saxpy (n - k, t, a (k + 1, k), 1, a (k + 1, j), 1)
   30    end do
   goto 50
   40    continue
   info = k
   50    continue
   60 end do
   70 continue
ipvt (n) = n
if (a (n, n) .eq.0.0e0) info = n

return
END SUBROUTINE sgefa
