SUBROUTINE sgedi (a, lda, n, ipvt, det, work, job)  
!     sgedi computes the determinant and inverse of a matrix
!     using the factors computed by sgeco or sgefa.
!
!     on entry
!        a       real(lda, n)
!                the output from sgeco or sgefa.
!        lda     integer
!                the leading dimension of the array  a .
!        n       integer
!                the order of the matrix  a .
!        ipvt    integer(n)
!                the pivot vector from sgeco or sgefa.
!        work    real(n)
!                work vector.  contents destroyed.
!        job     integer
!                = 11   both determinant and inverse.
!                = 01   inverse only.
!                = 10   determinant only.
!     on return
!        a       inverse of original matrix if requested.
!                otherwise unchanged.
!        det     real(2)
!                determinant of original matrix if requested.
!                otherwise not referenced.
!                determinant = det(1) * 10.0**det(2)
!                with  1.0 .le. abs(det(1)) .lt. 10.0
!                or  det(1) .eq. 0.0 .
!
!     error condition
!        a division by zero will occur if the input factor contains
!        a zero on the diagonal and the inverse is requested.
!        it will not occur if the subroutines are called correctly
!        and if sgeco has set rcond .gt. 0.0 or sgefa has set
!        info .eq. 0 .
!     linpack. this version dated 08/14/78 .
!     cleve moler, university of new mexico, argonne national lab.
!
!     subroutines and functions
!     blas saxpy,sscal,sswap
!     fortran abs,mod
implicit none

INTEGER, INTENT(IN)    :: n, job, lda, ipvt(1)
REAL   , INTENT(IN OUT):: a(lda,1)
REAL   , INTENT(OUT)   :: det(2), work(1)
!     internal variables
real    :: t, ten
integer :: i, j, k, kb, kp1, l, nm1
!     compute determinant
if (job / 10.eq.0) goto 70
det (1) = 1.0e0
det (2) = 0.0e0
ten = 10.0e0
do 50 i = 1, n
   if (ipvt (i) .ne.i) det (1) = - det (1)
   det (1) = a (i, i) * det (1)
!        ...exit
   if (det (1) .eq.0.0e0) goto 60
   10    if (abs (det (1) ) .ge.1.0e0) goto 20
   det (1) = ten * det (1)
   det (2) = det (2) - 1.0e0
   goto 10
   20    continue
   30    if (abs (det (1) ) .lt.ten) goto 40
   det (1) = det (1) / ten
   det (2) = det (2) + 1.0e0
   goto 30
   40    continue
   50 end do
   60 continue
   70 continue
!     compute inverse(u)
if (mod (job, 10) .eq.0) goto 150
do 100 k = 1, n
   a (k, k) = 1.0e0 / a (k, k)
   t = - a (k, k)
   call sscal (k - 1, t, a (1, k), 1)
   kp1 = k + 1
   if (n.lt.kp1) goto 90
   do 80 j = kp1, n
      t = a (k, j)
      a (k, j) = 0.0e0
      call saxpy (k, t, a (1, k), 1, a (1, j), 1)
   80    end do
   90    continue
  100 end do
!        form inverse(u)*inverse(l)
nm1 = n - 1
if (nm1.lt.1) goto 140
do 130 kb = 1, nm1
   k = n - kb
   kp1 = k + 1
   do 110 i = kp1, n
      work (i) = a (i, k)
      a (i, k) = 0.0e0
  110    end do
   do 120 j = kp1, n
      t = work (j)
      call saxpy (n, t, a (1, j), 1, a (1, k), 1)
  120    end do
   l = ipvt (k)
   if (l.ne.k) call sswap (n, a (1, k), 1, a (1, l), 1)
  130 end do
  140 continue
  150 continue

return
END SUBROUTINE sgedi
