LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dsbtrd.f
Go to the documentation of this file.
00001 *> \brief \b DSBTRD
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DSBTRD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsbtrd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsbtrd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsbtrd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DSBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
00022 *                          WORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          UPLO, VECT
00026 *       INTEGER            INFO, KD, LDAB, LDQ, N
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       DOUBLE PRECISION   AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ),
00030 *      $                   WORK( * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> DSBTRD reduces a real symmetric band matrix A to symmetric
00040 *> tridiagonal form T by an orthogonal similarity transformation:
00041 *> Q**T * A * Q = T.
00042 *> \endverbatim
00043 *
00044 *  Arguments:
00045 *  ==========
00046 *
00047 *> \param[in] VECT
00048 *> \verbatim
00049 *>          VECT is CHARACTER*1
00050 *>          = 'N':  do not form Q;
00051 *>          = 'V':  form Q;
00052 *>          = 'U':  update a matrix X, by forming X*Q.
00053 *> \endverbatim
00054 *>
00055 *> \param[in] UPLO
00056 *> \verbatim
00057 *>          UPLO is CHARACTER*1
00058 *>          = 'U':  Upper triangle of A is stored;
00059 *>          = 'L':  Lower triangle of A is stored.
00060 *> \endverbatim
00061 *>
00062 *> \param[in] N
00063 *> \verbatim
00064 *>          N is INTEGER
00065 *>          The order of the matrix A.  N >= 0.
00066 *> \endverbatim
00067 *>
00068 *> \param[in] KD
00069 *> \verbatim
00070 *>          KD is INTEGER
00071 *>          The number of superdiagonals of the matrix A if UPLO = 'U',
00072 *>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
00073 *> \endverbatim
00074 *>
00075 *> \param[in,out] AB
00076 *> \verbatim
00077 *>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
00078 *>          On entry, the upper or lower triangle of the symmetric band
00079 *>          matrix A, stored in the first KD+1 rows of the array.  The
00080 *>          j-th column of A is stored in the j-th column of the array AB
00081 *>          as follows:
00082 *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
00083 *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
00084 *>          On exit, the diagonal elements of AB are overwritten by the
00085 *>          diagonal elements of the tridiagonal matrix T; if KD > 0, the
00086 *>          elements on the first superdiagonal (if UPLO = 'U') or the
00087 *>          first subdiagonal (if UPLO = 'L') are overwritten by the
00088 *>          off-diagonal elements of T; the rest of AB is overwritten by
00089 *>          values generated during the reduction.
00090 *> \endverbatim
00091 *>
00092 *> \param[in] LDAB
00093 *> \verbatim
00094 *>          LDAB is INTEGER
00095 *>          The leading dimension of the array AB.  LDAB >= KD+1.
00096 *> \endverbatim
00097 *>
00098 *> \param[out] D
00099 *> \verbatim
00100 *>          D is DOUBLE PRECISION array, dimension (N)
00101 *>          The diagonal elements of the tridiagonal matrix T.
00102 *> \endverbatim
00103 *>
00104 *> \param[out] E
00105 *> \verbatim
00106 *>          E is DOUBLE PRECISION array, dimension (N-1)
00107 *>          The off-diagonal elements of the tridiagonal matrix T:
00108 *>          E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
00109 *> \endverbatim
00110 *>
00111 *> \param[in,out] Q
00112 *> \verbatim
00113 *>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
00114 *>          On entry, if VECT = 'U', then Q must contain an N-by-N
00115 *>          matrix X; if VECT = 'N' or 'V', then Q need not be set.
00116 *>
00117 *>          On exit:
00118 *>          if VECT = 'V', Q contains the N-by-N orthogonal matrix Q;
00119 *>          if VECT = 'U', Q contains the product X*Q;
00120 *>          if VECT = 'N', the array Q is not referenced.
00121 *> \endverbatim
00122 *>
00123 *> \param[in] LDQ
00124 *> \verbatim
00125 *>          LDQ is INTEGER
00126 *>          The leading dimension of the array Q.
00127 *>          LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.
00128 *> \endverbatim
00129 *>
00130 *> \param[out] WORK
00131 *> \verbatim
00132 *>          WORK is DOUBLE PRECISION array, dimension (N)
00133 *> \endverbatim
00134 *>
00135 *> \param[out] INFO
00136 *> \verbatim
00137 *>          INFO is INTEGER
00138 *>          = 0:  successful exit
00139 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00140 *> \endverbatim
00141 *
00142 *  Authors:
00143 *  ========
00144 *
00145 *> \author Univ. of Tennessee 
00146 *> \author Univ. of California Berkeley 
00147 *> \author Univ. of Colorado Denver 
00148 *> \author NAG Ltd. 
00149 *
00150 *> \date November 2011
00151 *
00152 *> \ingroup doubleOTHERcomputational
00153 *
00154 *> \par Further Details:
00155 *  =====================
00156 *>
00157 *> \verbatim
00158 *>
00159 *>  Modified by Linda Kaufman, Bell Labs.
00160 *> \endverbatim
00161 *>
00162 *  =====================================================================
00163       SUBROUTINE DSBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
00164      $                   WORK, INFO )
00165 *
00166 *  -- LAPACK computational routine (version 3.4.0) --
00167 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00168 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00169 *     November 2011
00170 *
00171 *     .. Scalar Arguments ..
00172       CHARACTER          UPLO, VECT
00173       INTEGER            INFO, KD, LDAB, LDQ, N
00174 *     ..
00175 *     .. Array Arguments ..
00176       DOUBLE PRECISION   AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ),
00177      $                   WORK( * )
00178 *     ..
00179 *
00180 *  =====================================================================
00181 *
00182 *     .. Parameters ..
00183       DOUBLE PRECISION   ZERO, ONE
00184       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00185 *     ..
00186 *     .. Local Scalars ..
00187       LOGICAL            INITQ, UPPER, WANTQ
00188       INTEGER            I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J,
00189      $                   J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1,
00190      $                   KDM1, KDN, L, LAST, LEND, NQ, NR, NRT
00191       DOUBLE PRECISION   TEMP
00192 *     ..
00193 *     .. External Subroutines ..
00194       EXTERNAL           DLAR2V, DLARGV, DLARTG, DLARTV, DLASET, DROT,
00195      $                   XERBLA
00196 *     ..
00197 *     .. Intrinsic Functions ..
00198       INTRINSIC          MAX, MIN
00199 *     ..
00200 *     .. External Functions ..
00201       LOGICAL            LSAME
00202       EXTERNAL           LSAME
00203 *     ..
00204 *     .. Executable Statements ..
00205 *
00206 *     Test the input parameters
00207 *
00208       INITQ = LSAME( VECT, 'V' )
00209       WANTQ = INITQ .OR. LSAME( VECT, 'U' )
00210       UPPER = LSAME( UPLO, 'U' )
00211       KD1 = KD + 1
00212       KDM1 = KD - 1
00213       INCX = LDAB - 1
00214       IQEND = 1
00215 *
00216       INFO = 0
00217       IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
00218          INFO = -1
00219       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00220          INFO = -2
00221       ELSE IF( N.LT.0 ) THEN
00222          INFO = -3
00223       ELSE IF( KD.LT.0 ) THEN
00224          INFO = -4
00225       ELSE IF( LDAB.LT.KD1 ) THEN
00226          INFO = -6
00227       ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN
00228          INFO = -10
00229       END IF
00230       IF( INFO.NE.0 ) THEN
00231          CALL XERBLA( 'DSBTRD', -INFO )
00232          RETURN
00233       END IF
00234 *
00235 *     Quick return if possible
00236 *
00237       IF( N.EQ.0 )
00238      $   RETURN
00239 *
00240 *     Initialize Q to the unit matrix, if needed
00241 *
00242       IF( INITQ )
00243      $   CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
00244 *
00245 *     Wherever possible, plane rotations are generated and applied in
00246 *     vector operations of length NR over the index set J1:J2:KD1.
00247 *
00248 *     The cosines and sines of the plane rotations are stored in the
00249 *     arrays D and WORK.
00250 *
00251       INCA = KD1*LDAB
00252       KDN = MIN( N-1, KD )
00253       IF( UPPER ) THEN
00254 *
00255          IF( KD.GT.1 ) THEN
00256 *
00257 *           Reduce to tridiagonal form, working with upper triangle
00258 *
00259             NR = 0
00260             J1 = KDN + 2
00261             J2 = 1
00262 *
00263             DO 90 I = 1, N - 2
00264 *
00265 *              Reduce i-th row of matrix to tridiagonal form
00266 *
00267                DO 80 K = KDN + 1, 2, -1
00268                   J1 = J1 + KDN
00269                   J2 = J2 + KDN
00270 *
00271                   IF( NR.GT.0 ) THEN
00272 *
00273 *                    generate plane rotations to annihilate nonzero
00274 *                    elements which have been created outside the band
00275 *
00276                      CALL DLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ),
00277      $                            KD1, D( J1 ), KD1 )
00278 *
00279 *                    apply rotations from the right
00280 *
00281 *
00282 *                    Dependent on the the number of diagonals either
00283 *                    DLARTV or DROT is used
00284 *
00285                      IF( NR.GE.2*KD-1 ) THEN
00286                         DO 10 L = 1, KD - 1
00287                            CALL DLARTV( NR, AB( L+1, J1-1 ), INCA,
00288      $                                  AB( L, J1 ), INCA, D( J1 ),
00289      $                                  WORK( J1 ), KD1 )
00290    10                   CONTINUE
00291 *
00292                      ELSE
00293                         JEND = J1 + ( NR-1 )*KD1
00294                         DO 20 JINC = J1, JEND, KD1
00295                            CALL DROT( KDM1, AB( 2, JINC-1 ), 1,
00296      $                                AB( 1, JINC ), 1, D( JINC ),
00297      $                                WORK( JINC ) )
00298    20                   CONTINUE
00299                      END IF
00300                   END IF
00301 *
00302 *
00303                   IF( K.GT.2 ) THEN
00304                      IF( K.LE.N-I+1 ) THEN
00305 *
00306 *                       generate plane rotation to annihilate a(i,i+k-1)
00307 *                       within the band
00308 *
00309                         CALL DLARTG( AB( KD-K+3, I+K-2 ),
00310      $                               AB( KD-K+2, I+K-1 ), D( I+K-1 ),
00311      $                               WORK( I+K-1 ), TEMP )
00312                         AB( KD-K+3, I+K-2 ) = TEMP
00313 *
00314 *                       apply rotation from the right
00315 *
00316                         CALL DROT( K-3, AB( KD-K+4, I+K-2 ), 1,
00317      $                             AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ),
00318      $                             WORK( I+K-1 ) )
00319                      END IF
00320                      NR = NR + 1
00321                      J1 = J1 - KDN - 1
00322                   END IF
00323 *
00324 *                 apply plane rotations from both sides to diagonal
00325 *                 blocks
00326 *
00327                   IF( NR.GT.0 )
00328      $               CALL DLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ),
00329      $                            AB( KD, J1 ), INCA, D( J1 ),
00330      $                            WORK( J1 ), KD1 )
00331 *
00332 *                 apply plane rotations from the left
00333 *
00334                   IF( NR.GT.0 ) THEN
00335                      IF( 2*KD-1.LT.NR ) THEN
00336 *
00337 *                    Dependent on the the number of diagonals either
00338 *                    DLARTV or DROT is used
00339 *
00340                         DO 30 L = 1, KD - 1
00341                            IF( J2+L.GT.N ) THEN
00342                               NRT = NR - 1
00343                            ELSE
00344                               NRT = NR
00345                            END IF
00346                            IF( NRT.GT.0 )
00347      $                        CALL DLARTV( NRT, AB( KD-L, J1+L ), INCA,
00348      $                                     AB( KD-L+1, J1+L ), INCA,
00349      $                                     D( J1 ), WORK( J1 ), KD1 )
00350    30                   CONTINUE
00351                      ELSE
00352                         J1END = J1 + KD1*( NR-2 )
00353                         IF( J1END.GE.J1 ) THEN
00354                            DO 40 JIN = J1, J1END, KD1
00355                               CALL DROT( KD-1, AB( KD-1, JIN+1 ), INCX,
00356      $                                   AB( KD, JIN+1 ), INCX,
00357      $                                   D( JIN ), WORK( JIN ) )
00358    40                      CONTINUE
00359                         END IF
00360                         LEND = MIN( KDM1, N-J2 )
00361                         LAST = J1END + KD1
00362                         IF( LEND.GT.0 )
00363      $                     CALL DROT( LEND, AB( KD-1, LAST+1 ), INCX,
00364      $                                AB( KD, LAST+1 ), INCX, D( LAST ),
00365      $                                WORK( LAST ) )
00366                      END IF
00367                   END IF
00368 *
00369                   IF( WANTQ ) THEN
00370 *
00371 *                    accumulate product of plane rotations in Q
00372 *
00373                      IF( INITQ ) THEN
00374 *
00375 *                 take advantage of the fact that Q was
00376 *                 initially the Identity matrix
00377 *
00378                         IQEND = MAX( IQEND, J2 )
00379                         I2 = MAX( 0, K-3 )
00380                         IQAEND = 1 + I*KD
00381                         IF( K.EQ.2 )
00382      $                     IQAEND = IQAEND + KD
00383                         IQAEND = MIN( IQAEND, IQEND )
00384                         DO 50 J = J1, J2, KD1
00385                            IBL = I - I2 / KDM1
00386                            I2 = I2 + 1
00387                            IQB = MAX( 1, J-IBL )
00388                            NQ = 1 + IQAEND - IQB
00389                            IQAEND = MIN( IQAEND+KD, IQEND )
00390                            CALL DROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
00391      $                                1, D( J ), WORK( J ) )
00392    50                   CONTINUE
00393                      ELSE
00394 *
00395                         DO 60 J = J1, J2, KD1
00396                            CALL DROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
00397      $                                D( J ), WORK( J ) )
00398    60                   CONTINUE
00399                      END IF
00400 *
00401                   END IF
00402 *
00403                   IF( J2+KDN.GT.N ) THEN
00404 *
00405 *                    adjust J2 to keep within the bounds of the matrix
00406 *
00407                      NR = NR - 1
00408                      J2 = J2 - KDN - 1
00409                   END IF
00410 *
00411                   DO 70 J = J1, J2, KD1
00412 *
00413 *                    create nonzero element a(j-1,j+kd) outside the band
00414 *                    and store it in WORK
00415 *
00416                      WORK( J+KD ) = WORK( J )*AB( 1, J+KD )
00417                      AB( 1, J+KD ) = D( J )*AB( 1, J+KD )
00418    70             CONTINUE
00419    80          CONTINUE
00420    90       CONTINUE
00421          END IF
00422 *
00423          IF( KD.GT.0 ) THEN
00424 *
00425 *           copy off-diagonal elements to E
00426 *
00427             DO 100 I = 1, N - 1
00428                E( I ) = AB( KD, I+1 )
00429   100       CONTINUE
00430          ELSE
00431 *
00432 *           set E to zero if original matrix was diagonal
00433 *
00434             DO 110 I = 1, N - 1
00435                E( I ) = ZERO
00436   110       CONTINUE
00437          END IF
00438 *
00439 *        copy diagonal elements to D
00440 *
00441          DO 120 I = 1, N
00442             D( I ) = AB( KD1, I )
00443   120    CONTINUE
00444 *
00445       ELSE
00446 *
00447          IF( KD.GT.1 ) THEN
00448 *
00449 *           Reduce to tridiagonal form, working with lower triangle
00450 *
00451             NR = 0
00452             J1 = KDN + 2
00453             J2 = 1
00454 *
00455             DO 210 I = 1, N - 2
00456 *
00457 *              Reduce i-th column of matrix to tridiagonal form
00458 *
00459                DO 200 K = KDN + 1, 2, -1
00460                   J1 = J1 + KDN
00461                   J2 = J2 + KDN
00462 *
00463                   IF( NR.GT.0 ) THEN
00464 *
00465 *                    generate plane rotations to annihilate nonzero
00466 *                    elements which have been created outside the band
00467 *
00468                      CALL DLARGV( NR, AB( KD1, J1-KD1 ), INCA,
00469      $                            WORK( J1 ), KD1, D( J1 ), KD1 )
00470 *
00471 *                    apply plane rotations from one side
00472 *
00473 *
00474 *                    Dependent on the the number of diagonals either
00475 *                    DLARTV or DROT is used
00476 *
00477                      IF( NR.GT.2*KD-1 ) THEN
00478                         DO 130 L = 1, KD - 1
00479                            CALL DLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA,
00480      $                                  AB( KD1-L+1, J1-KD1+L ), INCA,
00481      $                                  D( J1 ), WORK( J1 ), KD1 )
00482   130                   CONTINUE
00483                      ELSE
00484                         JEND = J1 + KD1*( NR-1 )
00485                         DO 140 JINC = J1, JEND, KD1
00486                            CALL DROT( KDM1, AB( KD, JINC-KD ), INCX,
00487      $                                AB( KD1, JINC-KD ), INCX,
00488      $                                D( JINC ), WORK( JINC ) )
00489   140                   CONTINUE
00490                      END IF
00491 *
00492                   END IF
00493 *
00494                   IF( K.GT.2 ) THEN
00495                      IF( K.LE.N-I+1 ) THEN
00496 *
00497 *                       generate plane rotation to annihilate a(i+k-1,i)
00498 *                       within the band
00499 *
00500                         CALL DLARTG( AB( K-1, I ), AB( K, I ),
00501      $                               D( I+K-1 ), WORK( I+K-1 ), TEMP )
00502                         AB( K-1, I ) = TEMP
00503 *
00504 *                       apply rotation from the left
00505 *
00506                         CALL DROT( K-3, AB( K-2, I+1 ), LDAB-1,
00507      $                             AB( K-1, I+1 ), LDAB-1, D( I+K-1 ),
00508      $                             WORK( I+K-1 ) )
00509                      END IF
00510                      NR = NR + 1
00511                      J1 = J1 - KDN - 1
00512                   END IF
00513 *
00514 *                 apply plane rotations from both sides to diagonal
00515 *                 blocks
00516 *
00517                   IF( NR.GT.0 )
00518      $               CALL DLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ),
00519      $                            AB( 2, J1-1 ), INCA, D( J1 ),
00520      $                            WORK( J1 ), KD1 )
00521 *
00522 *                 apply plane rotations from the right
00523 *
00524 *
00525 *                    Dependent on the the number of diagonals either
00526 *                    DLARTV or DROT is used
00527 *
00528                   IF( NR.GT.0 ) THEN
00529                      IF( NR.GT.2*KD-1 ) THEN
00530                         DO 150 L = 1, KD - 1
00531                            IF( J2+L.GT.N ) THEN
00532                               NRT = NR - 1
00533                            ELSE
00534                               NRT = NR
00535                            END IF
00536                            IF( NRT.GT.0 )
00537      $                        CALL DLARTV( NRT, AB( L+2, J1-1 ), INCA,
00538      $                                     AB( L+1, J1 ), INCA, D( J1 ),
00539      $                                     WORK( J1 ), KD1 )
00540   150                   CONTINUE
00541                      ELSE
00542                         J1END = J1 + KD1*( NR-2 )
00543                         IF( J1END.GE.J1 ) THEN
00544                            DO 160 J1INC = J1, J1END, KD1
00545                               CALL DROT( KDM1, AB( 3, J1INC-1 ), 1,
00546      $                                   AB( 2, J1INC ), 1, D( J1INC ),
00547      $                                   WORK( J1INC ) )
00548   160                      CONTINUE
00549                         END IF
00550                         LEND = MIN( KDM1, N-J2 )
00551                         LAST = J1END + KD1
00552                         IF( LEND.GT.0 )
00553      $                     CALL DROT( LEND, AB( 3, LAST-1 ), 1,
00554      $                                AB( 2, LAST ), 1, D( LAST ),
00555      $                                WORK( LAST ) )
00556                      END IF
00557                   END IF
00558 *
00559 *
00560 *
00561                   IF( WANTQ ) THEN
00562 *
00563 *                    accumulate product of plane rotations in Q
00564 *
00565                      IF( INITQ ) THEN
00566 *
00567 *                 take advantage of the fact that Q was
00568 *                 initially the Identity matrix
00569 *
00570                         IQEND = MAX( IQEND, J2 )
00571                         I2 = MAX( 0, K-3 )
00572                         IQAEND = 1 + I*KD
00573                         IF( K.EQ.2 )
00574      $                     IQAEND = IQAEND + KD
00575                         IQAEND = MIN( IQAEND, IQEND )
00576                         DO 170 J = J1, J2, KD1
00577                            IBL = I - I2 / KDM1
00578                            I2 = I2 + 1
00579                            IQB = MAX( 1, J-IBL )
00580                            NQ = 1 + IQAEND - IQB
00581                            IQAEND = MIN( IQAEND+KD, IQEND )
00582                            CALL DROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
00583      $                                1, D( J ), WORK( J ) )
00584   170                   CONTINUE
00585                      ELSE
00586 *
00587                         DO 180 J = J1, J2, KD1
00588                            CALL DROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
00589      $                                D( J ), WORK( J ) )
00590   180                   CONTINUE
00591                      END IF
00592                   END IF
00593 *
00594                   IF( J2+KDN.GT.N ) THEN
00595 *
00596 *                    adjust J2 to keep within the bounds of the matrix
00597 *
00598                      NR = NR - 1
00599                      J2 = J2 - KDN - 1
00600                   END IF
00601 *
00602                   DO 190 J = J1, J2, KD1
00603 *
00604 *                    create nonzero element a(j+kd,j-1) outside the
00605 *                    band and store it in WORK
00606 *
00607                      WORK( J+KD ) = WORK( J )*AB( KD1, J )
00608                      AB( KD1, J ) = D( J )*AB( KD1, J )
00609   190             CONTINUE
00610   200          CONTINUE
00611   210       CONTINUE
00612          END IF
00613 *
00614          IF( KD.GT.0 ) THEN
00615 *
00616 *           copy off-diagonal elements to E
00617 *
00618             DO 220 I = 1, N - 1
00619                E( I ) = AB( 2, I )
00620   220       CONTINUE
00621          ELSE
00622 *
00623 *           set E to zero if original matrix was diagonal
00624 *
00625             DO 230 I = 1, N - 1
00626                E( I ) = ZERO
00627   230       CONTINUE
00628          END IF
00629 *
00630 *        copy diagonal elements to D
00631 *
00632          DO 240 I = 1, N
00633             D( I ) = AB( 1, I )
00634   240    CONTINUE
00635       END IF
00636 *
00637       RETURN
00638 *
00639 *     End of DSBTRD
00640 *
00641       END
 All Files Functions