|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CGBTRF( M, N, KL, KU, AB, LDAB, 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, KL, KU, LDAB, M, N 00010 * .. 00011 * .. Array Arguments .. 00012 INTEGER IPIV( * ) 00013 COMPLEX AB( LDAB, * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * CGBTRF computes an LU factorization of a complex m-by-n band matrix A 00020 * using partial pivoting with row interchanges. 00021 * 00022 * This is the blocked version of the algorithm, calling Level 3 BLAS. 00023 * 00024 * Arguments 00025 * ========= 00026 * 00027 * M (input) INTEGER 00028 * The number of rows of the matrix A. M >= 0. 00029 * 00030 * N (input) INTEGER 00031 * The number of columns of the matrix A. N >= 0. 00032 * 00033 * KL (input) INTEGER 00034 * The number of subdiagonals within the band of A. KL >= 0. 00035 * 00036 * KU (input) INTEGER 00037 * The number of superdiagonals within the band of A. KU >= 0. 00038 * 00039 * AB (input/output) COMPLEX array, dimension (LDAB,N) 00040 * On entry, the matrix A in band storage, in rows KL+1 to 00041 * 2*KL+KU+1; rows 1 to KL of the array need not be set. 00042 * The j-th column of A is stored in the j-th column of the 00043 * array AB as follows: 00044 * AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl) 00045 * 00046 * On exit, details of the factorization: U is stored as an 00047 * upper triangular band matrix with KL+KU superdiagonals in 00048 * rows 1 to KL+KU+1, and the multipliers used during the 00049 * factorization are stored in rows KL+KU+2 to 2*KL+KU+1. 00050 * See below for further details. 00051 * 00052 * LDAB (input) INTEGER 00053 * The leading dimension of the array AB. LDAB >= 2*KL+KU+1. 00054 * 00055 * IPIV (output) INTEGER array, dimension (min(M,N)) 00056 * The pivot indices; for 1 <= i <= min(M,N), row i of the 00057 * matrix was interchanged with row IPIV(i). 00058 * 00059 * INFO (output) INTEGER 00060 * = 0: successful exit 00061 * < 0: if INFO = -i, the i-th argument had an illegal value 00062 * > 0: if INFO = +i, U(i,i) is exactly zero. The factorization 00063 * has been completed, but the factor U is exactly 00064 * singular, and division by zero will occur if it is used 00065 * to solve a system of equations. 00066 * 00067 * Further Details 00068 * =============== 00069 * 00070 * The band storage scheme is illustrated by the following example, when 00071 * M = N = 6, KL = 2, KU = 1: 00072 * 00073 * On entry: On exit: 00074 * 00075 * * * * + + + * * * u14 u25 u36 00076 * * * + + + + * * u13 u24 u35 u46 00077 * * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 00078 * a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 00079 * a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 * 00080 * a31 a42 a53 a64 * * m31 m42 m53 m64 * * 00081 * 00082 * Array elements marked * are not used by the routine; elements marked 00083 * + need not be set on entry, but are required by the routine to store 00084 * elements of U because of fill-in resulting from the row interchanges. 00085 * 00086 * ===================================================================== 00087 * 00088 * .. Parameters .. 00089 COMPLEX ONE, ZERO 00090 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ), 00091 $ ZERO = ( 0.0E+0, 0.0E+0 ) ) 00092 INTEGER NBMAX, LDWORK 00093 PARAMETER ( NBMAX = 64, LDWORK = NBMAX+1 ) 00094 * .. 00095 * .. Local Scalars .. 00096 INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP, 00097 $ JU, K2, KM, KV, NB, NW 00098 COMPLEX TEMP 00099 * .. 00100 * .. Local Arrays .. 00101 COMPLEX WORK13( LDWORK, NBMAX ), 00102 $ WORK31( LDWORK, NBMAX ) 00103 * .. 00104 * .. External Functions .. 00105 INTEGER ICAMAX, ILAENV 00106 EXTERNAL ICAMAX, ILAENV 00107 * .. 00108 * .. External Subroutines .. 00109 EXTERNAL CCOPY, CGBTF2, CGEMM, CGERU, CLASWP, CSCAL, 00110 $ CSWAP, CTRSM, XERBLA 00111 * .. 00112 * .. Intrinsic Functions .. 00113 INTRINSIC MAX, MIN 00114 * .. 00115 * .. Executable Statements .. 00116 * 00117 * KV is the number of superdiagonals in the factor U, allowing for 00118 * fill-in 00119 * 00120 KV = KU + KL 00121 * 00122 * Test the input parameters. 00123 * 00124 INFO = 0 00125 IF( M.LT.0 ) THEN 00126 INFO = -1 00127 ELSE IF( N.LT.0 ) THEN 00128 INFO = -2 00129 ELSE IF( KL.LT.0 ) THEN 00130 INFO = -3 00131 ELSE IF( KU.LT.0 ) THEN 00132 INFO = -4 00133 ELSE IF( LDAB.LT.KL+KV+1 ) THEN 00134 INFO = -6 00135 END IF 00136 IF( INFO.NE.0 ) THEN 00137 CALL XERBLA( 'CGBTRF', -INFO ) 00138 RETURN 00139 END IF 00140 * 00141 * Quick return if possible 00142 * 00143 IF( M.EQ.0 .OR. N.EQ.0 ) 00144 $ RETURN 00145 * 00146 * Determine the block size for this environment 00147 * 00148 NB = ILAENV( 1, 'CGBTRF', ' ', M, N, KL, KU ) 00149 * 00150 * The block size must not exceed the limit set by the size of the 00151 * local arrays WORK13 and WORK31. 00152 * 00153 NB = MIN( NB, NBMAX ) 00154 * 00155 IF( NB.LE.1 .OR. NB.GT.KL ) THEN 00156 * 00157 * Use unblocked code 00158 * 00159 CALL CGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) 00160 ELSE 00161 * 00162 * Use blocked code 00163 * 00164 * Zero the superdiagonal elements of the work array WORK13 00165 * 00166 DO 20 J = 1, NB 00167 DO 10 I = 1, J - 1 00168 WORK13( I, J ) = ZERO 00169 10 CONTINUE 00170 20 CONTINUE 00171 * 00172 * Zero the subdiagonal elements of the work array WORK31 00173 * 00174 DO 40 J = 1, NB 00175 DO 30 I = J + 1, NB 00176 WORK31( I, J ) = ZERO 00177 30 CONTINUE 00178 40 CONTINUE 00179 * 00180 * Gaussian elimination with partial pivoting 00181 * 00182 * Set fill-in elements in columns KU+2 to KV to zero 00183 * 00184 DO 60 J = KU + 2, MIN( KV, N ) 00185 DO 50 I = KV - J + 2, KL 00186 AB( I, J ) = ZERO 00187 50 CONTINUE 00188 60 CONTINUE 00189 * 00190 * JU is the index of the last column affected by the current 00191 * stage of the factorization 00192 * 00193 JU = 1 00194 * 00195 DO 180 J = 1, MIN( M, N ), NB 00196 JB = MIN( NB, MIN( M, N )-J+1 ) 00197 * 00198 * The active part of the matrix is partitioned 00199 * 00200 * A11 A12 A13 00201 * A21 A22 A23 00202 * A31 A32 A33 00203 * 00204 * Here A11, A21 and A31 denote the current block of JB columns 00205 * which is about to be factorized. The number of rows in the 00206 * partitioning are JB, I2, I3 respectively, and the numbers 00207 * of columns are JB, J2, J3. The superdiagonal elements of A13 00208 * and the subdiagonal elements of A31 lie outside the band. 00209 * 00210 I2 = MIN( KL-JB, M-J-JB+1 ) 00211 I3 = MIN( JB, M-J-KL+1 ) 00212 * 00213 * J2 and J3 are computed after JU has been updated. 00214 * 00215 * Factorize the current block of JB columns 00216 * 00217 DO 80 JJ = J, J + JB - 1 00218 * 00219 * Set fill-in elements in column JJ+KV to zero 00220 * 00221 IF( JJ+KV.LE.N ) THEN 00222 DO 70 I = 1, KL 00223 AB( I, JJ+KV ) = ZERO 00224 70 CONTINUE 00225 END IF 00226 * 00227 * Find pivot and test for singularity. KM is the number of 00228 * subdiagonal elements in the current column. 00229 * 00230 KM = MIN( KL, M-JJ ) 00231 JP = ICAMAX( KM+1, AB( KV+1, JJ ), 1 ) 00232 IPIV( JJ ) = JP + JJ - J 00233 IF( AB( KV+JP, JJ ).NE.ZERO ) THEN 00234 JU = MAX( JU, MIN( JJ+KU+JP-1, N ) ) 00235 IF( JP.NE.1 ) THEN 00236 * 00237 * Apply interchange to columns J to J+JB-1 00238 * 00239 IF( JP+JJ-1.LT.J+KL ) THEN 00240 * 00241 CALL CSWAP( JB, AB( KV+1+JJ-J, J ), LDAB-1, 00242 $ AB( KV+JP+JJ-J, J ), LDAB-1 ) 00243 ELSE 00244 * 00245 * The interchange affects columns J to JJ-1 of A31 00246 * which are stored in the work array WORK31 00247 * 00248 CALL CSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1, 00249 $ WORK31( JP+JJ-J-KL, 1 ), LDWORK ) 00250 CALL CSWAP( J+JB-JJ, AB( KV+1, JJ ), LDAB-1, 00251 $ AB( KV+JP, JJ ), LDAB-1 ) 00252 END IF 00253 END IF 00254 * 00255 * Compute multipliers 00256 * 00257 CALL CSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, JJ ), 00258 $ 1 ) 00259 * 00260 * Update trailing submatrix within the band and within 00261 * the current block. JM is the index of the last column 00262 * which needs to be updated. 00263 * 00264 JM = MIN( JU, J+JB-1 ) 00265 IF( JM.GT.JJ ) 00266 $ CALL CGERU( KM, JM-JJ, -ONE, AB( KV+2, JJ ), 1, 00267 $ AB( KV, JJ+1 ), LDAB-1, 00268 $ AB( KV+1, JJ+1 ), LDAB-1 ) 00269 ELSE 00270 * 00271 * If pivot is zero, set INFO to the index of the pivot 00272 * unless a zero pivot has already been found. 00273 * 00274 IF( INFO.EQ.0 ) 00275 $ INFO = JJ 00276 END IF 00277 * 00278 * Copy current column of A31 into the work array WORK31 00279 * 00280 NW = MIN( JJ-J+1, I3 ) 00281 IF( NW.GT.0 ) 00282 $ CALL CCOPY( NW, AB( KV+KL+1-JJ+J, JJ ), 1, 00283 $ WORK31( 1, JJ-J+1 ), 1 ) 00284 80 CONTINUE 00285 IF( J+JB.LE.N ) THEN 00286 * 00287 * Apply the row interchanges to the other blocks. 00288 * 00289 J2 = MIN( JU-J+1, KV ) - JB 00290 J3 = MAX( 0, JU-J-KV+1 ) 00291 * 00292 * Use CLASWP to apply the row interchanges to A12, A22, and 00293 * A32. 00294 * 00295 CALL CLASWP( J2, AB( KV+1-JB, J+JB ), LDAB-1, 1, JB, 00296 $ IPIV( J ), 1 ) 00297 * 00298 * Adjust the pivot indices. 00299 * 00300 DO 90 I = J, J + JB - 1 00301 IPIV( I ) = IPIV( I ) + J - 1 00302 90 CONTINUE 00303 * 00304 * Apply the row interchanges to A13, A23, and A33 00305 * columnwise. 00306 * 00307 K2 = J - 1 + JB + J2 00308 DO 110 I = 1, J3 00309 JJ = K2 + I 00310 DO 100 II = J + I - 1, J + JB - 1 00311 IP = IPIV( II ) 00312 IF( IP.NE.II ) THEN 00313 TEMP = AB( KV+1+II-JJ, JJ ) 00314 AB( KV+1+II-JJ, JJ ) = AB( KV+1+IP-JJ, JJ ) 00315 AB( KV+1+IP-JJ, JJ ) = TEMP 00316 END IF 00317 100 CONTINUE 00318 110 CONTINUE 00319 * 00320 * Update the relevant part of the trailing submatrix 00321 * 00322 IF( J2.GT.0 ) THEN 00323 * 00324 * Update A12 00325 * 00326 CALL CTRSM( 'Left', 'Lower', 'No transpose', 'Unit', 00327 $ JB, J2, ONE, AB( KV+1, J ), LDAB-1, 00328 $ AB( KV+1-JB, J+JB ), LDAB-1 ) 00329 * 00330 IF( I2.GT.0 ) THEN 00331 * 00332 * Update A22 00333 * 00334 CALL CGEMM( 'No transpose', 'No transpose', I2, J2, 00335 $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1, 00336 $ AB( KV+1-JB, J+JB ), LDAB-1, ONE, 00337 $ AB( KV+1, J+JB ), LDAB-1 ) 00338 END IF 00339 * 00340 IF( I3.GT.0 ) THEN 00341 * 00342 * Update A32 00343 * 00344 CALL CGEMM( 'No transpose', 'No transpose', I3, J2, 00345 $ JB, -ONE, WORK31, LDWORK, 00346 $ AB( KV+1-JB, J+JB ), LDAB-1, ONE, 00347 $ AB( KV+KL+1-JB, J+JB ), LDAB-1 ) 00348 END IF 00349 END IF 00350 * 00351 IF( J3.GT.0 ) THEN 00352 * 00353 * Copy the lower triangle of A13 into the work array 00354 * WORK13 00355 * 00356 DO 130 JJ = 1, J3 00357 DO 120 II = JJ, JB 00358 WORK13( II, JJ ) = AB( II-JJ+1, JJ+J+KV-1 ) 00359 120 CONTINUE 00360 130 CONTINUE 00361 * 00362 * Update A13 in the work array 00363 * 00364 CALL CTRSM( 'Left', 'Lower', 'No transpose', 'Unit', 00365 $ JB, J3, ONE, AB( KV+1, J ), LDAB-1, 00366 $ WORK13, LDWORK ) 00367 * 00368 IF( I2.GT.0 ) THEN 00369 * 00370 * Update A23 00371 * 00372 CALL CGEMM( 'No transpose', 'No transpose', I2, J3, 00373 $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1, 00374 $ WORK13, LDWORK, ONE, AB( 1+JB, J+KV ), 00375 $ LDAB-1 ) 00376 END IF 00377 * 00378 IF( I3.GT.0 ) THEN 00379 * 00380 * Update A33 00381 * 00382 CALL CGEMM( 'No transpose', 'No transpose', I3, J3, 00383 $ JB, -ONE, WORK31, LDWORK, WORK13, 00384 $ LDWORK, ONE, AB( 1+KL, J+KV ), LDAB-1 ) 00385 END IF 00386 * 00387 * Copy the lower triangle of A13 back into place 00388 * 00389 DO 150 JJ = 1, J3 00390 DO 140 II = JJ, JB 00391 AB( II-JJ+1, JJ+J+KV-1 ) = WORK13( II, JJ ) 00392 140 CONTINUE 00393 150 CONTINUE 00394 END IF 00395 ELSE 00396 * 00397 * Adjust the pivot indices. 00398 * 00399 DO 160 I = J, J + JB - 1 00400 IPIV( I ) = IPIV( I ) + J - 1 00401 160 CONTINUE 00402 END IF 00403 * 00404 * Partially undo the interchanges in the current block to 00405 * restore the upper triangular form of A31 and copy the upper 00406 * triangle of A31 back into place 00407 * 00408 DO 170 JJ = J + JB - 1, J, -1 00409 JP = IPIV( JJ ) - JJ + 1 00410 IF( JP.NE.1 ) THEN 00411 * 00412 * Apply interchange to columns J to JJ-1 00413 * 00414 IF( JP+JJ-1.LT.J+KL ) THEN 00415 * 00416 * The interchange does not affect A31 00417 * 00418 CALL CSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1, 00419 $ AB( KV+JP+JJ-J, J ), LDAB-1 ) 00420 ELSE 00421 * 00422 * The interchange does affect A31 00423 * 00424 CALL CSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1, 00425 $ WORK31( JP+JJ-J-KL, 1 ), LDWORK ) 00426 END IF 00427 END IF 00428 * 00429 * Copy the current column of A31 back into place 00430 * 00431 NW = MIN( I3, JJ-J+1 ) 00432 IF( NW.GT.0 ) 00433 $ CALL CCOPY( NW, WORK31( 1, JJ-J+1 ), 1, 00434 $ AB( KV+KL+1-JJ+J, JJ ), 1 ) 00435 170 CONTINUE 00436 180 CONTINUE 00437 END IF 00438 * 00439 RETURN 00440 * 00441 * End of CGBTRF 00442 * 00443 END
1.7.3 
