SRC/dsbtrd.f

   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