 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SGETRF( M, N, A, LDA, IPIV, INFO ) 00002 IMPLICIT NONE 00003 * 00004 * -- LAPACK routine (version 3.X) -- 00005 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00006 * May 2008 00007 * 00008 * .. Scalar Arguments .. 00009 INTEGER INFO, LDA, M, N 00010 * .. 00011 * .. Array Arguments .. 00012 INTEGER IPIV( * ) 00013 REAL A( LDA, * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * SGETRF 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 code implements an iterative version of Sivan Toledo's recursive 00029 * LU algorithm[1]. For square matrices, this iterative versions should 00030 * be within a factor of two of the optimum number of memory transfers. 00031 * 00032 * The pattern is as follows, with the large blocks of U being updated 00033 * in one call to STRSM, and the dotted lines denoting sections that 00034 * have had all pending permutations applied: 00035 * 00036 * 1 2 3 4 5 6 7 8 00037 * +-+-+---+-------+------ 00038 * | |1| | | 00039 * |.+-+ 2 | | 00040 * | | | | | 00041 * |.|.+-+-+ 4 | 00042 * | | | |1| | 00043 * | | |.+-+ | 00044 * | | | | | | 00045 * |.|.|.|.+-+-+---+ 8 00046 * | | | | | |1| | 00047 * | | | | |.+-+ 2 | 00048 * | | | | | | | | 00049 * | | | | |.|.+-+-+ 00050 * | | | | | | | |1| 00051 * | | | | | | |.+-+ 00052 * | | | | | | | | | 00053 * |.|.|.|.|.|.|.|.+----- 00054 * | | | | | | | | | 00055 * 00056 * The 1-2-1-4-1-2-1-8-... pattern is the position of the last 1 bit in 00057 * the binary expansion of the current column. Each Schur update is 00058 * applied as soon as the necessary portion of U is available. 00059 * 00060 * [1] Toledo, S. 1997. Locality of Reference in LU Decomposition with 00061 * Partial Pivoting. SIAM J. Matrix Anal. Appl. 18, 4 (Oct. 1997), 00062 * 1065-1081. http://dx.doi.org/10.1137/S0895479896297744 00063 * 00064 * Arguments 00065 * ========= 00066 * 00067 * M (input) INTEGER 00068 * The number of rows of the matrix A. M >= 0. 00069 * 00070 * N (input) INTEGER 00071 * The number of columns of the matrix A. N >= 0. 00072 * 00073 * A (input/output) REAL array, dimension (LDA,N) 00074 * On entry, the M-by-N matrix to be factored. 00075 * On exit, the factors L and U from the factorization 00076 * A = P*L*U; the unit diagonal elements of L are not stored. 00077 * 00078 * LDA (input) INTEGER 00079 * The leading dimension of the array A. LDA >= max(1,M). 00080 * 00081 * IPIV (output) INTEGER array, dimension (min(M,N)) 00082 * The pivot indices; for 1 <= i <= min(M,N), row i of the 00083 * matrix was interchanged with row IPIV(i). 00084 * 00085 * INFO (output) INTEGER 00086 * = 0: successful exit 00087 * < 0: if INFO = -i, the i-th argument had an illegal value 00088 * > 0: if INFO = i, U(i,i) is exactly zero. The factorization 00089 * has been completed, but the factor U is exactly 00090 * singular, and division by zero will occur if it is used 00091 * to solve a system of equations. 00092 * 00093 * ===================================================================== 00094 * 00095 * .. Parameters .. 00096 REAL ONE, ZERO, NEGONE 00097 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00098 PARAMETER ( NEGONE = -1.0E+0 ) 00099 * .. 00100 * .. Local Scalars .. 00101 REAL SFMIN, TMP 00102 INTEGER I, J, JP, NSTEP, NTOPIV, NPIVED, KAHEAD 00103 INTEGER KSTART, IPIVSTART, JPIVSTART, KCOLS 00104 * .. 00105 * .. External Functions .. 00106 REAL SLAMCH 00107 INTEGER ISAMAX 00108 LOGICAL SISNAN 00109 EXTERNAL SLAMCH, ISAMAX, SISNAN 00110 * .. 00111 * .. External Subroutines .. 00112 EXTERNAL STRSM, SSCAL, XERBLA, SLASWP 00113 * .. 00114 * .. Intrinsic Functions .. 00115 INTRINSIC MAX, MIN, IAND 00116 * .. 00117 * .. Executable Statements .. 00118 * 00119 * Test the input parameters. 00120 * 00121 INFO = 0 00122 IF( M.LT.0 ) THEN 00123 INFO = -1 00124 ELSE IF( N.LT.0 ) THEN 00125 INFO = -2 00126 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00127 INFO = -4 00128 END IF 00129 IF( INFO.NE.0 ) THEN 00130 CALL XERBLA( 'SGETRF', -INFO ) 00131 RETURN 00132 END IF 00133 * 00134 * Quick return if possible 00135 * 00136 IF( M.EQ.0 .OR. N.EQ.0 ) 00137 $ RETURN 00138 * 00139 * Compute machine safe minimum 00140 * 00141 SFMIN = SLAMCH( 'S' ) 00142 * 00143 NSTEP = MIN( M, N ) 00144 DO J = 1, NSTEP 00145 KAHEAD = IAND( J, -J ) 00146 KSTART = J + 1 - KAHEAD 00147 KCOLS = MIN( KAHEAD, M-J ) 00148 * 00149 * Find pivot. 00150 * 00151 JP = J - 1 + ISAMAX( M-J+1, A( J, J ), 1 ) 00152 IPIV( J ) = JP 00153 00154 ! Permute just this column. 00155 IF (JP .NE. J) THEN 00156 TMP = A( J, J ) 00157 A( J, J ) = A( JP, J ) 00158 A( JP, J ) = TMP 00159 END IF 00160 00161 ! Apply pending permutations to L 00162 NTOPIV = 1 00163 IPIVSTART = J 00164 JPIVSTART = J - NTOPIV 00165 DO WHILE ( NTOPIV .LT. KAHEAD ) 00166 CALL SLASWP( NTOPIV, A( 1, JPIVSTART ), LDA, IPIVSTART, J, 00167 $ IPIV, 1 ) 00168 IPIVSTART = IPIVSTART - NTOPIV; 00169 NTOPIV = NTOPIV * 2; 00170 JPIVSTART = JPIVSTART - NTOPIV; 00171 END DO 00172 00173 ! Permute U block to match L 00174 CALL SLASWP( KCOLS, A( 1,J+1 ), LDA, KSTART, J, IPIV, 1 ) 00175 00176 ! Factor the current column 00177 IF( A( J, J ).NE.ZERO .AND. .NOT.SISNAN( A( J, J ) ) ) THEN 00178 IF( ABS(A( J, J )) .GE. SFMIN ) THEN 00179 CALL SSCAL( M-J, ONE / A( J, J ), A( J+1, J ), 1 ) 00180 ELSE 00181 DO I = 1, M-J 00182 A( J+I, J ) = A( J+I, J ) / A( J, J ) 00183 END DO 00184 END IF 00185 ELSE IF( A( J,J ) .EQ. ZERO .AND. INFO .EQ. 0 ) THEN 00186 INFO = J 00187 END IF 00188 00189 ! Solve for U block. 00190 CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit', KAHEAD, 00191 $ KCOLS, ONE, A( KSTART, KSTART ), LDA, 00192 $ A( KSTART, J+1 ), LDA ) 00193 ! Schur complement. 00194 CALL SGEMM( 'No transpose', 'No transpose', M-J, 00195 $ KCOLS, KAHEAD, NEGONE, A( J+1, KSTART ), LDA, 00196 $ A( KSTART, J+1 ), LDA, ONE, A( J+1, J+1 ), LDA ) 00197 END DO 00198 00199 ! Handle pivot permutations on the way out of the recursion 00200 NPIVED = IAND( NSTEP, -NSTEP ) 00201 J = NSTEP - NPIVED 00202 DO WHILE ( J .GT. 0 ) 00203 NTOPIV = IAND( J, -J ) 00204 CALL SLASWP( NTOPIV, A( 1, J-NTOPIV+1 ), LDA, J+1, NSTEP, 00205 $ IPIV, 1 ) 00206 J = J - NTOPIV 00207 END DO 00208 00209 ! If short and wide, handle the rest of the columns. 00210 IF ( M .LT. N ) THEN 00211 CALL SLASWP( N-M, A( 1, M+KCOLS+1 ), LDA, 1, M, IPIV, 1 ) 00212 CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit', M, 00213 $ N-M, ONE, A, LDA, A( 1,M+KCOLS+1 ), LDA ) 00214 END IF 00215 00216 RETURN 00217 * 00218 * End of SGETRF 00219 * 00220 END
1.7.3 
