1*> \brief \b CHBTRD 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download CHBTRD + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chbtrd.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chbtrd.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chbtrd.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE CHBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, 22* WORK, INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER UPLO, VECT 26* INTEGER INFO, KD, LDAB, LDQ, N 27* .. 28* .. Array Arguments .. 29* REAL D( * ), E( * ) 30* COMPLEX AB( LDAB, * ), Q( LDQ, * ), WORK( * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> CHBTRD reduces a complex Hermitian band matrix A to real symmetric 40*> tridiagonal form T by a unitary similarity transformation: 41*> Q**H * A * Q = T. 42*> \endverbatim 43* 44* Arguments: 45* ========== 46* 47*> \param[in] VECT 48*> \verbatim 49*> VECT is CHARACTER*1 50*> = 'N': do not form Q; 51*> = 'V': form Q; 52*> = 'U': update a matrix X, by forming X*Q. 53*> \endverbatim 54*> 55*> \param[in] UPLO 56*> \verbatim 57*> UPLO is CHARACTER*1 58*> = 'U': Upper triangle of A is stored; 59*> = 'L': Lower triangle of A is stored. 60*> \endverbatim 61*> 62*> \param[in] N 63*> \verbatim 64*> N is INTEGER 65*> The order of the matrix A. N >= 0. 66*> \endverbatim 67*> 68*> \param[in] KD 69*> \verbatim 70*> KD is INTEGER 71*> The number of superdiagonals of the matrix A if UPLO = 'U', 72*> or the number of subdiagonals if UPLO = 'L'. KD >= 0. 73*> \endverbatim 74*> 75*> \param[in,out] AB 76*> \verbatim 77*> AB is COMPLEX array, dimension (LDAB,N) 78*> On entry, the upper or lower triangle of the Hermitian band 79*> matrix A, stored in the first KD+1 rows of the array. The 80*> j-th column of A is stored in the j-th column of the array AB 81*> as follows: 82*> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; 83*> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). 84*> On exit, the diagonal elements of AB are overwritten by the 85*> diagonal elements of the tridiagonal matrix T; if KD > 0, the 86*> elements on the first superdiagonal (if UPLO = 'U') or the 87*> first subdiagonal (if UPLO = 'L') are overwritten by the 88*> off-diagonal elements of T; the rest of AB is overwritten by 89*> values generated during the reduction. 90*> \endverbatim 91*> 92*> \param[in] LDAB 93*> \verbatim 94*> LDAB is INTEGER 95*> The leading dimension of the array AB. LDAB >= KD+1. 96*> \endverbatim 97*> 98*> \param[out] D 99*> \verbatim 100*> D is REAL array, dimension (N) 101*> The diagonal elements of the tridiagonal matrix T. 102*> \endverbatim 103*> 104*> \param[out] E 105*> \verbatim 106*> E is REAL array, dimension (N-1) 107*> The off-diagonal elements of the tridiagonal matrix T: 108*> E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'. 109*> \endverbatim 110*> 111*> \param[in,out] Q 112*> \verbatim 113*> Q is COMPLEX array, dimension (LDQ,N) 114*> On entry, if VECT = 'U', then Q must contain an N-by-N 115*> matrix X; if VECT = 'N' or 'V', then Q need not be set. 116*> 117*> On exit: 118*> if VECT = 'V', Q contains the N-by-N unitary matrix Q; 119*> if VECT = 'U', Q contains the product X*Q; 120*> if VECT = 'N', the array Q is not referenced. 121*> \endverbatim 122*> 123*> \param[in] LDQ 124*> \verbatim 125*> LDQ is INTEGER 126*> The leading dimension of the array Q. 127*> LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'. 128*> \endverbatim 129*> 130*> \param[out] WORK 131*> \verbatim 132*> WORK is COMPLEX array, dimension (N) 133*> \endverbatim 134*> 135*> \param[out] INFO 136*> \verbatim 137*> INFO is INTEGER 138*> = 0: successful exit 139*> < 0: if INFO = -i, the i-th argument had an illegal value 140*> \endverbatim 141* 142* Authors: 143* ======== 144* 145*> \author Univ. of Tennessee 146*> \author Univ. of California Berkeley 147*> \author Univ. of Colorado Denver 148*> \author NAG Ltd. 149* 150*> \ingroup complexOTHERcomputational 151* 152*> \par Further Details: 153* ===================== 154*> 155*> \verbatim 156*> 157*> Modified by Linda Kaufman, Bell Labs. 158*> \endverbatim 159*> 160* ===================================================================== 161 SUBROUTINE CHBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, 162 $ WORK, INFO ) 163* 164* -- LAPACK computational routine -- 165* -- LAPACK is a software package provided by Univ. of Tennessee, -- 166* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 167* 168* .. Scalar Arguments .. 169 CHARACTER UPLO, VECT 170 INTEGER INFO, KD, LDAB, LDQ, N 171* .. 172* .. Array Arguments .. 173 REAL D( * ), E( * ) 174 COMPLEX AB( LDAB, * ), Q( LDQ, * ), WORK( * ) 175* .. 176* 177* ===================================================================== 178* 179* .. Parameters .. 180 REAL ZERO 181 PARAMETER ( ZERO = 0.0E+0 ) 182 COMPLEX CZERO, CONE 183 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 184 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 185* .. 186* .. Local Scalars .. 187 LOGICAL INITQ, UPPER, WANTQ 188 INTEGER I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J, 189 $ J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1, 190 $ KDM1, KDN, L, LAST, LEND, NQ, NR, NRT 191 REAL ABST 192 COMPLEX T, TEMP 193* .. 194* .. External Subroutines .. 195 EXTERNAL CLACGV, CLAR2V, CLARGV, CLARTG, CLARTV, CLASET, 196 $ CROT, CSCAL, XERBLA 197* .. 198* .. Intrinsic Functions .. 199 INTRINSIC ABS, CONJG, MAX, MIN, REAL 200* .. 201* .. External Functions .. 202 LOGICAL LSAME 203 EXTERNAL LSAME 204* .. 205* .. Executable Statements .. 206* 207* Test the input parameters 208* 209 INITQ = LSAME( VECT, 'V' ) 210 WANTQ = INITQ .OR. LSAME( VECT, 'U' ) 211 UPPER = LSAME( UPLO, 'U' ) 212 KD1 = KD + 1 213 KDM1 = KD - 1 214 INCX = LDAB - 1 215 IQEND = 1 216* 217 INFO = 0 218 IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'N' ) ) THEN 219 INFO = -1 220 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 221 INFO = -2 222 ELSE IF( N.LT.0 ) THEN 223 INFO = -3 224 ELSE IF( KD.LT.0 ) THEN 225 INFO = -4 226 ELSE IF( LDAB.LT.KD1 ) THEN 227 INFO = -6 228 ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN 229 INFO = -10 230 END IF 231 IF( INFO.NE.0 ) THEN 232 CALL XERBLA( 'CHBTRD', -INFO ) 233 RETURN 234 END IF 235* 236* Quick return if possible 237* 238 IF( N.EQ.0 ) 239 $ RETURN 240* 241* Initialize Q to the unit matrix, if needed 242* 243 IF( INITQ ) 244 $ CALL CLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) 245* 246* Wherever possible, plane rotations are generated and applied in 247* vector operations of length NR over the index set J1:J2:KD1. 248* 249* The real cosines and complex sines of the plane rotations are 250* stored in the arrays D and WORK. 251* 252 INCA = KD1*LDAB 253 KDN = MIN( N-1, KD ) 254 IF( UPPER ) THEN 255* 256 IF( KD.GT.1 ) THEN 257* 258* Reduce to complex Hermitian tridiagonal form, working with 259* the upper triangle 260* 261 NR = 0 262 J1 = KDN + 2 263 J2 = 1 264* 265 AB( KD1, 1 ) = REAL( AB( KD1, 1 ) ) 266 DO 90 I = 1, N - 2 267* 268* Reduce i-th row of matrix to tridiagonal form 269* 270 DO 80 K = KDN + 1, 2, -1 271 J1 = J1 + KDN 272 J2 = J2 + KDN 273* 274 IF( NR.GT.0 ) THEN 275* 276* generate plane rotations to annihilate nonzero 277* elements which have been created outside the band 278* 279 CALL CLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ), 280 $ KD1, D( J1 ), KD1 ) 281* 282* apply rotations from the right 283* 284* 285* Dependent on the the number of diagonals either 286* CLARTV or CROT is used 287* 288 IF( NR.GE.2*KD-1 ) THEN 289 DO 10 L = 1, KD - 1 290 CALL CLARTV( NR, AB( L+1, J1-1 ), INCA, 291 $ AB( L, J1 ), INCA, D( J1 ), 292 $ WORK( J1 ), KD1 ) 293 10 CONTINUE 294* 295 ELSE 296 JEND = J1 + ( NR-1 )*KD1 297 DO 20 JINC = J1, JEND, KD1 298 CALL CROT( KDM1, AB( 2, JINC-1 ), 1, 299 $ AB( 1, JINC ), 1, D( JINC ), 300 $ WORK( JINC ) ) 301 20 CONTINUE 302 END IF 303 END IF 304* 305* 306 IF( K.GT.2 ) THEN 307 IF( K.LE.N-I+1 ) THEN 308* 309* generate plane rotation to annihilate a(i,i+k-1) 310* within the band 311* 312 CALL CLARTG( AB( KD-K+3, I+K-2 ), 313 $ AB( KD-K+2, I+K-1 ), D( I+K-1 ), 314 $ WORK( I+K-1 ), TEMP ) 315 AB( KD-K+3, I+K-2 ) = TEMP 316* 317* apply rotation from the right 318* 319 CALL CROT( K-3, AB( KD-K+4, I+K-2 ), 1, 320 $ AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ), 321 $ WORK( I+K-1 ) ) 322 END IF 323 NR = NR + 1 324 J1 = J1 - KDN - 1 325 END IF 326* 327* apply plane rotations from both sides to diagonal 328* blocks 329* 330 IF( NR.GT.0 ) 331 $ CALL CLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ), 332 $ AB( KD, J1 ), INCA, D( J1 ), 333 $ WORK( J1 ), KD1 ) 334* 335* apply plane rotations from the left 336* 337 IF( NR.GT.0 ) THEN 338 CALL CLACGV( NR, WORK( J1 ), KD1 ) 339 IF( 2*KD-1.LT.NR ) THEN 340* 341* Dependent on the the number of diagonals either 342* CLARTV or CROT is used 343* 344 DO 30 L = 1, KD - 1 345 IF( J2+L.GT.N ) THEN 346 NRT = NR - 1 347 ELSE 348 NRT = NR 349 END IF 350 IF( NRT.GT.0 ) 351 $ CALL CLARTV( NRT, AB( KD-L, J1+L ), INCA, 352 $ AB( KD-L+1, J1+L ), INCA, 353 $ D( J1 ), WORK( J1 ), KD1 ) 354 30 CONTINUE 355 ELSE 356 J1END = J1 + KD1*( NR-2 ) 357 IF( J1END.GE.J1 ) THEN 358 DO 40 JIN = J1, J1END, KD1 359 CALL CROT( KD-1, AB( KD-1, JIN+1 ), INCX, 360 $ AB( KD, JIN+1 ), INCX, 361 $ D( JIN ), WORK( JIN ) ) 362 40 CONTINUE 363 END IF 364 LEND = MIN( KDM1, N-J2 ) 365 LAST = J1END + KD1 366 IF( LEND.GT.0 ) 367 $ CALL CROT( LEND, AB( KD-1, LAST+1 ), INCX, 368 $ AB( KD, LAST+1 ), INCX, D( LAST ), 369 $ WORK( LAST ) ) 370 END IF 371 END IF 372* 373 IF( WANTQ ) THEN 374* 375* accumulate product of plane rotations in Q 376* 377 IF( INITQ ) THEN 378* 379* take advantage of the fact that Q was 380* initially the Identity matrix 381* 382 IQEND = MAX( IQEND, J2 ) 383 I2 = MAX( 0, K-3 ) 384 IQAEND = 1 + I*KD 385 IF( K.EQ.2 ) 386 $ IQAEND = IQAEND + KD 387 IQAEND = MIN( IQAEND, IQEND ) 388 DO 50 J = J1, J2, KD1 389 IBL = I - I2 / KDM1 390 I2 = I2 + 1 391 IQB = MAX( 1, J-IBL ) 392 NQ = 1 + IQAEND - IQB 393 IQAEND = MIN( IQAEND+KD, IQEND ) 394 CALL CROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ), 395 $ 1, D( J ), CONJG( WORK( J ) ) ) 396 50 CONTINUE 397 ELSE 398* 399 DO 60 J = J1, J2, KD1 400 CALL CROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1, 401 $ D( J ), CONJG( WORK( J ) ) ) 402 60 CONTINUE 403 END IF 404* 405 END IF 406* 407 IF( J2+KDN.GT.N ) THEN 408* 409* adjust J2 to keep within the bounds of the matrix 410* 411 NR = NR - 1 412 J2 = J2 - KDN - 1 413 END IF 414* 415 DO 70 J = J1, J2, KD1 416* 417* create nonzero element a(j-1,j+kd) outside the band 418* and store it in WORK 419* 420 WORK( J+KD ) = WORK( J )*AB( 1, J+KD ) 421 AB( 1, J+KD ) = D( J )*AB( 1, J+KD ) 422 70 CONTINUE 423 80 CONTINUE 424 90 CONTINUE 425 END IF 426* 427 IF( KD.GT.0 ) THEN 428* 429* make off-diagonal elements real and copy them to E 430* 431 DO 100 I = 1, N - 1 432 T = AB( KD, I+1 ) 433 ABST = ABS( T ) 434 AB( KD, I+1 ) = ABST 435 E( I ) = ABST 436 IF( ABST.NE.ZERO ) THEN 437 T = T / ABST 438 ELSE 439 T = CONE 440 END IF 441 IF( I.LT.N-1 ) 442 $ AB( KD, I+2 ) = AB( KD, I+2 )*T 443 IF( WANTQ ) THEN 444 CALL CSCAL( N, CONJG( T ), Q( 1, I+1 ), 1 ) 445 END IF 446 100 CONTINUE 447 ELSE 448* 449* set E to zero if original matrix was diagonal 450* 451 DO 110 I = 1, N - 1 452 E( I ) = ZERO 453 110 CONTINUE 454 END IF 455* 456* copy diagonal elements to D 457* 458 DO 120 I = 1, N 459 D( I ) = REAL( AB( KD1, I ) ) 460 120 CONTINUE 461* 462 ELSE 463* 464 IF( KD.GT.1 ) THEN 465* 466* Reduce to complex Hermitian tridiagonal form, working with 467* the lower triangle 468* 469 NR = 0 470 J1 = KDN + 2 471 J2 = 1 472* 473 AB( 1, 1 ) = REAL( AB( 1, 1 ) ) 474 DO 210 I = 1, N - 2 475* 476* Reduce i-th column of matrix to tridiagonal form 477* 478 DO 200 K = KDN + 1, 2, -1 479 J1 = J1 + KDN 480 J2 = J2 + KDN 481* 482 IF( NR.GT.0 ) THEN 483* 484* generate plane rotations to annihilate nonzero 485* elements which have been created outside the band 486* 487 CALL CLARGV( NR, AB( KD1, J1-KD1 ), INCA, 488 $ WORK( J1 ), KD1, D( J1 ), KD1 ) 489* 490* apply plane rotations from one side 491* 492* 493* Dependent on the the number of diagonals either 494* CLARTV or CROT is used 495* 496 IF( NR.GT.2*KD-1 ) THEN 497 DO 130 L = 1, KD - 1 498 CALL CLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA, 499 $ AB( KD1-L+1, J1-KD1+L ), INCA, 500 $ D( J1 ), WORK( J1 ), KD1 ) 501 130 CONTINUE 502 ELSE 503 JEND = J1 + KD1*( NR-1 ) 504 DO 140 JINC = J1, JEND, KD1 505 CALL CROT( KDM1, AB( KD, JINC-KD ), INCX, 506 $ AB( KD1, JINC-KD ), INCX, 507 $ D( JINC ), WORK( JINC ) ) 508 140 CONTINUE 509 END IF 510* 511 END IF 512* 513 IF( K.GT.2 ) THEN 514 IF( K.LE.N-I+1 ) THEN 515* 516* generate plane rotation to annihilate a(i+k-1,i) 517* within the band 518* 519 CALL CLARTG( AB( K-1, I ), AB( K, I ), 520 $ D( I+K-1 ), WORK( I+K-1 ), TEMP ) 521 AB( K-1, I ) = TEMP 522* 523* apply rotation from the left 524* 525 CALL CROT( K-3, AB( K-2, I+1 ), LDAB-1, 526 $ AB( K-1, I+1 ), LDAB-1, D( I+K-1 ), 527 $ WORK( I+K-1 ) ) 528 END IF 529 NR = NR + 1 530 J1 = J1 - KDN - 1 531 END IF 532* 533* apply plane rotations from both sides to diagonal 534* blocks 535* 536 IF( NR.GT.0 ) 537 $ CALL CLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ), 538 $ AB( 2, J1-1 ), INCA, D( J1 ), 539 $ WORK( J1 ), KD1 ) 540* 541* apply plane rotations from the right 542* 543* 544* Dependent on the the number of diagonals either 545* CLARTV or CROT is used 546* 547 IF( NR.GT.0 ) THEN 548 CALL CLACGV( NR, WORK( J1 ), KD1 ) 549 IF( NR.GT.2*KD-1 ) THEN 550 DO 150 L = 1, KD - 1 551 IF( J2+L.GT.N ) THEN 552 NRT = NR - 1 553 ELSE 554 NRT = NR 555 END IF 556 IF( NRT.GT.0 ) 557 $ CALL CLARTV( NRT, AB( L+2, J1-1 ), INCA, 558 $ AB( L+1, J1 ), INCA, D( J1 ), 559 $ WORK( J1 ), KD1 ) 560 150 CONTINUE 561 ELSE 562 J1END = J1 + KD1*( NR-2 ) 563 IF( J1END.GE.J1 ) THEN 564 DO 160 J1INC = J1, J1END, KD1 565 CALL CROT( KDM1, AB( 3, J1INC-1 ), 1, 566 $ AB( 2, J1INC ), 1, D( J1INC ), 567 $ WORK( J1INC ) ) 568 160 CONTINUE 569 END IF 570 LEND = MIN( KDM1, N-J2 ) 571 LAST = J1END + KD1 572 IF( LEND.GT.0 ) 573 $ CALL CROT( LEND, AB( 3, LAST-1 ), 1, 574 $ AB( 2, LAST ), 1, D( LAST ), 575 $ WORK( LAST ) ) 576 END IF 577 END IF 578* 579* 580* 581 IF( WANTQ ) THEN 582* 583* accumulate product of plane rotations in Q 584* 585 IF( INITQ ) THEN 586* 587* take advantage of the fact that Q was 588* initially the Identity matrix 589* 590 IQEND = MAX( IQEND, J2 ) 591 I2 = MAX( 0, K-3 ) 592 IQAEND = 1 + I*KD 593 IF( K.EQ.2 ) 594 $ IQAEND = IQAEND + KD 595 IQAEND = MIN( IQAEND, IQEND ) 596 DO 170 J = J1, J2, KD1 597 IBL = I - I2 / KDM1 598 I2 = I2 + 1 599 IQB = MAX( 1, J-IBL ) 600 NQ = 1 + IQAEND - IQB 601 IQAEND = MIN( IQAEND+KD, IQEND ) 602 CALL CROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ), 603 $ 1, D( J ), WORK( J ) ) 604 170 CONTINUE 605 ELSE 606* 607 DO 180 J = J1, J2, KD1 608 CALL CROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1, 609 $ D( J ), WORK( J ) ) 610 180 CONTINUE 611 END IF 612 END IF 613* 614 IF( J2+KDN.GT.N ) THEN 615* 616* adjust J2 to keep within the bounds of the matrix 617* 618 NR = NR - 1 619 J2 = J2 - KDN - 1 620 END IF 621* 622 DO 190 J = J1, J2, KD1 623* 624* create nonzero element a(j+kd,j-1) outside the 625* band and store it in WORK 626* 627 WORK( J+KD ) = WORK( J )*AB( KD1, J ) 628 AB( KD1, J ) = D( J )*AB( KD1, J ) 629 190 CONTINUE 630 200 CONTINUE 631 210 CONTINUE 632 END IF 633* 634 IF( KD.GT.0 ) THEN 635* 636* make off-diagonal elements real and copy them to E 637* 638 DO 220 I = 1, N - 1 639 T = AB( 2, I ) 640 ABST = ABS( T ) 641 AB( 2, I ) = ABST 642 E( I ) = ABST 643 IF( ABST.NE.ZERO ) THEN 644 T = T / ABST 645 ELSE 646 T = CONE 647 END IF 648 IF( I.LT.N-1 ) 649 $ AB( 2, I+1 ) = AB( 2, I+1 )*T 650 IF( WANTQ ) THEN 651 CALL CSCAL( N, T, Q( 1, I+1 ), 1 ) 652 END IF 653 220 CONTINUE 654 ELSE 655* 656* set E to zero if original matrix was diagonal 657* 658 DO 230 I = 1, N - 1 659 E( I ) = ZERO 660 230 CONTINUE 661 END IF 662* 663* copy diagonal elements to D 664* 665 DO 240 I = 1, N 666 D( I ) = REAL( AB( 1, I ) ) 667 240 CONTINUE 668 END IF 669* 670 RETURN 671* 672* End of CHBTRD 673* 674 END 675