sgefa Subroutine

  subroutine sgefa(a,lda,n,ipvt,info)

    !! Routine to factor a real matrix by Gaussian elimination
    !! 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, matrix to be factored.
    !! On exit, 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.
    !! lda : input integer : leading dimension of the array A
    !! n : input integer : order of the matrix A
    !! ipvt(n) : output integer array : pivot indices
    !! info : output integer : info flag
    !! = 0  normal completion
    !! = k  if  u(k,k) == 0.0
    !! This routine factors a real matrix by Gaussian elimination.
    !!     !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    use global_variables, only: verbose
    implicit none

    !  Arguments

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

    !  Local variables

    integer :: j,k,kp1,l,nm1
    real(dp) :: t

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

    if (verbose == 1) then
       do j = 1,n
          if (all(a(j,:) == 0.0D0)) then
             write(*,*) 'Line ',j, &
                  ' in matrix a in subroutine sgefa is all zero'
          end if
       end do
    end if

    info = 0
    nm1 = n - 1

    if (nm1 >= 1) then

       do 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) /= 0.0D0) then

             !  Interchange if necessary

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

             !  Compute multipliers

             t = -1.0D0/a(k,k)
             call sscal(n-k,t,a(k+1:n,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
                end if
                call saxpy(n-k,t,a(k+1:n,k),1,a(k+1:n,j),1)
             end do

          else
             info = k
             if (verbose == 1) then
                write(*,*) 'a(l,k) = 0.0D0 in subroutine sgefa'
                write(*,*) 'info=k=',info
             end if
          end if
       end do

    end if

    ipvt(n) = n
    if (a(n,n) == 0.0D0) then
       info = n
       if (verbose == 1) then
          write(*,*) 'Error: a(n,n) == 0.0D0 in subroutine sgefa'
       end if
    end if

  end subroutine sgefa