|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DGBTRF( 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 DOUBLE PRECISION AB( LDAB, * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * DGBTRF computes an LU factorization of a real 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) DOUBLE PRECISION 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 DOUBLE PRECISION ONE, ZERO 00090 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00091 INTEGER NBMAX, LDWORK 00092 PARAMETER ( NBMAX = 64, LDWORK = NBMAX+1 ) 00093 * .. 00094 * .. Local Scalars .. 00095 INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP, 00096 $ JU, K2, KM, KV, NB, NW 00097 DOUBLE PRECISION TEMP 00098 * .. 00099 * .. Local Arrays .. 00100 DOUBLE PRECISION WORK13( LDWORK, NBMAX ), 00101 $ WORK31( LDWORK, NBMAX ) 00102 * .. 00103 * .. External Functions .. 00104 INTEGER IDAMAX, ILAENV 00105 EXTERNAL IDAMAX, ILAENV 00106 * .. 00107 * .. External Subroutines .. 00108 EXTERNAL DCOPY, DGBTF2, DGEMM, DGER, DLASWP, DSCAL, 00109 $ DSWAP, DTRSM, XERBLA 00110 * .. 00111 * .. Intrinsic Functions .. 00112 INTRINSIC MAX, MIN 00113 * .. 00114 * .. Executable Statements .. 00115 * 00116 * KV is the number of superdiagonals in the factor U, allowing for 00117 * fill-in 00118 * 00119 KV = KU + KL 00120 * 00121 * Test the input parameters. 00122 * 00123 INFO = 0 00124 IF( M.LT.0 ) THEN 00125 INFO = -1 00126 ELSE IF( N.LT.0 ) THEN 00127 INFO = -2 00128 ELSE IF( KL.LT.0 ) THEN 00129 INFO = -3 00130 ELSE IF( KU.LT.0 ) THEN 00131 INFO = -4 00132 ELSE IF( LDAB.LT.KL+KV+1 ) THEN 00133 INFO = -6 00134 END IF 00135 IF( INFO.NE.0 ) THEN 00136 CALL XERBLA( 'DGBTRF', -INFO ) 00137 RETURN 00138 END IF 00139 * 00140 * Quick return if possible 00141 * 00142 IF( M.EQ.0 .OR. N.EQ.0 ) 00143 $ RETURN 00144 * 00145 * Determine the block size for this environment 00146 * 00147 NB = ILAENV( 1, 'DGBTRF', ' ', M, N, KL, KU ) 00148 * 00149 * The block size must not exceed the limit set by the size of the 00150 * local arrays WORK13 and WORK31. 00151 * 00152 NB = MIN( NB, NBMAX ) 00153 * 00154 IF( NB.LE.1 .OR. NB.GT.KL ) THEN 00155 * 00156 * Use unblocked code 00157 * 00158 CALL DGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) 00159 ELSE 00160 * 00161 * Use blocked code 00162 * 00163 * Zero the superdiagonal elements of the work array WORK13 00164 * 00165 DO 20 J = 1, NB 00166 DO 10 I = 1, J - 1 00167 WORK13( I, J ) = ZERO 00168 10 CONTINUE 00169 20 CONTINUE 00170 * 00171 * Zero the subdiagonal elements of the work array WORK31 00172 * 00173 DO 40 J = 1, NB 00174 DO 30 I = J + 1, NB 00175 WORK31( I, J ) = ZERO 00176 30 CONTINUE 00177 40 CONTINUE 00178 * 00179 * Gaussian elimination with partial pivoting 00180 * 00181 * Set fill-in elements in columns KU+2 to KV to zero 00182 * 00183 DO 60 J = KU + 2, MIN( KV, N ) 00184 DO 50 I = KV - J + 2, KL 00185 AB( I, J ) = ZERO 00186 50 CONTINUE 00187 60 CONTINUE 00188 * 00189 * JU is the index of the last column affected by the current 00190 * stage of the factorization 00191 * 00192 JU = 1 00193 * 00194 DO 180 J = 1, MIN( M, N ), NB 00195 JB = MIN( NB, MIN( M, N )-J+1 ) 00196 * 00197 * The active part of the matrix is partitioned 00198 * 00199 * A11 A12 A13 00200 * A21 A22 A23 00201 * A31 A32 A33 00202 * 00203 * Here A11, A21 and A31 denote the current block of JB columns 00204 * which is about to be factorized. The number of rows in the 00205 * partitioning are JB, I2, I3 respectively, and the numbers 00206 * of columns are JB, J2, J3. The superdiagonal elements of A13 00207 * and the subdiagonal elements of A31 lie outside the band. 00208 * 00209 I2 = MIN( KL-JB, M-J-JB+1 ) 00210 I3 = MIN( JB, M-J-KL+1 ) 00211 * 00212 * J2 and J3 are computed after JU has been updated. 00213 * 00214 * Factorize the current block of JB columns 00215 * 00216 DO 80 JJ = J, J + JB - 1 00217 * 00218 * Set fill-in elements in column JJ+KV to zero 00219 * 00220 IF( JJ+KV.LE.N ) THEN 00221 DO 70 I = 1, KL 00222 AB( I, JJ+KV ) = ZERO 00223 70 CONTINUE 00224 END IF 00225 * 00226 * Find pivot and test for singularity. KM is the number of 00227 * subdiagonal elements in the current column. 00228 * 00229 KM = MIN( KL, M-JJ ) 00230 JP = IDAMAX( KM+1, AB( KV+1, JJ ), 1 ) 00231 IPIV( JJ ) = JP + JJ - J 00232 IF( AB( KV+JP, JJ ).NE.ZERO ) THEN 00233 JU = MAX( JU, MIN( JJ+KU+JP-1, N ) ) 00234 IF( JP.NE.1 ) THEN 00235 * 00236 * Apply interchange to columns J to J+JB-1 00237 * 00238 IF( JP+JJ-1.LT.J+KL ) THEN 00239 * 00240 CALL DSWAP( JB, AB( KV+1+JJ-J, J ), LDAB-1, 00241 $ AB( KV+JP+JJ-J, J ), LDAB-1 ) 00242 ELSE 00243 * 00244 * The interchange affects columns J to JJ-1 of A31 00245 * which are stored in the work array WORK31 00246 * 00247 CALL DSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1, 00248 $ WORK31( JP+JJ-J-KL, 1 ), LDWORK ) 00249 CALL DSWAP( J+JB-JJ, AB( KV+1, JJ ), LDAB-1, 00250 $ AB( KV+JP, JJ ), LDAB-1 ) 00251 END IF 00252 END IF 00253 * 00254 * Compute multipliers 00255 * 00256 CALL DSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, JJ ), 00257 $ 1 ) 00258 * 00259 * Update trailing submatrix within the band and within 00260 * the current block. JM is the index of the last column 00261 * which needs to be updated. 00262 * 00263 JM = MIN( JU, J+JB-1 ) 00264 IF( JM.GT.JJ ) 00265 $ CALL DGER( KM, JM-JJ, -ONE, AB( KV+2, JJ ), 1, 00266 $ AB( KV, JJ+1 ), LDAB-1, 00267 $ AB( KV+1, JJ+1 ), LDAB-1 ) 00268 ELSE 00269 * 00270 * If pivot is zero, set INFO to the index of the pivot 00271 * unless a zero pivot has already been found. 00272 * 00273 IF( INFO.EQ.0 ) 00274 $ INFO = JJ 00275 END IF 00276 * 00277 * Copy current column of A31 into the work array WORK31 00278 * 00279 NW = MIN( JJ-J+1, I3 ) 00280 IF( NW.GT.0 ) 00281 $ CALL DCOPY( NW, AB( KV+KL+1-JJ+J, JJ ), 1, 00282 $ WORK31( 1, JJ-J+1 ), 1 ) 00283 80 CONTINUE 00284 IF( J+JB.LE.N ) THEN 00285 * 00286 * Apply the row interchanges to the other blocks. 00287 * 00288 J2 = MIN( JU-J+1, KV ) - JB 00289 J3 = MAX( 0, JU-J-KV+1 ) 00290 * 00291 * Use DLASWP to apply the row interchanges to A12, A22, and 00292 * A32. 00293 * 00294 CALL DLASWP( J2, AB( KV+1-JB, J+JB ), LDAB-1, 1, JB, 00295 $ IPIV( J ), 1 ) 00296 * 00297 * Adjust the pivot indices. 00298 * 00299 DO 90 I = J, J + JB - 1 00300 IPIV( I ) = IPIV( I ) + J - 1 00301 90 CONTINUE 00302 * 00303 * Apply the row interchanges to A13, A23, and A33 00304 * columnwise. 00305 * 00306 K2 = J - 1 + JB + J2 00307 DO 110 I = 1, J3 00308 JJ = K2 + I 00309 DO 100 II = J + I - 1, J + JB - 1 00310 IP = IPIV( II ) 00311 IF( IP.NE.II ) THEN 00312 TEMP = AB( KV+1+II-JJ, JJ ) 00313 AB( KV+1+II-JJ, JJ ) = AB( KV+1+IP-JJ, JJ ) 00314 AB( KV+1+IP-JJ, JJ ) = TEMP 00315 END IF 00316 100 CONTINUE 00317 110 CONTINUE 00318 * 00319 * Update the relevant part of the trailing submatrix 00320 * 00321 IF( J2.GT.0 ) THEN 00322 * 00323 * Update A12 00324 * 00325 CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', 00326 $ JB, J2, ONE, AB( KV+1, J ), LDAB-1, 00327 $ AB( KV+1-JB, J+JB ), LDAB-1 ) 00328 * 00329 IF( I2.GT.0 ) THEN 00330 * 00331 * Update A22 00332 * 00333 CALL DGEMM( 'No transpose', 'No transpose', I2, J2, 00334 $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1, 00335 $ AB( KV+1-JB, J+JB ), LDAB-1, ONE, 00336 $ AB( KV+1, J+JB ), LDAB-1 ) 00337 END IF 00338 * 00339 IF( I3.GT.0 ) THEN 00340 * 00341 * Update A32 00342 * 00343 CALL DGEMM( 'No transpose', 'No transpose', I3, J2, 00344 $ JB, -ONE, WORK31, LDWORK, 00345 $ AB( KV+1-JB, J+JB ), LDAB-1, ONE, 00346 $ AB( KV+KL+1-JB, J+JB ), LDAB-1 ) 00347 END IF 00348 END IF 00349 * 00350 IF( J3.GT.0 ) THEN 00351 * 00352 * Copy the lower triangle of A13 into the work array 00353 * WORK13 00354 * 00355 DO 130 JJ = 1, J3 00356 DO 120 II = JJ, JB 00357 WORK13( II, JJ ) = AB( II-JJ+1, JJ+J+KV-1 ) 00358 120 CONTINUE 00359 130 CONTINUE 00360 * 00361 * Update A13 in the work array 00362 * 00363 CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', 00364 $ JB, J3, ONE, AB( KV+1, J ), LDAB-1, 00365 $ WORK13, LDWORK ) 00366 * 00367 IF( I2.GT.0 ) THEN 00368 * 00369 * Update A23 00370 * 00371 CALL DGEMM( 'No transpose', 'No transpose', I2, J3, 00372 $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1, 00373 $ WORK13, LDWORK, ONE, AB( 1+JB, J+KV ), 00374 $ LDAB-1 ) 00375 END IF 00376 * 00377 IF( I3.GT.0 ) THEN 00378 * 00379 * Update A33 00380 * 00381 CALL DGEMM( 'No transpose', 'No transpose', I3, J3, 00382 $ JB, -ONE, WORK31, LDWORK, WORK13, 00383 $ LDWORK, ONE, AB( 1+KL, J+KV ), LDAB-1 ) 00384 END IF 00385 * 00386 * Copy the lower triangle of A13 back into place 00387 * 00388 DO 150 JJ = 1, J3 00389 DO 140 II = JJ, JB 00390 AB( II-JJ+1, JJ+J+KV-1 ) = WORK13( II, JJ ) 00391 140 CONTINUE 00392 150 CONTINUE 00393 END IF 00394 ELSE 00395 * 00396 * Adjust the pivot indices. 00397 * 00398 DO 160 I = J, J + JB - 1 00399 IPIV( I ) = IPIV( I ) + J - 1 00400 160 CONTINUE 00401 END IF 00402 * 00403 * Partially undo the interchanges in the current block to 00404 * restore the upper triangular form of A31 and copy the upper 00405 * triangle of A31 back into place 00406 * 00407 DO 170 JJ = J + JB - 1, J, -1 00408 JP = IPIV( JJ ) - JJ + 1 00409 IF( JP.NE.1 ) THEN 00410 * 00411 * Apply interchange to columns J to JJ-1 00412 * 00413 IF( JP+JJ-1.LT.J+KL ) THEN 00414 * 00415 * The interchange does not affect A31 00416 * 00417 CALL DSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1, 00418 $ AB( KV+JP+JJ-J, J ), LDAB-1 ) 00419 ELSE 00420 * 00421 * The interchange does affect A31 00422 * 00423 CALL DSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1, 00424 $ WORK31( JP+JJ-J-KL, 1 ), LDWORK ) 00425 END IF 00426 END IF 00427 * 00428 * Copy the current column of A31 back into place 00429 * 00430 NW = MIN( I3, JJ-J+1 ) 00431 IF( NW.GT.0 ) 00432 $ CALL DCOPY( NW, WORK31( 1, JJ-J+1 ), 1, 00433 $ AB( KV+KL+1-JJ+J, JJ ), 1 ) 00434 170 CONTINUE 00435 180 CONTINUE 00436 END IF 00437 * 00438 RETURN 00439 * 00440 * End of DGBTRF 00441 * 00442 END
1.7.3 
