00001 SUBROUTINE DGETF2( M, N, A, LDA, IPIV, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 INTEGER INFO, LDA, M, N 00010 * .. 00011 * .. Array Arguments .. 00012 INTEGER IPIV( * ) 00013 DOUBLE PRECISION A( LDA, * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * DGETF2 computes an LU factorization of a general m-by-n matrix A 00020 * using partial pivoting with row interchanges. 00021 * 00022 * The factorization has the form 00023 * A = P * L * U 00024 * where P is a permutation matrix, L is lower triangular with unit 00025 * diagonal elements (lower trapezoidal if m > n), and U is upper 00026 * triangular (upper trapezoidal if m < n). 00027 * 00028 * This is the right-looking Level 2 BLAS version of the algorithm. 00029 * 00030 * Arguments 00031 * ========= 00032 * 00033 * M (input) INTEGER 00034 * The number of rows of the matrix A. M >= 0. 00035 * 00036 * N (input) INTEGER 00037 * The number of columns of the matrix A. N >= 0. 00038 * 00039 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 00040 * On entry, the m by n matrix to be factored. 00041 * On exit, the factors L and U from the factorization 00042 * A = P*L*U; the unit diagonal elements of L are not stored. 00043 * 00044 * LDA (input) INTEGER 00045 * The leading dimension of the array A. LDA >= max(1,M). 00046 * 00047 * IPIV (output) INTEGER array, dimension (min(M,N)) 00048 * The pivot indices; for 1 <= i <= min(M,N), row i of the 00049 * matrix was interchanged with row IPIV(i). 00050 * 00051 * INFO (output) INTEGER 00052 * = 0: successful exit 00053 * < 0: if INFO = -k, the k-th argument had an illegal value 00054 * > 0: if INFO = k, U(k,k) is exactly zero. The factorization 00055 * has been completed, but the factor U is exactly 00056 * singular, and division by zero will occur if it is used 00057 * to solve a system of equations. 00058 * 00059 * ===================================================================== 00060 * 00061 * .. Parameters .. 00062 DOUBLE PRECISION ONE, ZERO 00063 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00064 * .. 00065 * .. Local Scalars .. 00066 DOUBLE PRECISION SFMIN 00067 INTEGER I, J, JP 00068 * .. 00069 * .. External Functions .. 00070 DOUBLE PRECISION DLAMCH 00071 INTEGER IDAMAX 00072 EXTERNAL DLAMCH, IDAMAX 00073 * .. 00074 * .. External Subroutines .. 00075 EXTERNAL DGER, DSCAL, DSWAP, XERBLA 00076 * .. 00077 * .. Intrinsic Functions .. 00078 INTRINSIC MAX, MIN 00079 * .. 00080 * .. Executable Statements .. 00081 * 00082 * Test the input parameters. 00083 * 00084 INFO = 0 00085 IF( M.LT.0 ) THEN 00086 INFO = -1 00087 ELSE IF( N.LT.0 ) THEN 00088 INFO = -2 00089 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00090 INFO = -4 00091 END IF 00092 IF( INFO.NE.0 ) THEN 00093 CALL XERBLA( 'DGETF2', -INFO ) 00094 RETURN 00095 END IF 00096 * 00097 * Quick return if possible 00098 * 00099 IF( M.EQ.0 .OR. N.EQ.0 ) 00100 $ RETURN 00101 * 00102 * Compute machine safe minimum 00103 * 00104 SFMIN = DLAMCH('S') 00105 * 00106 DO 10 J = 1, MIN( M, N ) 00107 * 00108 * Find pivot and test for singularity. 00109 * 00110 JP = J - 1 + IDAMAX( M-J+1, A( J, J ), 1 ) 00111 IPIV( J ) = JP 00112 IF( A( JP, J ).NE.ZERO ) THEN 00113 * 00114 * Apply the interchange to columns 1:N. 00115 * 00116 IF( JP.NE.J ) 00117 $ CALL DSWAP( N, A( J, 1 ), LDA, A( JP, 1 ), LDA ) 00118 * 00119 * Compute elements J+1:M of J-th column. 00120 * 00121 IF( J.LT.M ) THEN 00122 IF( ABS(A( J, J )) .GE. SFMIN ) THEN 00123 CALL DSCAL( M-J, ONE / A( J, J ), A( J+1, J ), 1 ) 00124 ELSE 00125 DO 20 I = 1, M-J 00126 A( J+I, J ) = A( J+I, J ) / A( J, J ) 00127 20 CONTINUE 00128 END IF 00129 END IF 00130 * 00131 ELSE IF( INFO.EQ.0 ) THEN 00132 * 00133 INFO = J 00134 END IF 00135 * 00136 IF( J.LT.MIN( M, N ) ) THEN 00137 * 00138 * Update trailing submatrix. 00139 * 00140 CALL DGER( M-J, N-J, -ONE, A( J+1, J ), 1, A( J, J+1 ), LDA, 00141 $ A( J+1, J+1 ), LDA ) 00142 END IF 00143 10 CONTINUE 00144 RETURN 00145 * 00146 * End of DGETF2 00147 * 00148 END
1.7.2