1*> \brief \b SPBTRF 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download SPBTRF + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/spbtrf.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/spbtrf.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/spbtrf.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SPBTRF( UPLO, N, KD, AB, LDAB, INFO ) 22* 23* .. Scalar Arguments .. 24* CHARACTER UPLO 25* INTEGER INFO, KD, LDAB, N 26* .. 27* .. Array Arguments .. 28* REAL AB( LDAB, * ) 29* .. 30* 31* 32*> \par Purpose: 33* ============= 34*> 35*> \verbatim 36*> 37*> SPBTRF computes the Cholesky factorization of a real symmetric 38*> positive definite band matrix A. 39*> 40*> The factorization has the form 41*> A = U**T * U, if UPLO = 'U', or 42*> A = L * L**T, if UPLO = 'L', 43*> where U is an upper triangular matrix and L is lower triangular. 44*> \endverbatim 45* 46* Arguments: 47* ========== 48* 49*> \param[in] UPLO 50*> \verbatim 51*> UPLO is CHARACTER*1 52*> = 'U': Upper triangle of A is stored; 53*> = 'L': Lower triangle of A is stored. 54*> \endverbatim 55*> 56*> \param[in] N 57*> \verbatim 58*> N is INTEGER 59*> The order of the matrix A. N >= 0. 60*> \endverbatim 61*> 62*> \param[in] KD 63*> \verbatim 64*> KD is INTEGER 65*> The number of superdiagonals of the matrix A if UPLO = 'U', 66*> or the number of subdiagonals if UPLO = 'L'. KD >= 0. 67*> \endverbatim 68*> 69*> \param[in,out] AB 70*> \verbatim 71*> AB is REAL array, dimension (LDAB,N) 72*> On entry, the upper or lower triangle of the symmetric band 73*> matrix A, stored in the first KD+1 rows of the array. The 74*> j-th column of A is stored in the j-th column of the array AB 75*> as follows: 76*> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; 77*> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). 78*> 79*> On exit, if INFO = 0, the triangular factor U or L from the 80*> Cholesky factorization A = U**T*U or A = L*L**T of the band 81*> matrix A, in the same storage format as A. 82*> \endverbatim 83*> 84*> \param[in] LDAB 85*> \verbatim 86*> LDAB is INTEGER 87*> The leading dimension of the array AB. LDAB >= KD+1. 88*> \endverbatim 89*> 90*> \param[out] INFO 91*> \verbatim 92*> INFO is INTEGER 93*> = 0: successful exit 94*> < 0: if INFO = -i, the i-th argument had an illegal value 95*> > 0: if INFO = i, the leading minor of order i is not 96*> positive definite, and the factorization could not be 97*> completed. 98*> \endverbatim 99* 100* Authors: 101* ======== 102* 103*> \author Univ. of Tennessee 104*> \author Univ. of California Berkeley 105*> \author Univ. of Colorado Denver 106*> \author NAG Ltd. 107* 108*> \ingroup realOTHERcomputational 109* 110*> \par Further Details: 111* ===================== 112*> 113*> \verbatim 114*> 115*> The band storage scheme is illustrated by the following example, when 116*> N = 6, KD = 2, and UPLO = 'U': 117*> 118*> On entry: On exit: 119*> 120*> * * a13 a24 a35 a46 * * u13 u24 u35 u46 121*> * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 122*> a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 123*> 124*> Similarly, if UPLO = 'L' the format of A is as follows: 125*> 126*> On entry: On exit: 127*> 128*> a11 a22 a33 a44 a55 a66 l11 l22 l33 l44 l55 l66 129*> a21 a32 a43 a54 a65 * l21 l32 l43 l54 l65 * 130*> a31 a42 a53 a64 * * l31 l42 l53 l64 * * 131*> 132*> Array elements marked * are not used by the routine. 133*> \endverbatim 134* 135*> \par Contributors: 136* ================== 137*> 138*> Peter Mayes and Giuseppe Radicati, IBM ECSEC, Rome, March 23, 1989 139* 140* ===================================================================== 141 SUBROUTINE SPBTRF( UPLO, N, KD, AB, LDAB, INFO ) 142* 143* -- LAPACK computational routine -- 144* -- LAPACK is a software package provided by Univ. of Tennessee, -- 145* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 146* 147* .. Scalar Arguments .. 148 CHARACTER UPLO 149 INTEGER INFO, KD, LDAB, N 150* .. 151* .. Array Arguments .. 152 REAL AB( LDAB, * ) 153* .. 154* 155* ===================================================================== 156* 157* .. Parameters .. 158 REAL ONE, ZERO 159 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 160 INTEGER NBMAX, LDWORK 161 PARAMETER ( NBMAX = 32, LDWORK = NBMAX+1 ) 162* .. 163* .. Local Scalars .. 164 INTEGER I, I2, I3, IB, II, J, JJ, NB 165* .. 166* .. Local Arrays .. 167 REAL WORK( LDWORK, NBMAX ) 168* .. 169* .. External Functions .. 170 LOGICAL LSAME 171 INTEGER ILAENV 172 EXTERNAL LSAME, ILAENV 173* .. 174* .. External Subroutines .. 175 EXTERNAL SGEMM, SPBTF2, SPOTF2, SSYRK, STRSM, XERBLA 176* .. 177* .. Intrinsic Functions .. 178 INTRINSIC MIN 179* .. 180* .. Executable Statements .. 181* 182* Test the input parameters. 183* 184 INFO = 0 185 IF( ( .NOT.LSAME( UPLO, 'U' ) ) .AND. 186 $ ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN 187 INFO = -1 188 ELSE IF( N.LT.0 ) THEN 189 INFO = -2 190 ELSE IF( KD.LT.0 ) THEN 191 INFO = -3 192 ELSE IF( LDAB.LT.KD+1 ) THEN 193 INFO = -5 194 END IF 195 IF( INFO.NE.0 ) THEN 196 CALL XERBLA( 'SPBTRF', -INFO ) 197 RETURN 198 END IF 199* 200* Quick return if possible 201* 202 IF( N.EQ.0 ) 203 $ RETURN 204* 205* Determine the block size for this environment 206* 207 NB = ILAENV( 1, 'SPBTRF', UPLO, N, KD, -1, -1 ) 208* 209* The block size must not exceed the semi-bandwidth KD, and must not 210* exceed the limit set by the size of the local array WORK. 211* 212 NB = MIN( NB, NBMAX ) 213* 214 IF( NB.LE.1 .OR. NB.GT.KD ) THEN 215* 216* Use unblocked code 217* 218 CALL SPBTF2( UPLO, N, KD, AB, LDAB, INFO ) 219 ELSE 220* 221* Use blocked code 222* 223 IF( LSAME( UPLO, 'U' ) ) THEN 224* 225* Compute the Cholesky factorization of a symmetric band 226* matrix, given the upper triangle of the matrix in band 227* storage. 228* 229* Zero the upper triangle of the work array. 230* 231 DO 20 J = 1, NB 232 DO 10 I = 1, J - 1 233 WORK( I, J ) = ZERO 234 10 CONTINUE 235 20 CONTINUE 236* 237* Process the band matrix one diagonal block at a time. 238* 239 DO 70 I = 1, N, NB 240 IB = MIN( NB, N-I+1 ) 241* 242* Factorize the diagonal block 243* 244 CALL SPOTF2( UPLO, IB, AB( KD+1, I ), LDAB-1, II ) 245 IF( II.NE.0 ) THEN 246 INFO = I + II - 1 247 GO TO 150 248 END IF 249 IF( I+IB.LE.N ) THEN 250* 251* Update the relevant part of the trailing submatrix. 252* If A11 denotes the diagonal block which has just been 253* factorized, then we need to update the remaining 254* blocks in the diagram: 255* 256* A11 A12 A13 257* A22 A23 258* A33 259* 260* The numbers of rows and columns in the partitioning 261* are IB, I2, I3 respectively. The blocks A12, A22 and 262* A23 are empty if IB = KD. The upper triangle of A13 263* lies outside the band. 264* 265 I2 = MIN( KD-IB, N-I-IB+1 ) 266 I3 = MIN( IB, N-I-KD+1 ) 267* 268 IF( I2.GT.0 ) THEN 269* 270* Update A12 271* 272 CALL STRSM( 'Left', 'Upper', 'Transpose', 273 $ 'Non-unit', IB, I2, ONE, AB( KD+1, I ), 274 $ LDAB-1, AB( KD+1-IB, I+IB ), LDAB-1 ) 275* 276* Update A22 277* 278 CALL SSYRK( 'Upper', 'Transpose', I2, IB, -ONE, 279 $ AB( KD+1-IB, I+IB ), LDAB-1, ONE, 280 $ AB( KD+1, I+IB ), LDAB-1 ) 281 END IF 282* 283 IF( I3.GT.0 ) THEN 284* 285* Copy the lower triangle of A13 into the work array. 286* 287 DO 40 JJ = 1, I3 288 DO 30 II = JJ, IB 289 WORK( II, JJ ) = AB( II-JJ+1, JJ+I+KD-1 ) 290 30 CONTINUE 291 40 CONTINUE 292* 293* Update A13 (in the work array). 294* 295 CALL STRSM( 'Left', 'Upper', 'Transpose', 296 $ 'Non-unit', IB, I3, ONE, AB( KD+1, I ), 297 $ LDAB-1, WORK, LDWORK ) 298* 299* Update A23 300* 301 IF( I2.GT.0 ) 302 $ CALL SGEMM( 'Transpose', 'No Transpose', I2, I3, 303 $ IB, -ONE, AB( KD+1-IB, I+IB ), 304 $ LDAB-1, WORK, LDWORK, ONE, 305 $ AB( 1+IB, I+KD ), LDAB-1 ) 306* 307* Update A33 308* 309 CALL SSYRK( 'Upper', 'Transpose', I3, IB, -ONE, 310 $ WORK, LDWORK, ONE, AB( KD+1, I+KD ), 311 $ LDAB-1 ) 312* 313* Copy the lower triangle of A13 back into place. 314* 315 DO 60 JJ = 1, I3 316 DO 50 II = JJ, IB 317 AB( II-JJ+1, JJ+I+KD-1 ) = WORK( II, JJ ) 318 50 CONTINUE 319 60 CONTINUE 320 END IF 321 END IF 322 70 CONTINUE 323 ELSE 324* 325* Compute the Cholesky factorization of a symmetric band 326* matrix, given the lower triangle of the matrix in band 327* storage. 328* 329* Zero the lower triangle of the work array. 330* 331 DO 90 J = 1, NB 332 DO 80 I = J + 1, NB 333 WORK( I, J ) = ZERO 334 80 CONTINUE 335 90 CONTINUE 336* 337* Process the band matrix one diagonal block at a time. 338* 339 DO 140 I = 1, N, NB 340 IB = MIN( NB, N-I+1 ) 341* 342* Factorize the diagonal block 343* 344 CALL SPOTF2( UPLO, IB, AB( 1, I ), LDAB-1, II ) 345 IF( II.NE.0 ) THEN 346 INFO = I + II - 1 347 GO TO 150 348 END IF 349 IF( I+IB.LE.N ) THEN 350* 351* Update the relevant part of the trailing submatrix. 352* If A11 denotes the diagonal block which has just been 353* factorized, then we need to update the remaining 354* blocks in the diagram: 355* 356* A11 357* A21 A22 358* A31 A32 A33 359* 360* The numbers of rows and columns in the partitioning 361* are IB, I2, I3 respectively. The blocks A21, A22 and 362* A32 are empty if IB = KD. The lower triangle of A31 363* lies outside the band. 364* 365 I2 = MIN( KD-IB, N-I-IB+1 ) 366 I3 = MIN( IB, N-I-KD+1 ) 367* 368 IF( I2.GT.0 ) THEN 369* 370* Update A21 371* 372 CALL STRSM( 'Right', 'Lower', 'Transpose', 373 $ 'Non-unit', I2, IB, ONE, AB( 1, I ), 374 $ LDAB-1, AB( 1+IB, I ), LDAB-1 ) 375* 376* Update A22 377* 378 CALL SSYRK( 'Lower', 'No Transpose', I2, IB, -ONE, 379 $ AB( 1+IB, I ), LDAB-1, ONE, 380 $ AB( 1, I+IB ), LDAB-1 ) 381 END IF 382* 383 IF( I3.GT.0 ) THEN 384* 385* Copy the upper triangle of A31 into the work array. 386* 387 DO 110 JJ = 1, IB 388 DO 100 II = 1, MIN( JJ, I3 ) 389 WORK( II, JJ ) = AB( KD+1-JJ+II, JJ+I-1 ) 390 100 CONTINUE 391 110 CONTINUE 392* 393* Update A31 (in the work array). 394* 395 CALL STRSM( 'Right', 'Lower', 'Transpose', 396 $ 'Non-unit', I3, IB, ONE, AB( 1, I ), 397 $ LDAB-1, WORK, LDWORK ) 398* 399* Update A32 400* 401 IF( I2.GT.0 ) 402 $ CALL SGEMM( 'No transpose', 'Transpose', I3, I2, 403 $ IB, -ONE, WORK, LDWORK, 404 $ AB( 1+IB, I ), LDAB-1, ONE, 405 $ AB( 1+KD-IB, I+IB ), LDAB-1 ) 406* 407* Update A33 408* 409 CALL SSYRK( 'Lower', 'No Transpose', I3, IB, -ONE, 410 $ WORK, LDWORK, ONE, AB( 1, I+KD ), 411 $ LDAB-1 ) 412* 413* Copy the upper triangle of A31 back into place. 414* 415 DO 130 JJ = 1, IB 416 DO 120 II = 1, MIN( JJ, I3 ) 417 AB( KD+1-JJ+II, JJ+I-1 ) = WORK( II, JJ ) 418 120 CONTINUE 419 130 CONTINUE 420 END IF 421 END IF 422 140 CONTINUE 423 END IF 424 END IF 425 RETURN 426* 427 150 CONTINUE 428 RETURN 429* 430* End of SPBTRF 431* 432 END 433