1*> \brief \b SGBTRF 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download SGBTRF + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgbtrf.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgbtrf.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgbtrf.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) 22* 23* .. Scalar Arguments .. 24* INTEGER INFO, KL, KU, LDAB, M, N 25* .. 26* .. Array Arguments .. 27* INTEGER IPIV( * ) 28* REAL AB( LDAB, * ) 29* .. 30* 31* 32*> \par Purpose: 33* ============= 34*> 35*> \verbatim 36*> 37*> SGBTRF computes an LU factorization of a real m-by-n band matrix A 38*> using partial pivoting with row interchanges. 39*> 40*> This is the blocked version of the algorithm, calling Level 3 BLAS. 41*> \endverbatim 42* 43* Arguments: 44* ========== 45* 46*> \param[in] M 47*> \verbatim 48*> M is INTEGER 49*> The number of rows of the matrix A. M >= 0. 50*> \endverbatim 51*> 52*> \param[in] N 53*> \verbatim 54*> N is INTEGER 55*> The number of columns of the matrix A. N >= 0. 56*> \endverbatim 57*> 58*> \param[in] KL 59*> \verbatim 60*> KL is INTEGER 61*> The number of subdiagonals within the band of A. KL >= 0. 62*> \endverbatim 63*> 64*> \param[in] KU 65*> \verbatim 66*> KU is INTEGER 67*> The number of superdiagonals within the band of A. KU >= 0. 68*> \endverbatim 69*> 70*> \param[in,out] AB 71*> \verbatim 72*> AB is REAL array, dimension (LDAB,N) 73*> On entry, the matrix A in band storage, in rows KL+1 to 74*> 2*KL+KU+1; rows 1 to KL of the array need not be set. 75*> The j-th column of A is stored in the j-th column of the 76*> array AB as follows: 77*> AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl) 78*> 79*> On exit, details of the factorization: U is stored as an 80*> upper triangular band matrix with KL+KU superdiagonals in 81*> rows 1 to KL+KU+1, and the multipliers used during the 82*> factorization are stored in rows KL+KU+2 to 2*KL+KU+1. 83*> See below for further details. 84*> \endverbatim 85*> 86*> \param[in] LDAB 87*> \verbatim 88*> LDAB is INTEGER 89*> The leading dimension of the array AB. LDAB >= 2*KL+KU+1. 90*> \endverbatim 91*> 92*> \param[out] IPIV 93*> \verbatim 94*> IPIV is INTEGER array, dimension (min(M,N)) 95*> The pivot indices; for 1 <= i <= min(M,N), row i of the 96*> matrix was interchanged with row IPIV(i). 97*> \endverbatim 98*> 99*> \param[out] INFO 100*> \verbatim 101*> INFO is INTEGER 102*> = 0: successful exit 103*> < 0: if INFO = -i, the i-th argument had an illegal value 104*> > 0: if INFO = +i, U(i,i) is exactly zero. The factorization 105*> has been completed, but the factor U is exactly 106*> singular, and division by zero will occur if it is used 107*> to solve a system of equations. 108*> \endverbatim 109* 110* Authors: 111* ======== 112* 113*> \author Univ. of Tennessee 114*> \author Univ. of California Berkeley 115*> \author Univ. of Colorado Denver 116*> \author NAG Ltd. 117* 118*> \ingroup realGBcomputational 119* 120*> \par Further Details: 121* ===================== 122*> 123*> \verbatim 124*> 125*> The band storage scheme is illustrated by the following example, when 126*> M = N = 6, KL = 2, KU = 1: 127*> 128*> On entry: On exit: 129*> 130*> * * * + + + * * * u14 u25 u36 131*> * * + + + + * * u13 u24 u35 u46 132*> * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 133*> a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 134*> a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 * 135*> a31 a42 a53 a64 * * m31 m42 m53 m64 * * 136*> 137*> Array elements marked * are not used by the routine; elements marked 138*> + need not be set on entry, but are required by the routine to store 139*> elements of U because of fill-in resulting from the row interchanges. 140*> \endverbatim 141*> 142* ===================================================================== 143 SUBROUTINE SGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) 144* 145* -- LAPACK computational routine -- 146* -- LAPACK is a software package provided by Univ. of Tennessee, -- 147* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 148* 149* .. Scalar Arguments .. 150 INTEGER INFO, KL, KU, LDAB, M, N 151* .. 152* .. Array Arguments .. 153 INTEGER IPIV( * ) 154 REAL AB( LDAB, * ) 155* .. 156* 157* ===================================================================== 158* 159* .. Parameters .. 160 REAL ONE, ZERO 161 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 162 INTEGER NBMAX, LDWORK 163 PARAMETER ( NBMAX = 64, LDWORK = NBMAX+1 ) 164* .. 165* .. Local Scalars .. 166 INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP, 167 $ JU, K2, KM, KV, NB, NW 168 REAL TEMP 169* .. 170* .. Local Arrays .. 171 REAL WORK13( LDWORK, NBMAX ), 172 $ WORK31( LDWORK, NBMAX ) 173* .. 174* .. External Functions .. 175 INTEGER ILAENV, ISAMAX 176 EXTERNAL ILAENV, ISAMAX 177* .. 178* .. External Subroutines .. 179 EXTERNAL SCOPY, SGBTF2, SGEMM, SGER, SLASWP, SSCAL, 180 $ SSWAP, STRSM, XERBLA 181* .. 182* .. Intrinsic Functions .. 183 INTRINSIC MAX, MIN 184* .. 185* .. Executable Statements .. 186* 187* KV is the number of superdiagonals in the factor U, allowing for 188* fill-in 189* 190 KV = KU + KL 191* 192* Test the input parameters. 193* 194 INFO = 0 195 IF( M.LT.0 ) THEN 196 INFO = -1 197 ELSE IF( N.LT.0 ) THEN 198 INFO = -2 199 ELSE IF( KL.LT.0 ) THEN 200 INFO = -3 201 ELSE IF( KU.LT.0 ) THEN 202 INFO = -4 203 ELSE IF( LDAB.LT.KL+KV+1 ) THEN 204 INFO = -6 205 END IF 206 IF( INFO.NE.0 ) THEN 207 CALL XERBLA( 'SGBTRF', -INFO ) 208 RETURN 209 END IF 210* 211* Quick return if possible 212* 213 IF( M.EQ.0 .OR. N.EQ.0 ) 214 $ RETURN 215* 216* Determine the block size for this environment 217* 218 NB = ILAENV( 1, 'SGBTRF', ' ', M, N, KL, KU ) 219* 220* The block size must not exceed the limit set by the size of the 221* local arrays WORK13 and WORK31. 222* 223 NB = MIN( NB, NBMAX ) 224* 225 IF( NB.LE.1 .OR. NB.GT.KL ) THEN 226* 227* Use unblocked code 228* 229 CALL SGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) 230 ELSE 231* 232* Use blocked code 233* 234* Zero the superdiagonal elements of the work array WORK13 235* 236 DO 20 J = 1, NB 237 DO 10 I = 1, J - 1 238 WORK13( I, J ) = ZERO 239 10 CONTINUE 240 20 CONTINUE 241* 242* Zero the subdiagonal elements of the work array WORK31 243* 244 DO 40 J = 1, NB 245 DO 30 I = J + 1, NB 246 WORK31( I, J ) = ZERO 247 30 CONTINUE 248 40 CONTINUE 249* 250* Gaussian elimination with partial pivoting 251* 252* Set fill-in elements in columns KU+2 to KV to zero 253* 254 DO 60 J = KU + 2, MIN( KV, N ) 255 DO 50 I = KV - J + 2, KL 256 AB( I, J ) = ZERO 257 50 CONTINUE 258 60 CONTINUE 259* 260* JU is the index of the last column affected by the current 261* stage of the factorization 262* 263 JU = 1 264* 265 DO 180 J = 1, MIN( M, N ), NB 266 JB = MIN( NB, MIN( M, N )-J+1 ) 267* 268* The active part of the matrix is partitioned 269* 270* A11 A12 A13 271* A21 A22 A23 272* A31 A32 A33 273* 274* Here A11, A21 and A31 denote the current block of JB columns 275* which is about to be factorized. The number of rows in the 276* partitioning are JB, I2, I3 respectively, and the numbers 277* of columns are JB, J2, J3. The superdiagonal elements of A13 278* and the subdiagonal elements of A31 lie outside the band. 279* 280 I2 = MIN( KL-JB, M-J-JB+1 ) 281 I3 = MIN( JB, M-J-KL+1 ) 282* 283* J2 and J3 are computed after JU has been updated. 284* 285* Factorize the current block of JB columns 286* 287 DO 80 JJ = J, J + JB - 1 288* 289* Set fill-in elements in column JJ+KV to zero 290* 291 IF( JJ+KV.LE.N ) THEN 292 DO 70 I = 1, KL 293 AB( I, JJ+KV ) = ZERO 294 70 CONTINUE 295 END IF 296* 297* Find pivot and test for singularity. KM is the number of 298* subdiagonal elements in the current column. 299* 300 KM = MIN( KL, M-JJ ) 301 JP = ISAMAX( KM+1, AB( KV+1, JJ ), 1 ) 302 IPIV( JJ ) = JP + JJ - J 303 IF( AB( KV+JP, JJ ).NE.ZERO ) THEN 304 JU = MAX( JU, MIN( JJ+KU+JP-1, N ) ) 305 IF( JP.NE.1 ) THEN 306* 307* Apply interchange to columns J to J+JB-1 308* 309 IF( JP+JJ-1.LT.J+KL ) THEN 310* 311 CALL SSWAP( JB, AB( KV+1+JJ-J, J ), LDAB-1, 312 $ AB( KV+JP+JJ-J, J ), LDAB-1 ) 313 ELSE 314* 315* The interchange affects columns J to JJ-1 of A31 316* which are stored in the work array WORK31 317* 318 CALL SSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1, 319 $ WORK31( JP+JJ-J-KL, 1 ), LDWORK ) 320 CALL SSWAP( J+JB-JJ, AB( KV+1, JJ ), LDAB-1, 321 $ AB( KV+JP, JJ ), LDAB-1 ) 322 END IF 323 END IF 324* 325* Compute multipliers 326* 327 CALL SSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, JJ ), 328 $ 1 ) 329* 330* Update trailing submatrix within the band and within 331* the current block. JM is the index of the last column 332* which needs to be updated. 333* 334 JM = MIN( JU, J+JB-1 ) 335 IF( JM.GT.JJ ) 336 $ CALL SGER( KM, JM-JJ, -ONE, AB( KV+2, JJ ), 1, 337 $ AB( KV, JJ+1 ), LDAB-1, 338 $ AB( KV+1, JJ+1 ), LDAB-1 ) 339 ELSE 340* 341* If pivot is zero, set INFO to the index of the pivot 342* unless a zero pivot has already been found. 343* 344 IF( INFO.EQ.0 ) 345 $ INFO = JJ 346 END IF 347* 348* Copy current column of A31 into the work array WORK31 349* 350 NW = MIN( JJ-J+1, I3 ) 351 IF( NW.GT.0 ) 352 $ CALL SCOPY( NW, AB( KV+KL+1-JJ+J, JJ ), 1, 353 $ WORK31( 1, JJ-J+1 ), 1 ) 354 80 CONTINUE 355 IF( J+JB.LE.N ) THEN 356* 357* Apply the row interchanges to the other blocks. 358* 359 J2 = MIN( JU-J+1, KV ) - JB 360 J3 = MAX( 0, JU-J-KV+1 ) 361* 362* Use SLASWP to apply the row interchanges to A12, A22, and 363* A32. 364* 365 CALL SLASWP( J2, AB( KV+1-JB, J+JB ), LDAB-1, 1, JB, 366 $ IPIV( J ), 1 ) 367* 368* Adjust the pivot indices. 369* 370 DO 90 I = J, J + JB - 1 371 IPIV( I ) = IPIV( I ) + J - 1 372 90 CONTINUE 373* 374* Apply the row interchanges to A13, A23, and A33 375* columnwise. 376* 377 K2 = J - 1 + JB + J2 378 DO 110 I = 1, J3 379 JJ = K2 + I 380 DO 100 II = J + I - 1, J + JB - 1 381 IP = IPIV( II ) 382 IF( IP.NE.II ) THEN 383 TEMP = AB( KV+1+II-JJ, JJ ) 384 AB( KV+1+II-JJ, JJ ) = AB( KV+1+IP-JJ, JJ ) 385 AB( KV+1+IP-JJ, JJ ) = TEMP 386 END IF 387 100 CONTINUE 388 110 CONTINUE 389* 390* Update the relevant part of the trailing submatrix 391* 392 IF( J2.GT.0 ) THEN 393* 394* Update A12 395* 396 CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit', 397 $ JB, J2, ONE, AB( KV+1, J ), LDAB-1, 398 $ AB( KV+1-JB, J+JB ), LDAB-1 ) 399* 400 IF( I2.GT.0 ) THEN 401* 402* Update A22 403* 404 CALL SGEMM( 'No transpose', 'No transpose', I2, J2, 405 $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1, 406 $ AB( KV+1-JB, J+JB ), LDAB-1, ONE, 407 $ AB( KV+1, J+JB ), LDAB-1 ) 408 END IF 409* 410 IF( I3.GT.0 ) THEN 411* 412* Update A32 413* 414 CALL SGEMM( 'No transpose', 'No transpose', I3, J2, 415 $ JB, -ONE, WORK31, LDWORK, 416 $ AB( KV+1-JB, J+JB ), LDAB-1, ONE, 417 $ AB( KV+KL+1-JB, J+JB ), LDAB-1 ) 418 END IF 419 END IF 420* 421 IF( J3.GT.0 ) THEN 422* 423* Copy the lower triangle of A13 into the work array 424* WORK13 425* 426 DO 130 JJ = 1, J3 427 DO 120 II = JJ, JB 428 WORK13( II, JJ ) = AB( II-JJ+1, JJ+J+KV-1 ) 429 120 CONTINUE 430 130 CONTINUE 431* 432* Update A13 in the work array 433* 434 CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit', 435 $ JB, J3, ONE, AB( KV+1, J ), LDAB-1, 436 $ WORK13, LDWORK ) 437* 438 IF( I2.GT.0 ) THEN 439* 440* Update A23 441* 442 CALL SGEMM( 'No transpose', 'No transpose', I2, J3, 443 $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1, 444 $ WORK13, LDWORK, ONE, AB( 1+JB, J+KV ), 445 $ LDAB-1 ) 446 END IF 447* 448 IF( I3.GT.0 ) THEN 449* 450* Update A33 451* 452 CALL SGEMM( 'No transpose', 'No transpose', I3, J3, 453 $ JB, -ONE, WORK31, LDWORK, WORK13, 454 $ LDWORK, ONE, AB( 1+KL, J+KV ), LDAB-1 ) 455 END IF 456* 457* Copy the lower triangle of A13 back into place 458* 459 DO 150 JJ = 1, J3 460 DO 140 II = JJ, JB 461 AB( II-JJ+1, JJ+J+KV-1 ) = WORK13( II, JJ ) 462 140 CONTINUE 463 150 CONTINUE 464 END IF 465 ELSE 466* 467* Adjust the pivot indices. 468* 469 DO 160 I = J, J + JB - 1 470 IPIV( I ) = IPIV( I ) + J - 1 471 160 CONTINUE 472 END IF 473* 474* Partially undo the interchanges in the current block to 475* restore the upper triangular form of A31 and copy the upper 476* triangle of A31 back into place 477* 478 DO 170 JJ = J + JB - 1, J, -1 479 JP = IPIV( JJ ) - JJ + 1 480 IF( JP.NE.1 ) THEN 481* 482* Apply interchange to columns J to JJ-1 483* 484 IF( JP+JJ-1.LT.J+KL ) THEN 485* 486* The interchange does not affect A31 487* 488 CALL SSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1, 489 $ AB( KV+JP+JJ-J, J ), LDAB-1 ) 490 ELSE 491* 492* The interchange does affect A31 493* 494 CALL SSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1, 495 $ WORK31( JP+JJ-J-KL, 1 ), LDWORK ) 496 END IF 497 END IF 498* 499* Copy the current column of A31 back into place 500* 501 NW = MIN( I3, JJ-J+1 ) 502 IF( NW.GT.0 ) 503 $ CALL SCOPY( NW, WORK31( 1, JJ-J+1 ), 1, 504 $ AB( KV+KL+1-JJ+J, JJ ), 1 ) 505 170 CONTINUE 506 180 CONTINUE 507 END IF 508* 509 RETURN 510* 511* End of SGBTRF 512* 513 END 514