![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SGBTRF 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SGBTRF + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgbtrf.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgbtrf.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgbtrf.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SGBTRF( 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 * REAL AB( LDAB, * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> SGBTRF computes an LU factorization of a real 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 REAL 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 realGBcomputational 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 SGBTRF( 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 REAL AB( LDAB, * ) 00158 * .. 00159 * 00160 * ===================================================================== 00161 * 00162 * .. Parameters .. 00163 REAL ONE, ZERO 00164 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00165 INTEGER NBMAX, LDWORK 00166 PARAMETER ( NBMAX = 64, LDWORK = NBMAX+1 ) 00167 * .. 00168 * .. Local Scalars .. 00169 INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP, 00170 $ JU, K2, KM, KV, NB, NW 00171 REAL TEMP 00172 * .. 00173 * .. Local Arrays .. 00174 REAL WORK13( LDWORK, NBMAX ), 00175 $ WORK31( LDWORK, NBMAX ) 00176 * .. 00177 * .. External Functions .. 00178 INTEGER ILAENV, ISAMAX 00179 EXTERNAL ILAENV, ISAMAX 00180 * .. 00181 * .. External Subroutines .. 00182 EXTERNAL SCOPY, SGBTF2, SGEMM, SGER, SLASWP, SSCAL, 00183 $ SSWAP, STRSM, XERBLA 00184 * .. 00185 * .. Intrinsic Functions .. 00186 INTRINSIC MAX, MIN 00187 * .. 00188 * .. Executable Statements .. 00189 * 00190 * KV is the number of superdiagonals in the factor U, allowing for 00191 * fill-in 00192 * 00193 KV = KU + KL 00194 * 00195 * Test the input parameters. 00196 * 00197 INFO = 0 00198 IF( M.LT.0 ) THEN 00199 INFO = -1 00200 ELSE IF( N.LT.0 ) THEN 00201 INFO = -2 00202 ELSE IF( KL.LT.0 ) THEN 00203 INFO = -3 00204 ELSE IF( KU.LT.0 ) THEN 00205 INFO = -4 00206 ELSE IF( LDAB.LT.KL+KV+1 ) THEN 00207 INFO = -6 00208 END IF 00209 IF( INFO.NE.0 ) THEN 00210 CALL XERBLA( 'SGBTRF', -INFO ) 00211 RETURN 00212 END IF 00213 * 00214 * Quick return if possible 00215 * 00216 IF( M.EQ.0 .OR. N.EQ.0 ) 00217 $ RETURN 00218 * 00219 * Determine the block size for this environment 00220 * 00221 NB = ILAENV( 1, 'SGBTRF', ' ', M, N, KL, KU ) 00222 * 00223 * The block size must not exceed the limit set by the size of the 00224 * local arrays WORK13 and WORK31. 00225 * 00226 NB = MIN( NB, NBMAX ) 00227 * 00228 IF( NB.LE.1 .OR. NB.GT.KL ) THEN 00229 * 00230 * Use unblocked code 00231 * 00232 CALL SGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) 00233 ELSE 00234 * 00235 * Use blocked code 00236 * 00237 * Zero the superdiagonal elements of the work array WORK13 00238 * 00239 DO 20 J = 1, NB 00240 DO 10 I = 1, J - 1 00241 WORK13( I, J ) = ZERO 00242 10 CONTINUE 00243 20 CONTINUE 00244 * 00245 * Zero the subdiagonal elements of the work array WORK31 00246 * 00247 DO 40 J = 1, NB 00248 DO 30 I = J + 1, NB 00249 WORK31( I, J ) = ZERO 00250 30 CONTINUE 00251 40 CONTINUE 00252 * 00253 * Gaussian elimination with partial pivoting 00254 * 00255 * Set fill-in elements in columns KU+2 to KV to zero 00256 * 00257 DO 60 J = KU + 2, MIN( KV, N ) 00258 DO 50 I = KV - J + 2, KL 00259 AB( I, J ) = ZERO 00260 50 CONTINUE 00261 60 CONTINUE 00262 * 00263 * JU is the index of the last column affected by the current 00264 * stage of the factorization 00265 * 00266 JU = 1 00267 * 00268 DO 180 J = 1, MIN( M, N ), NB 00269 JB = MIN( NB, MIN( M, N )-J+1 ) 00270 * 00271 * The active part of the matrix is partitioned 00272 * 00273 * A11 A12 A13 00274 * A21 A22 A23 00275 * A31 A32 A33 00276 * 00277 * Here A11, A21 and A31 denote the current block of JB columns 00278 * which is about to be factorized. The number of rows in the 00279 * partitioning are JB, I2, I3 respectively, and the numbers 00280 * of columns are JB, J2, J3. The superdiagonal elements of A13 00281 * and the subdiagonal elements of A31 lie outside the band. 00282 * 00283 I2 = MIN( KL-JB, M-J-JB+1 ) 00284 I3 = MIN( JB, M-J-KL+1 ) 00285 * 00286 * J2 and J3 are computed after JU has been updated. 00287 * 00288 * Factorize the current block of JB columns 00289 * 00290 DO 80 JJ = J, J + JB - 1 00291 * 00292 * Set fill-in elements in column JJ+KV to zero 00293 * 00294 IF( JJ+KV.LE.N ) THEN 00295 DO 70 I = 1, KL 00296 AB( I, JJ+KV ) = ZERO 00297 70 CONTINUE 00298 END IF 00299 * 00300 * Find pivot and test for singularity. KM is the number of 00301 * subdiagonal elements in the current column. 00302 * 00303 KM = MIN( KL, M-JJ ) 00304 JP = ISAMAX( KM+1, AB( KV+1, JJ ), 1 ) 00305 IPIV( JJ ) = JP + JJ - J 00306 IF( AB( KV+JP, JJ ).NE.ZERO ) THEN 00307 JU = MAX( JU, MIN( JJ+KU+JP-1, N ) ) 00308 IF( JP.NE.1 ) THEN 00309 * 00310 * Apply interchange to columns J to J+JB-1 00311 * 00312 IF( JP+JJ-1.LT.J+KL ) THEN 00313 * 00314 CALL SSWAP( JB, AB( KV+1+JJ-J, J ), LDAB-1, 00315 $ AB( KV+JP+JJ-J, J ), LDAB-1 ) 00316 ELSE 00317 * 00318 * The interchange affects columns J to JJ-1 of A31 00319 * which are stored in the work array WORK31 00320 * 00321 CALL SSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1, 00322 $ WORK31( JP+JJ-J-KL, 1 ), LDWORK ) 00323 CALL SSWAP( J+JB-JJ, AB( KV+1, JJ ), LDAB-1, 00324 $ AB( KV+JP, JJ ), LDAB-1 ) 00325 END IF 00326 END IF 00327 * 00328 * Compute multipliers 00329 * 00330 CALL SSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, JJ ), 00331 $ 1 ) 00332 * 00333 * Update trailing submatrix within the band and within 00334 * the current block. JM is the index of the last column 00335 * which needs to be updated. 00336 * 00337 JM = MIN( JU, J+JB-1 ) 00338 IF( JM.GT.JJ ) 00339 $ CALL SGER( KM, JM-JJ, -ONE, AB( KV+2, JJ ), 1, 00340 $ AB( KV, JJ+1 ), LDAB-1, 00341 $ AB( KV+1, JJ+1 ), LDAB-1 ) 00342 ELSE 00343 * 00344 * If pivot is zero, set INFO to the index of the pivot 00345 * unless a zero pivot has already been found. 00346 * 00347 IF( INFO.EQ.0 ) 00348 $ INFO = JJ 00349 END IF 00350 * 00351 * Copy current column of A31 into the work array WORK31 00352 * 00353 NW = MIN( JJ-J+1, I3 ) 00354 IF( NW.GT.0 ) 00355 $ CALL SCOPY( NW, AB( KV+KL+1-JJ+J, JJ ), 1, 00356 $ WORK31( 1, JJ-J+1 ), 1 ) 00357 80 CONTINUE 00358 IF( J+JB.LE.N ) THEN 00359 * 00360 * Apply the row interchanges to the other blocks. 00361 * 00362 J2 = MIN( JU-J+1, KV ) - JB 00363 J3 = MAX( 0, JU-J-KV+1 ) 00364 * 00365 * Use SLASWP to apply the row interchanges to A12, A22, and 00366 * A32. 00367 * 00368 CALL SLASWP( J2, AB( KV+1-JB, J+JB ), LDAB-1, 1, JB, 00369 $ IPIV( J ), 1 ) 00370 * 00371 * Adjust the pivot indices. 00372 * 00373 DO 90 I = J, J + JB - 1 00374 IPIV( I ) = IPIV( I ) + J - 1 00375 90 CONTINUE 00376 * 00377 * Apply the row interchanges to A13, A23, and A33 00378 * columnwise. 00379 * 00380 K2 = J - 1 + JB + J2 00381 DO 110 I = 1, J3 00382 JJ = K2 + I 00383 DO 100 II = J + I - 1, J + JB - 1 00384 IP = IPIV( II ) 00385 IF( IP.NE.II ) THEN 00386 TEMP = AB( KV+1+II-JJ, JJ ) 00387 AB( KV+1+II-JJ, JJ ) = AB( KV+1+IP-JJ, JJ ) 00388 AB( KV+1+IP-JJ, JJ ) = TEMP 00389 END IF 00390 100 CONTINUE 00391 110 CONTINUE 00392 * 00393 * Update the relevant part of the trailing submatrix 00394 * 00395 IF( J2.GT.0 ) THEN 00396 * 00397 * Update A12 00398 * 00399 CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit', 00400 $ JB, J2, ONE, AB( KV+1, J ), LDAB-1, 00401 $ AB( KV+1-JB, J+JB ), LDAB-1 ) 00402 * 00403 IF( I2.GT.0 ) THEN 00404 * 00405 * Update A22 00406 * 00407 CALL SGEMM( 'No transpose', 'No transpose', I2, J2, 00408 $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1, 00409 $ AB( KV+1-JB, J+JB ), LDAB-1, ONE, 00410 $ AB( KV+1, J+JB ), LDAB-1 ) 00411 END IF 00412 * 00413 IF( I3.GT.0 ) THEN 00414 * 00415 * Update A32 00416 * 00417 CALL SGEMM( 'No transpose', 'No transpose', I3, J2, 00418 $ JB, -ONE, WORK31, LDWORK, 00419 $ AB( KV+1-JB, J+JB ), LDAB-1, ONE, 00420 $ AB( KV+KL+1-JB, J+JB ), LDAB-1 ) 00421 END IF 00422 END IF 00423 * 00424 IF( J3.GT.0 ) THEN 00425 * 00426 * Copy the lower triangle of A13 into the work array 00427 * WORK13 00428 * 00429 DO 130 JJ = 1, J3 00430 DO 120 II = JJ, JB 00431 WORK13( II, JJ ) = AB( II-JJ+1, JJ+J+KV-1 ) 00432 120 CONTINUE 00433 130 CONTINUE 00434 * 00435 * Update A13 in the work array 00436 * 00437 CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit', 00438 $ JB, J3, ONE, AB( KV+1, J ), LDAB-1, 00439 $ WORK13, LDWORK ) 00440 * 00441 IF( I2.GT.0 ) THEN 00442 * 00443 * Update A23 00444 * 00445 CALL SGEMM( 'No transpose', 'No transpose', I2, J3, 00446 $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1, 00447 $ WORK13, LDWORK, ONE, AB( 1+JB, J+KV ), 00448 $ LDAB-1 ) 00449 END IF 00450 * 00451 IF( I3.GT.0 ) THEN 00452 * 00453 * Update A33 00454 * 00455 CALL SGEMM( 'No transpose', 'No transpose', I3, J3, 00456 $ JB, -ONE, WORK31, LDWORK, WORK13, 00457 $ LDWORK, ONE, AB( 1+KL, J+KV ), LDAB-1 ) 00458 END IF 00459 * 00460 * Copy the lower triangle of A13 back into place 00461 * 00462 DO 150 JJ = 1, J3 00463 DO 140 II = JJ, JB 00464 AB( II-JJ+1, JJ+J+KV-1 ) = WORK13( II, JJ ) 00465 140 CONTINUE 00466 150 CONTINUE 00467 END IF 00468 ELSE 00469 * 00470 * Adjust the pivot indices. 00471 * 00472 DO 160 I = J, J + JB - 1 00473 IPIV( I ) = IPIV( I ) + J - 1 00474 160 CONTINUE 00475 END IF 00476 * 00477 * Partially undo the interchanges in the current block to 00478 * restore the upper triangular form of A31 and copy the upper 00479 * triangle of A31 back into place 00480 * 00481 DO 170 JJ = J + JB - 1, J, -1 00482 JP = IPIV( JJ ) - JJ + 1 00483 IF( JP.NE.1 ) THEN 00484 * 00485 * Apply interchange to columns J to JJ-1 00486 * 00487 IF( JP+JJ-1.LT.J+KL ) THEN 00488 * 00489 * The interchange does not affect A31 00490 * 00491 CALL SSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1, 00492 $ AB( KV+JP+JJ-J, J ), LDAB-1 ) 00493 ELSE 00494 * 00495 * The interchange does affect A31 00496 * 00497 CALL SSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1, 00498 $ WORK31( JP+JJ-J-KL, 1 ), LDWORK ) 00499 END IF 00500 END IF 00501 * 00502 * Copy the current column of A31 back into place 00503 * 00504 NW = MIN( I3, JJ-J+1 ) 00505 IF( NW.GT.0 ) 00506 $ CALL SCOPY( NW, WORK31( 1, JJ-J+1 ), 1, 00507 $ AB( KV+KL+1-JJ+J, JJ ), 1 ) 00508 170 CONTINUE 00509 180 CONTINUE 00510 END IF 00511 * 00512 RETURN 00513 * 00514 * End of SGBTRF 00515 * 00516 END