1*> \brief \b CPBTRF 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download CPBTRF + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cpbtrf.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cpbtrf.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cpbtrf.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE CPBTRF( UPLO, N, KD, AB, LDAB, INFO ) 22* 23* .. Scalar Arguments .. 24* CHARACTER UPLO 25* INTEGER INFO, KD, LDAB, N 26* .. 27* .. Array Arguments .. 28* COMPLEX AB( LDAB, * ) 29* .. 30* 31* 32*> \par Purpose: 33* ============= 34*> 35*> \verbatim 36*> 37*> CPBTRF computes the Cholesky factorization of a complex Hermitian 38*> positive definite band matrix A. 39*> 40*> The factorization has the form 41*> A = U**H * U, if UPLO = 'U', or 42*> A = L * L**H, 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 COMPLEX array, dimension (LDAB,N) 72*> On entry, the upper or lower triangle of the Hermitian 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**H*U or A = L*L**H 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 complexOTHERcomputational 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 CPBTRF( 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 COMPLEX AB( LDAB, * ) 153* .. 154* 155* ===================================================================== 156* 157* .. Parameters .. 158 REAL ONE, ZERO 159 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 160 COMPLEX CONE 161 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 162 INTEGER NBMAX, LDWORK 163 PARAMETER ( NBMAX = 32, LDWORK = NBMAX+1 ) 164* .. 165* .. Local Scalars .. 166 INTEGER I, I2, I3, IB, II, J, JJ, NB 167* .. 168* .. Local Arrays .. 169 COMPLEX WORK( LDWORK, NBMAX ) 170* .. 171* .. External Functions .. 172 LOGICAL LSAME 173 INTEGER ILAENV 174 EXTERNAL LSAME, ILAENV 175* .. 176* .. External Subroutines .. 177 EXTERNAL CGEMM, CHERK, CPBTF2, CPOTF2, CTRSM, XERBLA 178* .. 179* .. Intrinsic Functions .. 180 INTRINSIC MIN 181* .. 182* .. Executable Statements .. 183* 184* Test the input parameters. 185* 186 INFO = 0 187 IF( ( .NOT.LSAME( UPLO, 'U' ) ) .AND. 188 $ ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN 189 INFO = -1 190 ELSE IF( N.LT.0 ) THEN 191 INFO = -2 192 ELSE IF( KD.LT.0 ) THEN 193 INFO = -3 194 ELSE IF( LDAB.LT.KD+1 ) THEN 195 INFO = -5 196 END IF 197 IF( INFO.NE.0 ) THEN 198 CALL XERBLA( 'CPBTRF', -INFO ) 199 RETURN 200 END IF 201* 202* Quick return if possible 203* 204 IF( N.EQ.0 ) 205 $ RETURN 206* 207* Determine the block size for this environment 208* 209 NB = ILAENV( 1, 'CPBTRF', UPLO, N, KD, -1, -1 ) 210* 211* The block size must not exceed the semi-bandwidth KD, and must not 212* exceed the limit set by the size of the local array WORK. 213* 214 NB = MIN( NB, NBMAX ) 215* 216 IF( NB.LE.1 .OR. NB.GT.KD ) THEN 217* 218* Use unblocked code 219* 220 CALL CPBTF2( UPLO, N, KD, AB, LDAB, INFO ) 221 ELSE 222* 223* Use blocked code 224* 225 IF( LSAME( UPLO, 'U' ) ) THEN 226* 227* Compute the Cholesky factorization of a Hermitian band 228* matrix, given the upper triangle of the matrix in band 229* storage. 230* 231* Zero the upper triangle of the work array. 232* 233 DO 20 J = 1, NB 234 DO 10 I = 1, J - 1 235 WORK( I, J ) = ZERO 236 10 CONTINUE 237 20 CONTINUE 238* 239* Process the band matrix one diagonal block at a time. 240* 241 DO 70 I = 1, N, NB 242 IB = MIN( NB, N-I+1 ) 243* 244* Factorize the diagonal block 245* 246 CALL CPOTF2( UPLO, IB, AB( KD+1, I ), LDAB-1, II ) 247 IF( II.NE.0 ) THEN 248 INFO = I + II - 1 249 GO TO 150 250 END IF 251 IF( I+IB.LE.N ) THEN 252* 253* Update the relevant part of the trailing submatrix. 254* If A11 denotes the diagonal block which has just been 255* factorized, then we need to update the remaining 256* blocks in the diagram: 257* 258* A11 A12 A13 259* A22 A23 260* A33 261* 262* The numbers of rows and columns in the partitioning 263* are IB, I2, I3 respectively. The blocks A12, A22 and 264* A23 are empty if IB = KD. The upper triangle of A13 265* lies outside the band. 266* 267 I2 = MIN( KD-IB, N-I-IB+1 ) 268 I3 = MIN( IB, N-I-KD+1 ) 269* 270 IF( I2.GT.0 ) THEN 271* 272* Update A12 273* 274 CALL CTRSM( 'Left', 'Upper', 'Conjugate transpose', 275 $ 'Non-unit', IB, I2, CONE, 276 $ AB( KD+1, I ), LDAB-1, 277 $ AB( KD+1-IB, I+IB ), LDAB-1 ) 278* 279* Update A22 280* 281 CALL CHERK( 'Upper', 'Conjugate transpose', I2, IB, 282 $ -ONE, AB( KD+1-IB, I+IB ), LDAB-1, ONE, 283 $ AB( KD+1, I+IB ), LDAB-1 ) 284 END IF 285* 286 IF( I3.GT.0 ) THEN 287* 288* Copy the lower triangle of A13 into the work array. 289* 290 DO 40 JJ = 1, I3 291 DO 30 II = JJ, IB 292 WORK( II, JJ ) = AB( II-JJ+1, JJ+I+KD-1 ) 293 30 CONTINUE 294 40 CONTINUE 295* 296* Update A13 (in the work array). 297* 298 CALL CTRSM( 'Left', 'Upper', 'Conjugate transpose', 299 $ 'Non-unit', IB, I3, CONE, 300 $ AB( KD+1, I ), LDAB-1, WORK, LDWORK ) 301* 302* Update A23 303* 304 IF( I2.GT.0 ) 305 $ CALL CGEMM( 'Conjugate transpose', 306 $ 'No transpose', I2, I3, IB, -CONE, 307 $ AB( KD+1-IB, I+IB ), LDAB-1, WORK, 308 $ LDWORK, CONE, AB( 1+IB, I+KD ), 309 $ LDAB-1 ) 310* 311* Update A33 312* 313 CALL CHERK( 'Upper', 'Conjugate transpose', I3, IB, 314 $ -ONE, WORK, LDWORK, ONE, 315 $ AB( KD+1, I+KD ), LDAB-1 ) 316* 317* Copy the lower triangle of A13 back into place. 318* 319 DO 60 JJ = 1, I3 320 DO 50 II = JJ, IB 321 AB( II-JJ+1, JJ+I+KD-1 ) = WORK( II, JJ ) 322 50 CONTINUE 323 60 CONTINUE 324 END IF 325 END IF 326 70 CONTINUE 327 ELSE 328* 329* Compute the Cholesky factorization of a Hermitian band 330* matrix, given the lower triangle of the matrix in band 331* storage. 332* 333* Zero the lower triangle of the work array. 334* 335 DO 90 J = 1, NB 336 DO 80 I = J + 1, NB 337 WORK( I, J ) = ZERO 338 80 CONTINUE 339 90 CONTINUE 340* 341* Process the band matrix one diagonal block at a time. 342* 343 DO 140 I = 1, N, NB 344 IB = MIN( NB, N-I+1 ) 345* 346* Factorize the diagonal block 347* 348 CALL CPOTF2( UPLO, IB, AB( 1, I ), LDAB-1, II ) 349 IF( II.NE.0 ) THEN 350 INFO = I + II - 1 351 GO TO 150 352 END IF 353 IF( I+IB.LE.N ) THEN 354* 355* Update the relevant part of the trailing submatrix. 356* If A11 denotes the diagonal block which has just been 357* factorized, then we need to update the remaining 358* blocks in the diagram: 359* 360* A11 361* A21 A22 362* A31 A32 A33 363* 364* The numbers of rows and columns in the partitioning 365* are IB, I2, I3 respectively. The blocks A21, A22 and 366* A32 are empty if IB = KD. The lower triangle of A31 367* lies outside the band. 368* 369 I2 = MIN( KD-IB, N-I-IB+1 ) 370 I3 = MIN( IB, N-I-KD+1 ) 371* 372 IF( I2.GT.0 ) THEN 373* 374* Update A21 375* 376 CALL CTRSM( 'Right', 'Lower', 377 $ 'Conjugate transpose', 'Non-unit', I2, 378 $ IB, CONE, AB( 1, I ), LDAB-1, 379 $ AB( 1+IB, I ), LDAB-1 ) 380* 381* Update A22 382* 383 CALL CHERK( 'Lower', 'No transpose', I2, IB, -ONE, 384 $ AB( 1+IB, I ), LDAB-1, ONE, 385 $ AB( 1, I+IB ), LDAB-1 ) 386 END IF 387* 388 IF( I3.GT.0 ) THEN 389* 390* Copy the upper triangle of A31 into the work array. 391* 392 DO 110 JJ = 1, IB 393 DO 100 II = 1, MIN( JJ, I3 ) 394 WORK( II, JJ ) = AB( KD+1-JJ+II, JJ+I-1 ) 395 100 CONTINUE 396 110 CONTINUE 397* 398* Update A31 (in the work array). 399* 400 CALL CTRSM( 'Right', 'Lower', 401 $ 'Conjugate transpose', 'Non-unit', I3, 402 $ IB, CONE, AB( 1, I ), LDAB-1, WORK, 403 $ LDWORK ) 404* 405* Update A32 406* 407 IF( I2.GT.0 ) 408 $ CALL CGEMM( 'No transpose', 409 $ 'Conjugate transpose', I3, I2, IB, 410 $ -CONE, WORK, LDWORK, AB( 1+IB, I ), 411 $ LDAB-1, CONE, AB( 1+KD-IB, I+IB ), 412 $ LDAB-1 ) 413* 414* Update A33 415* 416 CALL CHERK( 'Lower', 'No transpose', I3, IB, -ONE, 417 $ WORK, LDWORK, ONE, AB( 1, I+KD ), 418 $ LDAB-1 ) 419* 420* Copy the upper triangle of A31 back into place. 421* 422 DO 130 JJ = 1, IB 423 DO 120 II = 1, MIN( JJ, I3 ) 424 AB( KD+1-JJ+II, JJ+I-1 ) = WORK( II, JJ ) 425 120 CONTINUE 426 130 CONTINUE 427 END IF 428 END IF 429 140 CONTINUE 430 END IF 431 END IF 432 RETURN 433* 434 150 CONTINUE 435 RETURN 436* 437* End of CPBTRF 438* 439 END 440