00001 SUBROUTINE CGETRF( 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 COMPLEX A( LDA, * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * CGETRF 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 3 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) COMPLEX 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 = -i, the i-th argument had an illegal value 00054 * > 0: if INFO = i, U(i,i) 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 COMPLEX ONE 00063 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) ) 00064 * .. 00065 * .. Local Scalars .. 00066 INTEGER I, IINFO, J, JB, NB 00067 * .. 00068 * .. External Subroutines .. 00069 EXTERNAL CGEMM, CGETF2, CLASWP, CTRSM, XERBLA 00070 * .. 00071 * .. External Functions .. 00072 INTEGER ILAENV 00073 EXTERNAL ILAENV 00074 * .. 00075 * .. Intrinsic Functions .. 00076 INTRINSIC MAX, MIN 00077 * .. 00078 * .. Executable Statements .. 00079 * 00080 * Test the input parameters. 00081 * 00082 INFO = 0 00083 IF( M.LT.0 ) THEN 00084 INFO = -1 00085 ELSE IF( N.LT.0 ) THEN 00086 INFO = -2 00087 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00088 INFO = -4 00089 END IF 00090 IF( INFO.NE.0 ) THEN 00091 CALL XERBLA( 'CGETRF', -INFO ) 00092 RETURN 00093 END IF 00094 * 00095 * Quick return if possible 00096 * 00097 IF( M.EQ.0 .OR. N.EQ.0 ) 00098 $ RETURN 00099 * 00100 * Determine the block size for this environment. 00101 * 00102 NB = ILAENV( 1, 'CGETRF', ' ', M, N, -1, -1 ) 00103 IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN 00104 * 00105 * Use unblocked code. 00106 * 00107 CALL CGETF2( M, N, A, LDA, IPIV, INFO ) 00108 ELSE 00109 * 00110 * Use blocked code. 00111 * 00112 DO 20 J = 1, MIN( M, N ), NB 00113 JB = MIN( MIN( M, N )-J+1, NB ) 00114 * 00115 * Factor diagonal and subdiagonal blocks and test for exact 00116 * singularity. 00117 * 00118 CALL CGETF2( M-J+1, JB, A( J, J ), LDA, IPIV( J ), IINFO ) 00119 * 00120 * Adjust INFO and the pivot indices. 00121 * 00122 IF( INFO.EQ.0 .AND. IINFO.GT.0 ) 00123 $ INFO = IINFO + J - 1 00124 DO 10 I = J, MIN( M, J+JB-1 ) 00125 IPIV( I ) = J - 1 + IPIV( I ) 00126 10 CONTINUE 00127 * 00128 * Apply interchanges to columns 1:J-1. 00129 * 00130 CALL CLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 ) 00131 * 00132 IF( J+JB.LE.N ) THEN 00133 * 00134 * Apply interchanges to columns J+JB:N. 00135 * 00136 CALL CLASWP( N-J-JB+1, A( 1, J+JB ), LDA, J, J+JB-1, 00137 $ IPIV, 1 ) 00138 * 00139 * Compute block row of U. 00140 * 00141 CALL CTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB, 00142 $ N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ), 00143 $ LDA ) 00144 IF( J+JB.LE.M ) THEN 00145 * 00146 * Update trailing submatrix. 00147 * 00148 CALL CGEMM( 'No transpose', 'No transpose', M-J-JB+1, 00149 $ N-J-JB+1, JB, -ONE, A( J+JB, J ), LDA, 00150 $ A( J, J+JB ), LDA, ONE, A( J+JB, J+JB ), 00151 $ LDA ) 00152 END IF 00153 END IF 00154 20 CONTINUE 00155 END IF 00156 RETURN 00157 * 00158 * End of CGETRF 00159 * 00160 END
1.7.2