      SUBROUTINE SLUBL( M, N, A, LDA, IPIV, INFO )
*
*  -- LAPACK routine --
*     Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab,
*     Courant Institute, NAG Ltd., and Rice University
*     August 16, 1989
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, M, N
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * )
      REAL               A( LDA, * )
*     ..
*
*  Purpose
*  =======
*
*  SLUBL 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 left-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 x n matrix to be factored.
*          On exit, the factors L and 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.  Row i of the matrix was interchanged with
*          row IPIV(i).
*
*  INFO    (output) INTEGER
*          = 0: successful exit
*          < 0: if INFO = -k, the k-th argument had an illegal value
*          > 0: if INFO = k, U(k,k) 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 or to compute the inverse
*               of A.
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ONE
      PARAMETER          ( ONE = 1.0E+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, J, JB, NB
*     ..
*     .. External Subroutines ..
      EXTERNAL           ENVIR, SGEMM, SLASWP, SLUL, STRSM, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX, MIN
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      INFO = 0
      IF( M.LT.0 ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -4
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SLUBL', -INFO )
         RETURN
      END IF
*
*     Quick return if possible.
*
      IF( M.EQ.0 .OR. N.EQ.0 )
     $   RETURN
*
*     Determine the block size for this environment.
*
      CALL ENVIR( 'Block', NB )
      IF( NB.EQ.1 ) THEN
         CALL SLUL( M, N, A, LDA, IPIV, INFO )
      ELSE
*
         DO 20 J = 1, MIN( M, N ), NB
            JB = MIN( MIN( M, N )-J+1, NB )
*
*           Apply previous interchanges to current block.
*
            CALL SLASWP( JB, A( 1, J ), LDA, 1, J-1, IPIV, 1 )
*
*           Compute superdiagonal block of U.
*
            CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit', J-1,
     $                  JB, ONE, A, LDA, A( 1, J ), LDA )
*
*           Update diagonal and subdiagonal blocks.
*
            CALL SGEMM( 'No transpose', 'No transpose', M-J+1, JB, J-1,
     $                  -ONE, A( J, 1 ), LDA, A( 1, J ), LDA, ONE,
     $                  A( J, J ), LDA )
*
*           Factorize diagonal and subdiagonal blocks and test for exact
*           singularity.
*
            CALL SLUL( M-J+1, JB, A( J, J ), LDA, IPIV( J ), INFO )
            DO 10 I = J, J + JB - 1
               IPIV( I ) = J - 1 + IPIV( I )
   10       CONTINUE
            IF( INFO.EQ.0 ) THEN
*
*              Apply interchanges in current block to previous blocks.
*
               CALL SLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 )
            ELSE
*
*              If INFO is not zero, a zero pivot was found in SLUL.
*              Correct the index returned from SLUL and go on.
*
               INFO = INFO + J - 1
            END IF
   20    CONTINUE
*
         IF( N.GT.M ) THEN
*
*           Apply interchanges from columns 1 to m to columns m+1 to n.
*
            CALL SLASWP( N-M, A( 1, M+1 ), LDA, 1, M, IPIV, 1 )
*
*           Complete columns m+1 to n of U.
*
            CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit', M, N-M,
     $                  ONE, A, LDA, A( 1, M+1 ), LDA )
         END IF
      END IF
      RETURN
*
*     End of SLUBL
*
      END
      SUBROUTINE SLUL( M, N, A, LDA, IPIV, INFO )
*
*  -- LAPACK routine --
*     Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab,
*     Courant Institute, NAG Ltd., and Rice University
*     August 16, 1989
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, M, N
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * )
      REAL               A( LDA, * )
*     ..
*
*  Purpose
*  =======
*
*  SLUL 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 left-looking Level 2 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 x n matrix to be factored.
*          On exit, the factors L and 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.  Row i of the matrix was interchanged with
*          row IPIV(i).
*
*  INFO    (output) INTEGER
*          = 0: successful exit
*          < 0: if INFO = -k, the k-th argument had an illegal value
*          > 0: if INFO = k, U(k,k) 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 or to compute the inverse
*               of A.
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ONE, ZERO
      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, IP, J, JP, LEN
      REAL               T
*     ..
*     .. External Functions ..
      INTEGER            ISAMAX
      EXTERNAL           ISAMAX
*     ..
*     .. External Subroutines ..
      EXTERNAL           SGEMV, SSCAL, SSWAP, STRSV, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX, MIN
*     ..
*     .. Executable Statements ..
*
*
*     Test the input parameters.
*
      INFO = 0
      IF( M.LT.0 ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -4
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SLUL', -INFO )
         RETURN
      END IF
*
*     Quick return if possible.
*
      IF( M.EQ.0 .OR. N.EQ.0 )
     $   RETURN
*
      DO 20 J = 1, N
*
*        Apply previous interchanges to J-th column.
*
         LEN = MIN( J-1, M )
         DO 10 I = 1, LEN
            IP = IPIV( I )
            IF( IP.NE.I ) THEN
               T = A( I, J )
               A( I, J ) = A( IP, J )
               A( IP, J ) = T
            END IF
   10    CONTINUE
*
*        Compute elements 1:LEN of J-th column.
*
         CALL STRSV( 'Lower', 'No transpose', 'Unit', LEN, A, LDA,
     $               A( 1, J ), 1 )
*
*        Update elements J:M of J-th column.
*
         IF( J.LE.M ) THEN
            CALL SGEMV( 'No transpose', M-J+1, J-1, -ONE, A( J, 1 ),
     $                  LDA, A( 1, J ), 1, ONE, A( J, J ), 1 )
*
*           Find pivot and test for singularity.
*
            JP = J - 1 + ISAMAX( M-J+1, A( J, J ), 1 )
            IPIV( J ) = JP
            IF( A( JP, J ).NE.ZERO ) THEN
*
*              Apply interchange to columns 1:J.
*
               IF( JP.NE.J )
     $            CALL SSWAP( J, A( J, 1 ), LDA, A( JP, 1 ), LDA )
*
*              Compute elements J+1:M of J-th column.
*
               IF( J.LT.M )
     $            CALL SSCAL( M-J, ONE/A( J, J ), A( J+1, J ), 1 )
            ELSE
*
*              If A( JP, J ) is zero, set INFO to indicate that a zero
*              pivot has been found.  INFO returns the index of the
*              last zero pivot to the calling routine if there is more
*              than one.
*
               INFO = J
            END IF
         END IF
   20 CONTINUE
      RETURN
*
*     End of SLUL
*
      END
      SUBROUTINE SLUBC( M, N, A, LDA, IPIV, INFO )
*
*  -- LAPACK routine --
*     Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab,
*     Courant Institute, NAG Ltd., and Rice University
*     August 16, 1989
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, M, N
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * )
      REAL               A( LDA, * )
*     ..
*
*  Purpose
*  =======
*
*  SLUBC 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 Level 3 BLAS version of the Crout 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 x n matrix to be factored.
*          On exit, the factors L and 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.  Row i of the matrix was interchanged with
*          row IPIV(i).
*
*  INFO    (output) INTEGER
*          = 0: successful exit
*          < 0: if INFO = -k, the k-th argument had an illegal value
*          > 0: if INFO = k, U(k,k) 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 or to compute the inverse
*               of A.
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ONE
      PARAMETER          ( ONE = 1.0E+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, J, JB, NB
*     ..
*     .. External Subroutines ..
      EXTERNAL           ENVIR, SGEMM, SLASWP, SLUC, STRSM
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MIN
*     ..
*     .. Executable Statements ..
*
      INFO = 0
*
*     Determine the block size for this environment.
*
      CALL ENVIR( 'Block', NB )
      IF( NB.EQ.1 ) THEN
         CALL SLUC( M, N, A, LDA, IPIV, INFO )
      ELSE
*
         DO 20 J = 1, MIN( M, N ), NB
            JB = MIN( MIN( M, N )-J+1, NB )
*
*           Update diagonal and subdiagonal blocks.
*
            CALL SGEMM( 'No transpose', 'No transpose', M-J+1, JB, J-1,
     $                  -ONE, A( J, 1 ), LDA, A( 1, J ), LDA, ONE,
     $                  A( J, J ), LDA )
*
*           Factorize diagonal and subdiagonal blocks and test for
*           singularity.
*
            CALL SLUC( M-J+1, JB, A( J, J ), LDA, IPIV( J ), INFO )
*
*           Update pivot indices and apply the interchanges to the
*           columns on either side of the current block.
*
            DO 10 I = J, MIN( M, J+JB-1 )
               IPIV( I ) = J - 1 + IPIV( I )
   10       CONTINUE
            CALL SLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 )
            CALL SLASWP( N-J-JB+1, A( 1, J+JB ), LDA, J, J+JB-1, IPIV,
     $                   1 )
*
            IF( INFO.NE.0 )
     $         GO TO 30
*
            IF( J+JB.LE.N ) THEN
*
*              Compute block row of U.
*
               CALL SGEMM( 'No transpose', 'No transpose', JB, N-J-JB+1,
     $                     J-1, -ONE, A( J, 1 ), LDA, A( 1, J+JB ), LDA,
     $                     ONE, A( J, J+JB ), LDA )
               CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB,
     $                     N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ),
     $                     LDA )
            END IF
   20    CONTINUE
      END IF
      RETURN
*
   30 CONTINUE
      INFO = INFO + J - 1
      RETURN
*
*     End of SLUBC
*
      END
      SUBROUTINE SLUC( M, N, A, LDA, IPIV, INFO )
*
*  -- LAPACK routine --
*     Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab,
*     Courant Institute, NAG Ltd., and Rice University
*     August 16, 1989
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, M, N
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * )
      REAL               A( LDA, * )
*     ..
*
*  Purpose
*  =======
*
*  SLUC 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 Level 2 BLAS version of the Crout 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 x n matrix to be factored.
*          On exit, the factors L and 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.  Row i of the matrix was interchanged with
*          row IPIV(i).
*
*  INFO    (output) INTEGER
*          = 0: successful exit
*          < 0: if INFO = -k, the k-th argument had an illegal value
*          > 0: if INFO = k, U(k,k) 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 or to compute the inverse
*               of A.
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ONE, ZERO
      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            J, JP
*     ..
*     .. External Functions ..
      INTEGER            ISAMAX
      EXTERNAL           ISAMAX
*     ..
*     .. External Subroutines ..
      EXTERNAL           SGEMV, SSCAL, SSWAP
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MIN
*     ..
*     .. Executable Statements ..
*
      INFO = 0
*
      DO 10 J = 1, MIN( M, N )
*
*        Update diagonal and subdiagonal elements in column J.
*
         CALL SGEMV( 'No transpose', M-J+1, J-1, -ONE, A( J, 1 ), LDA,
     $               A( 1, J ), 1, ONE, A( J, J ), 1 )
*
*        Find pivot and test for singularity.
*
         JP = J - 1 + ISAMAX( M-J+1, A( J, J ), 1 )
         IPIV( J ) = JP
         IF( A( JP, J ).NE.ZERO ) THEN
*
*           Apply interchange to columns 1:J.
*
            IF( JP.NE.J )
     $         CALL SSWAP( N, A( J, 1 ), LDA, A( JP, 1 ), LDA )
*
*           Compute elements J+1:M of J-th column.
*
            IF( J.LT.M )
     $         CALL SSCAL( M-J, ONE/A( J, J ), A( J+1, J ), 1 )
         ELSE
*
*           If A( JP, J ) is zero, set INFO to indicate that a zero
*           pivot has been found.  INFO returns the index of the
*           last zero pivot to the calling routine if there is more
*           than one.
*
            INFO = J
         END IF
*
         IF( J+1.LE.N ) THEN
*
*           Compute block row of U.
*
            CALL SGEMV( 'Transpose', J-1, N-J, -ONE, A( 1, J+1 ), LDA,
     $                  A( J, 1 ), LDA, ONE, A( J, J+1 ), LDA )
         END IF
   10 CONTINUE
      RETURN
*
*     End of SLUC
*
      END
      SUBROUTINE SLUBR( M, N, A, LDA, IPIV, INFO )
*
*  -- LAPACK routine --
*     Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab,
*     Courant Institute, NAG Ltd., and Rice University
*     August 16, 1989
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, M, N
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * )
      REAL               A( LDA, * )
*     ..
*
*  Purpose
*  =======
*
*  SLUBR 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 x n matrix to be factored.
*          On exit, the factors L and 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.  Row i of the matrix was interchanged with
*          row IPIV(i).
*
*  INFO    (output) INTEGER
*          = 0: successful exit
*          < 0: if INFO = -k, the k-th argument had an illegal value
*          > 0: if INFO = k, U(k,k) 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 or to compute the inverse
*               of A.
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ONE
      PARAMETER          ( ONE = 1.0E+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, J, JB, NB
*     ..
*     .. External Subroutines ..
      EXTERNAL           ENVIR, SGEMM, SLASWP, SLUR, STRSM
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MIN
*     ..
*     .. Executable Statements ..
*
      INFO = 0
*
*     Determine the block size for this environment.
*
      CALL ENVIR( 'Block', NB )
      IF( NB.EQ.1 ) THEN
         CALL SLUR( M, N, A, LDA, IPIV, INFO )
      ELSE
*
         DO 20 J = 1, MIN( M, N ), NB
            JB = MIN( MIN( M, N )-J+1, NB )
*
*           Factorize diagonal and subdiagonal blocks and test for
*           singularity.
*
            CALL SLUR( M-J+1, JB, A( J, J ), LDA, IPIV( J ), INFO )
            DO 10 I = J, J + JB - 1
               IPIV( I ) = J - 1 + IPIV( I )
   10       CONTINUE
            IF( INFO.NE.0 )
     $         GO TO 30
*
*           Apply interchanges to previous blocks.
*
            CALL SLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 )
*
*           Apply interchanges to following blocks.
*
            CALL SLASWP( N-J-JB+1, A( 1, J+JB ), LDA, J, J+JB-1, IPIV,
     $                   1 )
            IF( J+JB.LE.N ) THEN
*
*              Compute block row of U.
*
               CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB,
     $                     N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ),
     $                     LDA )
               IF( J+JB.LE.M ) THEN
*
*                 Update trailing submatrix.
*
                  CALL SGEMM( '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
   20    CONTINUE
      END IF
      RETURN
*
   30 CONTINUE
      INFO = INFO + J - 1
      RETURN
*
*     End of SLUBR
*
      END
      SUBROUTINE SLUR( M, N, A, LDA, IPIV, INFO )
*
*  -- LAPACK routine --
*     Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab,
*     Courant Institute, NAG Ltd., and Rice University
*     August 16, 1989
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, M, N
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * )
      REAL               A( LDA, * )
*     ..
*
*  Purpose
*  =======
*
*  SLUR 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 2 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 x n matrix to be factored.
*          On exit, the factors L and 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.  Row i of the matrix was interchanged with
*          row IPIV(i).
*
*  INFO    (output) INTEGER
*          = 0: successful exit
*          < 0: if INFO = -k, the k-th argument had an illegal value
*          > 0: if INFO = k, U(k,k) 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 or to compute the inverse
*               of A.
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ONE, ZERO
      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            J, JP
*     ..
*     .. External Functions ..
      INTEGER            ISAMAX
      EXTERNAL           ISAMAX
*     ..
*     .. External Subroutines ..
      EXTERNAL           SGER, SSCAL, SSWAP
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MIN
*     ..
*     .. Executable Statements ..
*
      INFO = 0
*
      DO 10 J = 1, MIN( M, N )
*
*        Find pivot and test for singularity.
*
         JP = J - 1 + ISAMAX( M-J+1, A( J, J ), 1 )
         IPIV( J ) = JP
         IF( A( JP, J ).NE.ZERO ) THEN
*
*           Interchange rows J and JP.
*
            IF( JP.NE.J )
     $         CALL SSWAP( N, A( J, 1 ), LDA, A( JP, 1 ), LDA )
*
*           Compute elements J+1:M of J-th column.
*
            IF( J.LT.M )
     $         CALL SSCAL( M-J, ONE/A( J, J ), A( J+1, J ), 1 )
         ELSE
*
*           If A( JP, J ) is zero, set INFO to indicate that a zero
*           pivot has been found.  INFO returns the index of the last
*           zero pivot to the calling routine if there is more than one.
*
            INFO = J
         END IF
*
         IF( J+1.LE.N ) THEN
*
*           Update trailing submatrix.
*
            CALL SGER( M-J, N-J, -ONE, A( J+1, J ), 1, A( J, J+1 ), LDA,
     $                 A( J+1, J+1 ), LDA )
         END IF
   10 CONTINUE
      RETURN
*
*     End of SLUR
*
      END
      SUBROUTINE ENVIR( WHAT, NB )
*
*  -- LAPACK auxiliary routine --
*     Argonne National Lab, Courant Institute, and N.A.G. Ltd.
*     April 1, 1989
*
*     .. Scalar Arguments ..
      CHARACTER         WHAT
      INTEGER           NB
*     ..
*
*  Purpose
*  =======
*
*     ENVIR returns the block size to be used by the block algorithms.
*
*  Arguments
*  =========
*
*  WHAT   - CHARACTER
*           WHAT is a character code for the value to be returned.
*
*  NB     - INTEGER
*           The parameter NB returns an integer value from common
*           storage.
*
*     .. External Functions ..
      LOGICAL           LSAME
      EXTERNAL          LSAME
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC         MAX
*     ..
*     .. Scalars in Common ..
      INTEGER           NBLOCK, NPROC, NSHIFT
*     ..
*     .. Common blocks ..
      COMMON            / CENVIR / NBLOCK, NPROC, NSHIFT
*     ..
*     .. Executable Statements ..
      IF( LSAME( WHAT, 'B' ) ) THEN
         NB = MAX( 1, NBLOCK )
      ELSE IF( LSAME( WHAT, 'P' ) ) THEN
         NB = NPROC
      ELSE IF( LSAME( WHAT, 'S' ) ) THEN
         NB = NSHIFT
      END IF
      RETURN
*
*     End of ENVIR
*
      END