1*> \brief \b DSBTRD 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DSBTRD + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsbtrd.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsbtrd.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsbtrd.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DSBTRD( 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* DOUBLE PRECISION AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ), 30* $ WORK( * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> DSBTRD reduces a real symmetric band matrix A to symmetric 40*> tridiagonal form T by an orthogonal similarity transformation: 41*> Q**T * 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 DOUBLE PRECISION array, dimension (LDAB,N) 78*> On entry, the upper or lower triangle of the symmetric 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 DOUBLE PRECISION array, dimension (N) 101*> The diagonal elements of the tridiagonal matrix T. 102*> \endverbatim 103*> 104*> \param[out] E 105*> \verbatim 106*> E is DOUBLE PRECISION 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 DOUBLE PRECISION 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 orthogonal 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 DOUBLE PRECISION 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*> \date November 2011 151* 152*> \ingroup doubleOTHERcomputational 153* 154*> \par Further Details: 155* ===================== 156*> 157*> \verbatim 158*> 159*> Modified by Linda Kaufman, Bell Labs. 160*> \endverbatim 161*> 162* ===================================================================== 163 SUBROUTINE DSBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, 164 $ WORK, INFO ) 165* 166* -- LAPACK computational routine (version 3.4.0) -- 167* -- LAPACK is a software package provided by Univ. of Tennessee, -- 168* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 169* November 2011 170* 171* .. Scalar Arguments .. 172 CHARACTER UPLO, VECT 173 INTEGER INFO, KD, LDAB, LDQ, N 174* .. 175* .. Array Arguments .. 176 DOUBLE PRECISION AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ), 177 $ WORK( * ) 178* .. 179* 180* ===================================================================== 181* 182* .. Parameters .. 183 DOUBLE PRECISION ZERO, ONE 184 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+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 DOUBLE PRECISION TEMP 192* .. 193* .. External Subroutines .. 194 EXTERNAL DLAR2V, DLARGV, DLARTG, DLARTV, DLASET, DROT, 195 $ XERBLA 196* .. 197* .. Intrinsic Functions .. 198 INTRINSIC MAX, MIN 199* .. 200* .. External Functions .. 201 LOGICAL LSAME 202 EXTERNAL LSAME 203* .. 204* .. Executable Statements .. 205* 206* Test the input parameters 207* 208 INITQ = LSAME( VECT, 'V' ) 209 WANTQ = INITQ .OR. LSAME( VECT, 'U' ) 210 UPPER = LSAME( UPLO, 'U' ) 211 KD1 = KD + 1 212 KDM1 = KD - 1 213 INCX = LDAB - 1 214 IQEND = 1 215* 216 INFO = 0 217 IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'N' ) ) THEN 218 INFO = -1 219 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 220 INFO = -2 221 ELSE IF( N.LT.0 ) THEN 222 INFO = -3 223 ELSE IF( KD.LT.0 ) THEN 224 INFO = -4 225 ELSE IF( LDAB.LT.KD1 ) THEN 226 INFO = -6 227 ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN 228 INFO = -10 229 END IF 230 IF( INFO.NE.0 ) THEN 231 CALL XERBLA( 'DSBTRD', -INFO ) 232 RETURN 233 END IF 234* 235* Quick return if possible 236* 237 IF( N.EQ.0 ) 238 $ RETURN 239* 240* Initialize Q to the unit matrix, if needed 241* 242 IF( INITQ ) 243 $ CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ ) 244* 245* Wherever possible, plane rotations are generated and applied in 246* vector operations of length NR over the index set J1:J2:KD1. 247* 248* The cosines and sines of the plane rotations are stored in the 249* arrays D and WORK. 250* 251 INCA = KD1*LDAB 252 KDN = MIN( N-1, KD ) 253 IF( UPPER ) THEN 254* 255 IF( KD.GT.1 ) THEN 256* 257* Reduce to tridiagonal form, working with upper triangle 258* 259 NR = 0 260 J1 = KDN + 2 261 J2 = 1 262* 263 DO 90 I = 1, N - 2 264* 265* Reduce i-th row of matrix to tridiagonal form 266* 267 DO 80 K = KDN + 1, 2, -1 268 J1 = J1 + KDN 269 J2 = J2 + KDN 270* 271 IF( NR.GT.0 ) THEN 272* 273* generate plane rotations to annihilate nonzero 274* elements which have been created outside the band 275* 276 CALL DLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ), 277 $ KD1, D( J1 ), KD1 ) 278* 279* apply rotations from the right 280* 281* 282* Dependent on the the number of diagonals either 283* DLARTV or DROT is used 284* 285 IF( NR.GE.2*KD-1 ) THEN 286 DO 10 L = 1, KD - 1 287 CALL DLARTV( NR, AB( L+1, J1-1 ), INCA, 288 $ AB( L, J1 ), INCA, D( J1 ), 289 $ WORK( J1 ), KD1 ) 290 10 CONTINUE 291* 292 ELSE 293 JEND = J1 + ( NR-1 )*KD1 294 DO 20 JINC = J1, JEND, KD1 295 CALL DROT( KDM1, AB( 2, JINC-1 ), 1, 296 $ AB( 1, JINC ), 1, D( JINC ), 297 $ WORK( JINC ) ) 298 20 CONTINUE 299 END IF 300 END IF 301* 302* 303 IF( K.GT.2 ) THEN 304 IF( K.LE.N-I+1 ) THEN 305* 306* generate plane rotation to annihilate a(i,i+k-1) 307* within the band 308* 309 CALL DLARTG( AB( KD-K+3, I+K-2 ), 310 $ AB( KD-K+2, I+K-1 ), D( I+K-1 ), 311 $ WORK( I+K-1 ), TEMP ) 312 AB( KD-K+3, I+K-2 ) = TEMP 313* 314* apply rotation from the right 315* 316 CALL DROT( K-3, AB( KD-K+4, I+K-2 ), 1, 317 $ AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ), 318 $ WORK( I+K-1 ) ) 319 END IF 320 NR = NR + 1 321 J1 = J1 - KDN - 1 322 END IF 323* 324* apply plane rotations from both sides to diagonal 325* blocks 326* 327 IF( NR.GT.0 ) 328 $ CALL DLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ), 329 $ AB( KD, J1 ), INCA, D( J1 ), 330 $ WORK( J1 ), KD1 ) 331* 332* apply plane rotations from the left 333* 334 IF( NR.GT.0 ) THEN 335 IF( 2*KD-1.LT.NR ) THEN 336* 337* Dependent on the the number of diagonals either 338* DLARTV or DROT is used 339* 340 DO 30 L = 1, KD - 1 341 IF( J2+L.GT.N ) THEN 342 NRT = NR - 1 343 ELSE 344 NRT = NR 345 END IF 346 IF( NRT.GT.0 ) 347 $ CALL DLARTV( NRT, AB( KD-L, J1+L ), INCA, 348 $ AB( KD-L+1, J1+L ), INCA, 349 $ D( J1 ), WORK( J1 ), KD1 ) 350 30 CONTINUE 351 ELSE 352 J1END = J1 + KD1*( NR-2 ) 353 IF( J1END.GE.J1 ) THEN 354 DO 40 JIN = J1, J1END, KD1 355 CALL DROT( KD-1, AB( KD-1, JIN+1 ), INCX, 356 $ AB( KD, JIN+1 ), INCX, 357 $ D( JIN ), WORK( JIN ) ) 358 40 CONTINUE 359 END IF 360 LEND = MIN( KDM1, N-J2 ) 361 LAST = J1END + KD1 362 IF( LEND.GT.0 ) 363 $ CALL DROT( LEND, AB( KD-1, LAST+1 ), INCX, 364 $ AB( KD, LAST+1 ), INCX, D( LAST ), 365 $ WORK( LAST ) ) 366 END IF 367 END IF 368* 369 IF( WANTQ ) THEN 370* 371* accumulate product of plane rotations in Q 372* 373 IF( INITQ ) THEN 374* 375* take advantage of the fact that Q was 376* initially the Identity matrix 377* 378 IQEND = MAX( IQEND, J2 ) 379 I2 = MAX( 0, K-3 ) 380 IQAEND = 1 + I*KD 381 IF( K.EQ.2 ) 382 $ IQAEND = IQAEND + KD 383 IQAEND = MIN( IQAEND, IQEND ) 384 DO 50 J = J1, J2, KD1 385 IBL = I - I2 / KDM1 386 I2 = I2 + 1 387 IQB = MAX( 1, J-IBL ) 388 NQ = 1 + IQAEND - IQB 389 IQAEND = MIN( IQAEND+KD, IQEND ) 390 CALL DROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ), 391 $ 1, D( J ), WORK( J ) ) 392 50 CONTINUE 393 ELSE 394* 395 DO 60 J = J1, J2, KD1 396 CALL DROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1, 397 $ D( J ), WORK( J ) ) 398 60 CONTINUE 399 END IF 400* 401 END IF 402* 403 IF( J2+KDN.GT.N ) THEN 404* 405* adjust J2 to keep within the bounds of the matrix 406* 407 NR = NR - 1 408 J2 = J2 - KDN - 1 409 END IF 410* 411 DO 70 J = J1, J2, KD1 412* 413* create nonzero element a(j-1,j+kd) outside the band 414* and store it in WORK 415* 416 WORK( J+KD ) = WORK( J )*AB( 1, J+KD ) 417 AB( 1, J+KD ) = D( J )*AB( 1, J+KD ) 418 70 CONTINUE 419 80 CONTINUE 420 90 CONTINUE 421 END IF 422* 423 IF( KD.GT.0 ) THEN 424* 425* copy off-diagonal elements to E 426* 427 DO 100 I = 1, N - 1 428 E( I ) = AB( KD, I+1 ) 429 100 CONTINUE 430 ELSE 431* 432* set E to zero if original matrix was diagonal 433* 434 DO 110 I = 1, N - 1 435 E( I ) = ZERO 436 110 CONTINUE 437 END IF 438* 439* copy diagonal elements to D 440* 441 DO 120 I = 1, N 442 D( I ) = AB( KD1, I ) 443 120 CONTINUE 444* 445 ELSE 446* 447 IF( KD.GT.1 ) THEN 448* 449* Reduce to tridiagonal form, working with lower triangle 450* 451 NR = 0 452 J1 = KDN + 2 453 J2 = 1 454* 455 DO 210 I = 1, N - 2 456* 457* Reduce i-th column of matrix to tridiagonal form 458* 459 DO 200 K = KDN + 1, 2, -1 460 J1 = J1 + KDN 461 J2 = J2 + KDN 462* 463 IF( NR.GT.0 ) THEN 464* 465* generate plane rotations to annihilate nonzero 466* elements which have been created outside the band 467* 468 CALL DLARGV( NR, AB( KD1, J1-KD1 ), INCA, 469 $ WORK( J1 ), KD1, D( J1 ), KD1 ) 470* 471* apply plane rotations from one side 472* 473* 474* Dependent on the the number of diagonals either 475* DLARTV or DROT is used 476* 477 IF( NR.GT.2*KD-1 ) THEN 478 DO 130 L = 1, KD - 1 479 CALL DLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA, 480 $ AB( KD1-L+1, J1-KD1+L ), INCA, 481 $ D( J1 ), WORK( J1 ), KD1 ) 482 130 CONTINUE 483 ELSE 484 JEND = J1 + KD1*( NR-1 ) 485 DO 140 JINC = J1, JEND, KD1 486 CALL DROT( KDM1, AB( KD, JINC-KD ), INCX, 487 $ AB( KD1, JINC-KD ), INCX, 488 $ D( JINC ), WORK( JINC ) ) 489 140 CONTINUE 490 END IF 491* 492 END IF 493* 494 IF( K.GT.2 ) THEN 495 IF( K.LE.N-I+1 ) THEN 496* 497* generate plane rotation to annihilate a(i+k-1,i) 498* within the band 499* 500 CALL DLARTG( AB( K-1, I ), AB( K, I ), 501 $ D( I+K-1 ), WORK( I+K-1 ), TEMP ) 502 AB( K-1, I ) = TEMP 503* 504* apply rotation from the left 505* 506 CALL DROT( K-3, AB( K-2, I+1 ), LDAB-1, 507 $ AB( K-1, I+1 ), LDAB-1, D( I+K-1 ), 508 $ WORK( I+K-1 ) ) 509 END IF 510 NR = NR + 1 511 J1 = J1 - KDN - 1 512 END IF 513* 514* apply plane rotations from both sides to diagonal 515* blocks 516* 517 IF( NR.GT.0 ) 518 $ CALL DLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ), 519 $ AB( 2, J1-1 ), INCA, D( J1 ), 520 $ WORK( J1 ), KD1 ) 521* 522* apply plane rotations from the right 523* 524* 525* Dependent on the the number of diagonals either 526* DLARTV or DROT is used 527* 528 IF( NR.GT.0 ) THEN 529 IF( NR.GT.2*KD-1 ) THEN 530 DO 150 L = 1, KD - 1 531 IF( J2+L.GT.N ) THEN 532 NRT = NR - 1 533 ELSE 534 NRT = NR 535 END IF 536 IF( NRT.GT.0 ) 537 $ CALL DLARTV( NRT, AB( L+2, J1-1 ), INCA, 538 $ AB( L+1, J1 ), INCA, D( J1 ), 539 $ WORK( J1 ), KD1 ) 540 150 CONTINUE 541 ELSE 542 J1END = J1 + KD1*( NR-2 ) 543 IF( J1END.GE.J1 ) THEN 544 DO 160 J1INC = J1, J1END, KD1 545 CALL DROT( KDM1, AB( 3, J1INC-1 ), 1, 546 $ AB( 2, J1INC ), 1, D( J1INC ), 547 $ WORK( J1INC ) ) 548 160 CONTINUE 549 END IF 550 LEND = MIN( KDM1, N-J2 ) 551 LAST = J1END + KD1 552 IF( LEND.GT.0 ) 553 $ CALL DROT( LEND, AB( 3, LAST-1 ), 1, 554 $ AB( 2, LAST ), 1, D( LAST ), 555 $ WORK( LAST ) ) 556 END IF 557 END IF 558* 559* 560* 561 IF( WANTQ ) THEN 562* 563* accumulate product of plane rotations in Q 564* 565 IF( INITQ ) THEN 566* 567* take advantage of the fact that Q was 568* initially the Identity matrix 569* 570 IQEND = MAX( IQEND, J2 ) 571 I2 = MAX( 0, K-3 ) 572 IQAEND = 1 + I*KD 573 IF( K.EQ.2 ) 574 $ IQAEND = IQAEND + KD 575 IQAEND = MIN( IQAEND, IQEND ) 576 DO 170 J = J1, J2, KD1 577 IBL = I - I2 / KDM1 578 I2 = I2 + 1 579 IQB = MAX( 1, J-IBL ) 580 NQ = 1 + IQAEND - IQB 581 IQAEND = MIN( IQAEND+KD, IQEND ) 582 CALL DROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ), 583 $ 1, D( J ), WORK( J ) ) 584 170 CONTINUE 585 ELSE 586* 587 DO 180 J = J1, J2, KD1 588 CALL DROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1, 589 $ D( J ), WORK( J ) ) 590 180 CONTINUE 591 END IF 592 END IF 593* 594 IF( J2+KDN.GT.N ) THEN 595* 596* adjust J2 to keep within the bounds of the matrix 597* 598 NR = NR - 1 599 J2 = J2 - KDN - 1 600 END IF 601* 602 DO 190 J = J1, J2, KD1 603* 604* create nonzero element a(j+kd,j-1) outside the 605* band and store it in WORK 606* 607 WORK( J+KD ) = WORK( J )*AB( KD1, J ) 608 AB( KD1, J ) = D( J )*AB( KD1, J ) 609 190 CONTINUE 610 200 CONTINUE 611 210 CONTINUE 612 END IF 613* 614 IF( KD.GT.0 ) THEN 615* 616* copy off-diagonal elements to E 617* 618 DO 220 I = 1, N - 1 619 E( I ) = AB( 2, I ) 620 220 CONTINUE 621 ELSE 622* 623* set E to zero if original matrix was diagonal 624* 625 DO 230 I = 1, N - 1 626 E( I ) = ZERO 627 230 CONTINUE 628 END IF 629* 630* copy diagonal elements to D 631* 632 DO 240 I = 1, N 633 D( I ) = AB( 1, I ) 634 240 CONTINUE 635 END IF 636* 637 RETURN 638* 639* End of DSBTRD 640* 641 END 642