001:       SUBROUTINE DGETF2( M, N, A, LDA, IPIV, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       INTEGER            INFO, LDA, M, N
009: *     ..
010: *     .. Array Arguments ..
011:       INTEGER            IPIV( * )
012:       DOUBLE PRECISION   A( LDA, * )
013: *     ..
014: *
015: *  Purpose
016: *  =======
017: *
018: *  DGETF2 computes an LU factorization of a general m-by-n matrix A
019: *  using partial pivoting with row interchanges.
020: *
021: *  The factorization has the form
022: *     A = P * L * U
023: *  where P is a permutation matrix, L is lower triangular with unit
024: *  diagonal elements (lower trapezoidal if m > n), and U is upper
025: *  triangular (upper trapezoidal if m < n).
026: *
027: *  This is the right-looking Level 2 BLAS version of the algorithm.
028: *
029: *  Arguments
030: *  =========
031: *
032: *  M       (input) INTEGER
033: *          The number of rows of the matrix A.  M >= 0.
034: *
035: *  N       (input) INTEGER
036: *          The number of columns of the matrix A.  N >= 0.
037: *
038: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
039: *          On entry, the m by n matrix to be factored.
040: *          On exit, the factors L and U from the factorization
041: *          A = P*L*U; the unit diagonal elements of L are not stored.
042: *
043: *  LDA     (input) INTEGER
044: *          The leading dimension of the array A.  LDA >= max(1,M).
045: *
046: *  IPIV    (output) INTEGER array, dimension (min(M,N))
047: *          The pivot indices; for 1 <= i <= min(M,N), row i of the
048: *          matrix was interchanged with row IPIV(i).
049: *
050: *  INFO    (output) INTEGER
051: *          = 0: successful exit
052: *          < 0: if INFO = -k, the k-th argument had an illegal value
053: *          > 0: if INFO = k, U(k,k) is exactly zero. The factorization
054: *               has been completed, but the factor U is exactly
055: *               singular, and division by zero will occur if it is used
056: *               to solve a system of equations.
057: *
058: *  =====================================================================
059: *
060: *     .. Parameters ..
061:       DOUBLE PRECISION   ONE, ZERO
062:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
063: *     ..
064: *     .. Local Scalars ..
065:       DOUBLE PRECISION   SFMIN 
066:       INTEGER            I, J, JP
067: *     ..
068: *     .. External Functions ..
069:       DOUBLE PRECISION   DLAMCH      
070:       INTEGER            IDAMAX
071:       EXTERNAL           DLAMCH, IDAMAX
072: *     ..
073: *     .. External Subroutines ..
074:       EXTERNAL           DGER, DSCAL, DSWAP, XERBLA
075: *     ..
076: *     .. Intrinsic Functions ..
077:       INTRINSIC          MAX, MIN
078: *     ..
079: *     .. Executable Statements ..
080: *
081: *     Test the input parameters.
082: *
083:       INFO = 0
084:       IF( M.LT.0 ) THEN
085:          INFO = -1
086:       ELSE IF( N.LT.0 ) THEN
087:          INFO = -2
088:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
089:          INFO = -4
090:       END IF
091:       IF( INFO.NE.0 ) THEN
092:          CALL XERBLA( 'DGETF2', -INFO )
093:          RETURN
094:       END IF
095: *
096: *     Quick return if possible
097: *
098:       IF( M.EQ.0 .OR. N.EQ.0 )
099:      $   RETURN
100: *
101: *     Compute machine safe minimum 
102: * 
103:       SFMIN = DLAMCH('S')  
104: *
105:       DO 10 J = 1, MIN( M, N )
106: *
107: *        Find pivot and test for singularity.
108: *
109:          JP = J - 1 + IDAMAX( M-J+1, A( J, J ), 1 )
110:          IPIV( J ) = JP
111:          IF( A( JP, J ).NE.ZERO ) THEN
112: *
113: *           Apply the interchange to columns 1:N.
114: *
115:             IF( JP.NE.J )
116:      $         CALL DSWAP( N, A( J, 1 ), LDA, A( JP, 1 ), LDA )
117: *
118: *           Compute elements J+1:M of J-th column.
119: *
120:             IF( J.LT.M ) THEN 
121:                IF( ABS(A( J, J )) .GE. SFMIN ) THEN 
122:                   CALL DSCAL( M-J, ONE / A( J, J ), A( J+1, J ), 1 ) 
123:                ELSE 
124:                  DO 20 I = 1, M-J 
125:                     A( J+I, J ) = A( J+I, J ) / A( J, J ) 
126:    20            CONTINUE 
127:                END IF 
128:             END IF 
129: *
130:          ELSE IF( INFO.EQ.0 ) THEN
131: *
132:             INFO = J
133:          END IF
134: *
135:          IF( J.LT.MIN( M, N ) ) THEN
136: *
137: *           Update trailing submatrix.
138: *
139:             CALL DGER( M-J, N-J, -ONE, A( J+1, J ), 1, A( J, J+1 ), LDA,
140:      $                 A( J+1, J+1 ), LDA )
141:          END IF
142:    10 CONTINUE
143:       RETURN
144: *
145: *     End of DGETF2
146: *
147:       END
148: