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