Annotation of rpl/lapack/lapack/dsbtrd.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DSBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
! 2: $ WORK, INFO )
! 3: *
! 4: * -- LAPACK routine (version 3.2) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 8: *
! 9: * .. Scalar Arguments ..
! 10: CHARACTER UPLO, VECT
! 11: INTEGER INFO, KD, LDAB, LDQ, N
! 12: * ..
! 13: * .. Array Arguments ..
! 14: DOUBLE PRECISION AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ),
! 15: $ WORK( * )
! 16: * ..
! 17: *
! 18: * Purpose
! 19: * =======
! 20: *
! 21: * DSBTRD reduces a real symmetric band matrix A to symmetric
! 22: * tridiagonal form T by an orthogonal similarity transformation:
! 23: * Q**T * A * Q = T.
! 24: *
! 25: * Arguments
! 26: * =========
! 27: *
! 28: * VECT (input) CHARACTER*1
! 29: * = 'N': do not form Q;
! 30: * = 'V': form Q;
! 31: * = 'U': update a matrix X, by forming X*Q.
! 32: *
! 33: * UPLO (input) CHARACTER*1
! 34: * = 'U': Upper triangle of A is stored;
! 35: * = 'L': Lower triangle of A is stored.
! 36: *
! 37: * N (input) INTEGER
! 38: * The order of the matrix A. N >= 0.
! 39: *
! 40: * KD (input) INTEGER
! 41: * The number of superdiagonals of the matrix A if UPLO = 'U',
! 42: * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
! 43: *
! 44: * AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
! 45: * On entry, the upper or lower triangle of the symmetric band
! 46: * matrix A, stored in the first KD+1 rows of the array. The
! 47: * j-th column of A is stored in the j-th column of the array AB
! 48: * as follows:
! 49: * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
! 50: * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
! 51: * On exit, the diagonal elements of AB are overwritten by the
! 52: * diagonal elements of the tridiagonal matrix T; if KD > 0, the
! 53: * elements on the first superdiagonal (if UPLO = 'U') or the
! 54: * first subdiagonal (if UPLO = 'L') are overwritten by the
! 55: * off-diagonal elements of T; the rest of AB is overwritten by
! 56: * values generated during the reduction.
! 57: *
! 58: * LDAB (input) INTEGER
! 59: * The leading dimension of the array AB. LDAB >= KD+1.
! 60: *
! 61: * D (output) DOUBLE PRECISION array, dimension (N)
! 62: * The diagonal elements of the tridiagonal matrix T.
! 63: *
! 64: * E (output) DOUBLE PRECISION array, dimension (N-1)
! 65: * The off-diagonal elements of the tridiagonal matrix T:
! 66: * E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
! 67: *
! 68: * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
! 69: * On entry, if VECT = 'U', then Q must contain an N-by-N
! 70: * matrix X; if VECT = 'N' or 'V', then Q need not be set.
! 71: *
! 72: * On exit:
! 73: * if VECT = 'V', Q contains the N-by-N orthogonal matrix Q;
! 74: * if VECT = 'U', Q contains the product X*Q;
! 75: * if VECT = 'N', the array Q is not referenced.
! 76: *
! 77: * LDQ (input) INTEGER
! 78: * The leading dimension of the array Q.
! 79: * LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.
! 80: *
! 81: * WORK (workspace) DOUBLE PRECISION array, dimension (N)
! 82: *
! 83: * INFO (output) INTEGER
! 84: * = 0: successful exit
! 85: * < 0: if INFO = -i, the i-th argument had an illegal value
! 86: *
! 87: * Further Details
! 88: * ===============
! 89: *
! 90: * Modified by Linda Kaufman, Bell Labs.
! 91: *
! 92: * =====================================================================
! 93: *
! 94: * .. Parameters ..
! 95: DOUBLE PRECISION ZERO, ONE
! 96: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
! 97: * ..
! 98: * .. Local Scalars ..
! 99: LOGICAL INITQ, UPPER, WANTQ
! 100: INTEGER I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J,
! 101: $ J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1,
! 102: $ KDM1, KDN, L, LAST, LEND, NQ, NR, NRT
! 103: DOUBLE PRECISION TEMP
! 104: * ..
! 105: * .. External Subroutines ..
! 106: EXTERNAL DLAR2V, DLARGV, DLARTG, DLARTV, DLASET, DROT,
! 107: $ XERBLA
! 108: * ..
! 109: * .. Intrinsic Functions ..
! 110: INTRINSIC MAX, MIN
! 111: * ..
! 112: * .. External Functions ..
! 113: LOGICAL LSAME
! 114: EXTERNAL LSAME
! 115: * ..
! 116: * .. Executable Statements ..
! 117: *
! 118: * Test the input parameters
! 119: *
! 120: INITQ = LSAME( VECT, 'V' )
! 121: WANTQ = INITQ .OR. LSAME( VECT, 'U' )
! 122: UPPER = LSAME( UPLO, 'U' )
! 123: KD1 = KD + 1
! 124: KDM1 = KD - 1
! 125: INCX = LDAB - 1
! 126: IQEND = 1
! 127: *
! 128: INFO = 0
! 129: IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
! 130: INFO = -1
! 131: ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 132: INFO = -2
! 133: ELSE IF( N.LT.0 ) THEN
! 134: INFO = -3
! 135: ELSE IF( KD.LT.0 ) THEN
! 136: INFO = -4
! 137: ELSE IF( LDAB.LT.KD1 ) THEN
! 138: INFO = -6
! 139: ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN
! 140: INFO = -10
! 141: END IF
! 142: IF( INFO.NE.0 ) THEN
! 143: CALL XERBLA( 'DSBTRD', -INFO )
! 144: RETURN
! 145: END IF
! 146: *
! 147: * Quick return if possible
! 148: *
! 149: IF( N.EQ.0 )
! 150: $ RETURN
! 151: *
! 152: * Initialize Q to the unit matrix, if needed
! 153: *
! 154: IF( INITQ )
! 155: $ CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
! 156: *
! 157: * Wherever possible, plane rotations are generated and applied in
! 158: * vector operations of length NR over the index set J1:J2:KD1.
! 159: *
! 160: * The cosines and sines of the plane rotations are stored in the
! 161: * arrays D and WORK.
! 162: *
! 163: INCA = KD1*LDAB
! 164: KDN = MIN( N-1, KD )
! 165: IF( UPPER ) THEN
! 166: *
! 167: IF( KD.GT.1 ) THEN
! 168: *
! 169: * Reduce to tridiagonal form, working with upper triangle
! 170: *
! 171: NR = 0
! 172: J1 = KDN + 2
! 173: J2 = 1
! 174: *
! 175: DO 90 I = 1, N - 2
! 176: *
! 177: * Reduce i-th row of matrix to tridiagonal form
! 178: *
! 179: DO 80 K = KDN + 1, 2, -1
! 180: J1 = J1 + KDN
! 181: J2 = J2 + KDN
! 182: *
! 183: IF( NR.GT.0 ) THEN
! 184: *
! 185: * generate plane rotations to annihilate nonzero
! 186: * elements which have been created outside the band
! 187: *
! 188: CALL DLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ),
! 189: $ KD1, D( J1 ), KD1 )
! 190: *
! 191: * apply rotations from the right
! 192: *
! 193: *
! 194: * Dependent on the the number of diagonals either
! 195: * DLARTV or DROT is used
! 196: *
! 197: IF( NR.GE.2*KD-1 ) THEN
! 198: DO 10 L = 1, KD - 1
! 199: CALL DLARTV( NR, AB( L+1, J1-1 ), INCA,
! 200: $ AB( L, J1 ), INCA, D( J1 ),
! 201: $ WORK( J1 ), KD1 )
! 202: 10 CONTINUE
! 203: *
! 204: ELSE
! 205: JEND = J1 + ( NR-1 )*KD1
! 206: DO 20 JINC = J1, JEND, KD1
! 207: CALL DROT( KDM1, AB( 2, JINC-1 ), 1,
! 208: $ AB( 1, JINC ), 1, D( JINC ),
! 209: $ WORK( JINC ) )
! 210: 20 CONTINUE
! 211: END IF
! 212: END IF
! 213: *
! 214: *
! 215: IF( K.GT.2 ) THEN
! 216: IF( K.LE.N-I+1 ) THEN
! 217: *
! 218: * generate plane rotation to annihilate a(i,i+k-1)
! 219: * within the band
! 220: *
! 221: CALL DLARTG( AB( KD-K+3, I+K-2 ),
! 222: $ AB( KD-K+2, I+K-1 ), D( I+K-1 ),
! 223: $ WORK( I+K-1 ), TEMP )
! 224: AB( KD-K+3, I+K-2 ) = TEMP
! 225: *
! 226: * apply rotation from the right
! 227: *
! 228: CALL DROT( K-3, AB( KD-K+4, I+K-2 ), 1,
! 229: $ AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ),
! 230: $ WORK( I+K-1 ) )
! 231: END IF
! 232: NR = NR + 1
! 233: J1 = J1 - KDN - 1
! 234: END IF
! 235: *
! 236: * apply plane rotations from both sides to diagonal
! 237: * blocks
! 238: *
! 239: IF( NR.GT.0 )
! 240: $ CALL DLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ),
! 241: $ AB( KD, J1 ), INCA, D( J1 ),
! 242: $ WORK( J1 ), KD1 )
! 243: *
! 244: * apply plane rotations from the left
! 245: *
! 246: IF( NR.GT.0 ) THEN
! 247: IF( 2*KD-1.LT.NR ) THEN
! 248: *
! 249: * Dependent on the the number of diagonals either
! 250: * DLARTV or DROT is used
! 251: *
! 252: DO 30 L = 1, KD - 1
! 253: IF( J2+L.GT.N ) THEN
! 254: NRT = NR - 1
! 255: ELSE
! 256: NRT = NR
! 257: END IF
! 258: IF( NRT.GT.0 )
! 259: $ CALL DLARTV( NRT, AB( KD-L, J1+L ), INCA,
! 260: $ AB( KD-L+1, J1+L ), INCA,
! 261: $ D( J1 ), WORK( J1 ), KD1 )
! 262: 30 CONTINUE
! 263: ELSE
! 264: J1END = J1 + KD1*( NR-2 )
! 265: IF( J1END.GE.J1 ) THEN
! 266: DO 40 JIN = J1, J1END, KD1
! 267: CALL DROT( KD-1, AB( KD-1, JIN+1 ), INCX,
! 268: $ AB( KD, JIN+1 ), INCX,
! 269: $ D( JIN ), WORK( JIN ) )
! 270: 40 CONTINUE
! 271: END IF
! 272: LEND = MIN( KDM1, N-J2 )
! 273: LAST = J1END + KD1
! 274: IF( LEND.GT.0 )
! 275: $ CALL DROT( LEND, AB( KD-1, LAST+1 ), INCX,
! 276: $ AB( KD, LAST+1 ), INCX, D( LAST ),
! 277: $ WORK( LAST ) )
! 278: END IF
! 279: END IF
! 280: *
! 281: IF( WANTQ ) THEN
! 282: *
! 283: * accumulate product of plane rotations in Q
! 284: *
! 285: IF( INITQ ) THEN
! 286: *
! 287: * take advantage of the fact that Q was
! 288: * initially the Identity matrix
! 289: *
! 290: IQEND = MAX( IQEND, J2 )
! 291: I2 = MAX( 0, K-3 )
! 292: IQAEND = 1 + I*KD
! 293: IF( K.EQ.2 )
! 294: $ IQAEND = IQAEND + KD
! 295: IQAEND = MIN( IQAEND, IQEND )
! 296: DO 50 J = J1, J2, KD1
! 297: IBL = I - I2 / KDM1
! 298: I2 = I2 + 1
! 299: IQB = MAX( 1, J-IBL )
! 300: NQ = 1 + IQAEND - IQB
! 301: IQAEND = MIN( IQAEND+KD, IQEND )
! 302: CALL DROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
! 303: $ 1, D( J ), WORK( J ) )
! 304: 50 CONTINUE
! 305: ELSE
! 306: *
! 307: DO 60 J = J1, J2, KD1
! 308: CALL DROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
! 309: $ D( J ), WORK( J ) )
! 310: 60 CONTINUE
! 311: END IF
! 312: *
! 313: END IF
! 314: *
! 315: IF( J2+KDN.GT.N ) THEN
! 316: *
! 317: * adjust J2 to keep within the bounds of the matrix
! 318: *
! 319: NR = NR - 1
! 320: J2 = J2 - KDN - 1
! 321: END IF
! 322: *
! 323: DO 70 J = J1, J2, KD1
! 324: *
! 325: * create nonzero element a(j-1,j+kd) outside the band
! 326: * and store it in WORK
! 327: *
! 328: WORK( J+KD ) = WORK( J )*AB( 1, J+KD )
! 329: AB( 1, J+KD ) = D( J )*AB( 1, J+KD )
! 330: 70 CONTINUE
! 331: 80 CONTINUE
! 332: 90 CONTINUE
! 333: END IF
! 334: *
! 335: IF( KD.GT.0 ) THEN
! 336: *
! 337: * copy off-diagonal elements to E
! 338: *
! 339: DO 100 I = 1, N - 1
! 340: E( I ) = AB( KD, I+1 )
! 341: 100 CONTINUE
! 342: ELSE
! 343: *
! 344: * set E to zero if original matrix was diagonal
! 345: *
! 346: DO 110 I = 1, N - 1
! 347: E( I ) = ZERO
! 348: 110 CONTINUE
! 349: END IF
! 350: *
! 351: * copy diagonal elements to D
! 352: *
! 353: DO 120 I = 1, N
! 354: D( I ) = AB( KD1, I )
! 355: 120 CONTINUE
! 356: *
! 357: ELSE
! 358: *
! 359: IF( KD.GT.1 ) THEN
! 360: *
! 361: * Reduce to tridiagonal form, working with lower triangle
! 362: *
! 363: NR = 0
! 364: J1 = KDN + 2
! 365: J2 = 1
! 366: *
! 367: DO 210 I = 1, N - 2
! 368: *
! 369: * Reduce i-th column of matrix to tridiagonal form
! 370: *
! 371: DO 200 K = KDN + 1, 2, -1
! 372: J1 = J1 + KDN
! 373: J2 = J2 + KDN
! 374: *
! 375: IF( NR.GT.0 ) THEN
! 376: *
! 377: * generate plane rotations to annihilate nonzero
! 378: * elements which have been created outside the band
! 379: *
! 380: CALL DLARGV( NR, AB( KD1, J1-KD1 ), INCA,
! 381: $ WORK( J1 ), KD1, D( J1 ), KD1 )
! 382: *
! 383: * apply plane rotations from one side
! 384: *
! 385: *
! 386: * Dependent on the the number of diagonals either
! 387: * DLARTV or DROT is used
! 388: *
! 389: IF( NR.GT.2*KD-1 ) THEN
! 390: DO 130 L = 1, KD - 1
! 391: CALL DLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA,
! 392: $ AB( KD1-L+1, J1-KD1+L ), INCA,
! 393: $ D( J1 ), WORK( J1 ), KD1 )
! 394: 130 CONTINUE
! 395: ELSE
! 396: JEND = J1 + KD1*( NR-1 )
! 397: DO 140 JINC = J1, JEND, KD1
! 398: CALL DROT( KDM1, AB( KD, JINC-KD ), INCX,
! 399: $ AB( KD1, JINC-KD ), INCX,
! 400: $ D( JINC ), WORK( JINC ) )
! 401: 140 CONTINUE
! 402: END IF
! 403: *
! 404: END IF
! 405: *
! 406: IF( K.GT.2 ) THEN
! 407: IF( K.LE.N-I+1 ) THEN
! 408: *
! 409: * generate plane rotation to annihilate a(i+k-1,i)
! 410: * within the band
! 411: *
! 412: CALL DLARTG( AB( K-1, I ), AB( K, I ),
! 413: $ D( I+K-1 ), WORK( I+K-1 ), TEMP )
! 414: AB( K-1, I ) = TEMP
! 415: *
! 416: * apply rotation from the left
! 417: *
! 418: CALL DROT( K-3, AB( K-2, I+1 ), LDAB-1,
! 419: $ AB( K-1, I+1 ), LDAB-1, D( I+K-1 ),
! 420: $ WORK( I+K-1 ) )
! 421: END IF
! 422: NR = NR + 1
! 423: J1 = J1 - KDN - 1
! 424: END IF
! 425: *
! 426: * apply plane rotations from both sides to diagonal
! 427: * blocks
! 428: *
! 429: IF( NR.GT.0 )
! 430: $ CALL DLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ),
! 431: $ AB( 2, J1-1 ), INCA, D( J1 ),
! 432: $ WORK( J1 ), KD1 )
! 433: *
! 434: * apply plane rotations from the right
! 435: *
! 436: *
! 437: * Dependent on the the number of diagonals either
! 438: * DLARTV or DROT is used
! 439: *
! 440: IF( NR.GT.0 ) THEN
! 441: IF( NR.GT.2*KD-1 ) THEN
! 442: DO 150 L = 1, KD - 1
! 443: IF( J2+L.GT.N ) THEN
! 444: NRT = NR - 1
! 445: ELSE
! 446: NRT = NR
! 447: END IF
! 448: IF( NRT.GT.0 )
! 449: $ CALL DLARTV( NRT, AB( L+2, J1-1 ), INCA,
! 450: $ AB( L+1, J1 ), INCA, D( J1 ),
! 451: $ WORK( J1 ), KD1 )
! 452: 150 CONTINUE
! 453: ELSE
! 454: J1END = J1 + KD1*( NR-2 )
! 455: IF( J1END.GE.J1 ) THEN
! 456: DO 160 J1INC = J1, J1END, KD1
! 457: CALL DROT( KDM1, AB( 3, J1INC-1 ), 1,
! 458: $ AB( 2, J1INC ), 1, D( J1INC ),
! 459: $ WORK( J1INC ) )
! 460: 160 CONTINUE
! 461: END IF
! 462: LEND = MIN( KDM1, N-J2 )
! 463: LAST = J1END + KD1
! 464: IF( LEND.GT.0 )
! 465: $ CALL DROT( LEND, AB( 3, LAST-1 ), 1,
! 466: $ AB( 2, LAST ), 1, D( LAST ),
! 467: $ WORK( LAST ) )
! 468: END IF
! 469: END IF
! 470: *
! 471: *
! 472: *
! 473: IF( WANTQ ) THEN
! 474: *
! 475: * accumulate product of plane rotations in Q
! 476: *
! 477: IF( INITQ ) THEN
! 478: *
! 479: * take advantage of the fact that Q was
! 480: * initially the Identity matrix
! 481: *
! 482: IQEND = MAX( IQEND, J2 )
! 483: I2 = MAX( 0, K-3 )
! 484: IQAEND = 1 + I*KD
! 485: IF( K.EQ.2 )
! 486: $ IQAEND = IQAEND + KD
! 487: IQAEND = MIN( IQAEND, IQEND )
! 488: DO 170 J = J1, J2, KD1
! 489: IBL = I - I2 / KDM1
! 490: I2 = I2 + 1
! 491: IQB = MAX( 1, J-IBL )
! 492: NQ = 1 + IQAEND - IQB
! 493: IQAEND = MIN( IQAEND+KD, IQEND )
! 494: CALL DROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
! 495: $ 1, D( J ), WORK( J ) )
! 496: 170 CONTINUE
! 497: ELSE
! 498: *
! 499: DO 180 J = J1, J2, KD1
! 500: CALL DROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
! 501: $ D( J ), WORK( J ) )
! 502: 180 CONTINUE
! 503: END IF
! 504: END IF
! 505: *
! 506: IF( J2+KDN.GT.N ) THEN
! 507: *
! 508: * adjust J2 to keep within the bounds of the matrix
! 509: *
! 510: NR = NR - 1
! 511: J2 = J2 - KDN - 1
! 512: END IF
! 513: *
! 514: DO 190 J = J1, J2, KD1
! 515: *
! 516: * create nonzero element a(j+kd,j-1) outside the
! 517: * band and store it in WORK
! 518: *
! 519: WORK( J+KD ) = WORK( J )*AB( KD1, J )
! 520: AB( KD1, J ) = D( J )*AB( KD1, J )
! 521: 190 CONTINUE
! 522: 200 CONTINUE
! 523: 210 CONTINUE
! 524: END IF
! 525: *
! 526: IF( KD.GT.0 ) THEN
! 527: *
! 528: * copy off-diagonal elements to E
! 529: *
! 530: DO 220 I = 1, N - 1
! 531: E( I ) = AB( 2, I )
! 532: 220 CONTINUE
! 533: ELSE
! 534: *
! 535: * set E to zero if original matrix was diagonal
! 536: *
! 537: DO 230 I = 1, N - 1
! 538: E( I ) = ZERO
! 539: 230 CONTINUE
! 540: END IF
! 541: *
! 542: * copy diagonal elements to D
! 543: *
! 544: DO 240 I = 1, N
! 545: D( I ) = AB( 1, I )
! 546: 240 CONTINUE
! 547: END IF
! 548: *
! 549: RETURN
! 550: *
! 551: * End of DSBTRD
! 552: *
! 553: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsbtrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack