sgedi Subroutine

  subroutine sgedi(a,lda,n,ipvt,det,job)

    !! Routine to compute the determinant and inverse of a matrix
    !! author: Cleve Moler, University of New Mexico, Argonne National Lab.
    !! author: P J Knight, CCFE, Culham Science Centre
    !! a(lda,n) : input/output real array :
    !! On entry, output from <A HREF="sgefa.html">sgefa</A>.
    !! On exit, the inverse if requested, otherwise unchanged
    !! lda      : input integer : leading dimension of the array A
    !! n        : input integer : order of the matrix A
    !! ipvt(n)  : input integer array : pivot vector from sgefa
    !! det(2)   : output real array : 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 .
    !! job : input integer : switch for required outputs
    !! = 11   both determinant and inverse.
    !! = 01   inverse only.
    !! = 10   determinant only.
    !! This routine computes the determinant and inverse of a matrix
    !! using the factors computed by (SGECO or) <A HREF="sgefa.html">SGEFA</A>.
    !! <P>A division by zero will occur if the input factor contains
    !! a zero on the diagonal and the inverse is requested.
    !!     !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    implicit none

    !  Arguments

    integer, intent(in) :: lda,n,job
    integer, dimension(n), intent(in) :: ipvt
    real(dp), dimension(lda,n), intent(inout) :: a

    !  Local variables

    integer :: i,j,k,kk,kb,kp1,l,nm1
    real(dp), parameter :: ten = 10.0D0
    real(dp) :: t
    real(dp), dimension(2) :: det
    real(dp), dimension(n) :: work

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    if ((job/10) /= 0) then  !  Compute determinant

       det(1) = 1.0D0
       det(2) = 0.0D0

       do i = 1, n
          if (ipvt(i) /=  i) det(1) = -det(1)
          det(1) = a(i,i)*det(1)

          if (det(1) == 0.0D0) exit

          do
             if (abs(det(1)) >= 1.0D0) exit
             det(1) = ten*det(1)
             det(2) = det(2) - 1.0D0
          end do

          do
             if (abs(det(1)) < ten) exit
             det(1) = det(1)/ten
             det(2) = det(2) + 1.0D0
          end do
       end do

    end if

    !  Compute inverse(u)

    if (mod(job,10) /= 0) then

       do k = 1, n
          a(k,k) = 1.0D0/a(k,k)
          t = -a(k,k)

          call sscal(k-1,t,a(1:n,k),1)
          kp1 = k + 1
          if (n >= kp1) then
             do j = kp1, n
                t = a(k,j)
                a(k,j) = 0.0D0
                kk = k
                call saxpy(kk,t,a(1:n,k),1,a(1:n,j),1)
             end do
          end if
       end do

       !  Form inverse(u)*inverse(l)

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

             do i = kp1, n
                work(i) = a(i,k)
                a(i,k) = 0.0D0
             end do

             do j = kp1, n
                t = work(j)
                call saxpy(n,t,a(1:n,j),1,a(1:n,k),1)
             end do

             l = ipvt(k)

             if (l /= k) call sswap(n,a(1:n,k),1,a(1:n,l),1)
          end do

       end if
    end if

  end subroutine sgedi