#include "lapacknames.inc"
SUBROUTINE SGETRF( M, N, A, LDA, IPIV, INFO )
  USE LA_CONSTANTS
  USE LA_AUXILIARY, ONLY: ILAENV, XERBLA, LA_LASWP
  USE LA_BLAS3, ONLY: LA_GEMM, LA_TRSM
  USE LA_XGETF2

!
!  -- LAPACK routine (version 3.0) --
!     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
!     Courant Institute, Argonne National Lab, and Rice University
!     March 31, 1993
!     04-09-02:  LAPACK 3E version (eca)
!
!     .. Scalar Arguments ..
      INTEGER            INFO, LDA, M, N
!     ..
!     .. Array Arguments ..
      INTEGER            IPIV( * )
      REAL(WP)           A( LDA, * )
!     ..
!
!  Purpose
!  =======
!
!  SGETRF computes an LU factorization of a general M-by-N matrix A
!  using partial pivoting with row interchanges.
!
!  The factorization has the form
!     A = P * L * U
!  where P is a permutation matrix, L is lower triangular with unit
!  diagonal elements (lower trapezoidal if m > n), and U is upper
!  triangular (upper trapezoidal if m < n).
!
!  This is the right-looking Level 3 BLAS version of the algorithm.
!
!  Arguments
!  =========
!
!  M       (input) INTEGER
!          The number of rows of the matrix A.  M >= 0.
!
!  N       (input) INTEGER
!          The number of columns of the matrix A.  N >= 0.
!
!  A       (input/output) REAL array, dimension (LDA,N)
!          On entry, the M-by-N matrix to be factored.
!          On exit, the factors L and U from the factorization
!          A = P*L*U; the unit diagonal elements of L are not stored.
!
!  LDA     (input) INTEGER
!          The leading dimension of the array A.  LDA >= max(1,M).
!
!  IPIV    (output) INTEGER array, dimension (min(M,N))
!          The pivot indices; for 1 <= i <= min(M,N), row i of the
!          matrix was interchanged with row IPIV(i).
!
!  INFO    (output) INTEGER
!          = 0:  successful exit
!          < 0:  if INFO = -i, the i-th argument had an illegal value
!          > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
!                has been completed, but the factor U is exactly
!                singular, and division by zero will occur if it is used
!                to solve a system of equations.
!
!  =====================================================================
!
!     .. Local Scalars ..
      INTEGER            I, IINFO, J, JB, NB
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC          MAX, MIN
!     ..
!     .. Executable Statements ..
!
!     Test the input parameters.
      INFO = 0
      IF( M < 0 ) THEN
         INFO = -1
      ELSE IF( N < 0 ) THEN
         INFO = -2
      ELSE IF( LDA < MAX(1,M) ) THEN
         INFO = -4
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( SPREFIX // 'GETRF', -INFO )
         RETURN
      END IF
      
      !     Determine the block size for this environment.
      NB = ILAENV( 1, SPREFIX // 'GETRF', ' ', M, N, -1, -1 )
      IF( NB <= 1 .OR. NB >= MIN(M,N) ) THEN
         !        Use unblocked code.
         CALL LA_GETF2( M, N, A, LDA, IPIV, INFO )
      ELSE
         !        Use blocked code.
         DO J = 1, MIN(M,N), NB
            JB = MIN(MIN(M,N)-J+1, NB)
            !           Factor diagonal and subdiagonal blocks and test for exact
            !           singularity.  
            CALL LA_GETF2( M-J+1, JB, A(J,J), LDA, IPIV(J), IINFO )
            !           Adjust INFO and the pivot indices.
            IF( INFO == 0 .AND. IINFO > 0 ) &
                 INFO = IINFO + J - 1
            DO I = J, MIN(M,J+JB-1)
               IPIV(I) = J - 1 + IPIV(I)
            END DO
            !           Apply interchanges to columns 1:J-1.
            CALL LA_LASWP( J-1, A(1,1), LDA, J, J+JB-1, IPIV(1), 1 )
            !
            IF( J+JB <= N ) THEN
               !        Apply interchanges to columns J+JB:N.
               CALL LA_LASWP( N-J-JB+1, A(1,J+JB), LDA, J, J+JB-1, &
                    IPIV(1), 1 )
               !        Compute block row of U.
               CALL LA_TRSM( 'Left', 'Lower', 'No transpose', 'Unit', &
                    JB, N-J-JB+1, ONE, A(J,J), LDA, &
                    A(J,J+JB), LDA )
               IF( J+JB <= M ) THEN
                  !     Update trailing submatrix.
                  CALL LA_GEMM( 'No transpose', 'No transpose', &
                       M-J-JB+1, N-J-JB+1, JB, -ONE, &
                       A(J+JB,J), LDA, A(J,J+JB), LDA, &
                       ONE, A(J+JB,J+JB), LDA )
               END IF
            END IF
         END DO
         
      END IF
      RETURN
      
      
!
!     End of SGETRF
!
      END