![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SSBTRD 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SSBTRD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssbtrd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssbtrd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssbtrd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SSBTRD( 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 * REAL AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ), 00030 * $ WORK( * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> SSBTRD 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 REAL 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 REAL array, dimension (N) 00101 *> The diagonal elements of the tridiagonal matrix T. 00102 *> \endverbatim 00103 *> 00104 *> \param[out] E 00105 *> \verbatim 00106 *> E is REAL 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 REAL 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 REAL 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 realOTHERcomputational 00153 * 00154 *> \par Further Details: 00155 * ===================== 00156 *> 00157 *> \verbatim 00158 *> 00159 *> Modified by Linda Kaufman, Bell Labs. 00160 *> \endverbatim 00161 *> 00162 * ===================================================================== 00163 SUBROUTINE SSBTRD( 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 REAL AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ), 00177 $ WORK( * ) 00178 * .. 00179 * 00180 * ===================================================================== 00181 * 00182 * .. Parameters .. 00183 REAL ZERO, ONE 00184 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+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 REAL TEMP 00192 * .. 00193 * .. External Subroutines .. 00194 EXTERNAL SLAR2V, SLARGV, SLARTG, SLARTV, SLASET, SROT, 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( 'SSBTRD', -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 SLASET( '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 SLARGV( 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 * SLARTV or SROT is used 00284 * 00285 IF( NR.GE.2*KD-1 ) THEN 00286 DO 10 L = 1, KD - 1 00287 CALL SLARTV( 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 SROT( 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 SLARTG( 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 SROT( 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 SLAR2V( 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 * SLARTV or SROT 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 SLARTV( 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 SROT( 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 SROT( 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 SROT( 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 SROT( 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 SLARGV( 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 * SLARTV or SROT is used 00476 * 00477 IF( NR.GT.2*KD-1 ) THEN 00478 DO 130 L = 1, KD - 1 00479 CALL SLARTV( 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 SROT( 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 SLARTG( 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 SROT( 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 SLAR2V( 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 * SLARTV or SROT 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 SLARTV( 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 SROT( 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 SROT( 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 SROT( 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 SROT( 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 SSBTRD 00640 * 00641 END