1 /* ./src_f77/dsbtrd.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static doublereal c_b9 = 0.;
11 static doublereal c_b10 = 1.;
12 static integer c__1 = 1;
13
dsbtrd_(char * vect,char * uplo,integer * n,integer * kd,doublereal * ab,integer * ldab,doublereal * d__,doublereal * e,doublereal * q,integer * ldq,doublereal * work,integer * info,ftnlen vect_len,ftnlen uplo_len)14 /* Subroutine */ int dsbtrd_(char *vect, char *uplo, integer *n, integer *kd,
15 doublereal *ab, integer *ldab, doublereal *d__, doublereal *e,
16 doublereal *q, integer *ldq, doublereal *work, integer *info, ftnlen
17 vect_len, ftnlen uplo_len)
18 {
19 /* System generated locals */
20 integer ab_dim1, ab_offset, q_dim1, q_offset, i__1, i__2, i__3, i__4,
21 i__5;
22
23 /* Local variables */
24 static integer i__, j, k, l, i2, j1, j2, nq, nr, kd1, ibl, iqb, kdn, jin,
25 nrt, kdm1, inca, jend, lend, jinc, incx, last;
26 static doublereal temp;
27 extern /* Subroutine */ int drot_(integer *, doublereal *, integer *,
28 doublereal *, integer *, doublereal *, doublereal *);
29 static integer j1end, j1inc, iqend;
30 extern logical lsame_(char *, char *, ftnlen, ftnlen);
31 static logical initq, wantq, upper;
32 extern /* Subroutine */ int dlar2v_(integer *, doublereal *, doublereal *,
33 doublereal *, integer *, doublereal *, doublereal *, integer *);
34 static integer iqaend;
35 extern /* Subroutine */ int dlaset_(char *, integer *, integer *,
36 doublereal *, doublereal *, doublereal *, integer *, ftnlen),
37 dlartg_(doublereal *, doublereal *, doublereal *, doublereal *,
38 doublereal *), xerbla_(char *, integer *, ftnlen), dlargv_(
39 integer *, doublereal *, integer *, doublereal *, integer *,
40 doublereal *, integer *), dlartv_(integer *, doublereal *,
41 integer *, doublereal *, integer *, doublereal *, doublereal *,
42 integer *);
43
44
45 /* -- LAPACK routine (version 3.0) -- */
46 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
47 /* Courant Institute, Argonne National Lab, and Rice University */
48 /* June 30, 1999 */
49
50 /* .. Scalar Arguments .. */
51 /* .. */
52 /* .. Array Arguments .. */
53 /* .. */
54
55 /* Purpose */
56 /* ======= */
57
58 /* DSBTRD reduces a real symmetric band matrix A to symmetric */
59 /* tridiagonal form T by an orthogonal similarity transformation: */
60 /* Q**T * A * Q = T. */
61
62 /* Arguments */
63 /* ========= */
64
65 /* VECT (input) CHARACTER*1 */
66 /* = 'N': do not form Q; */
67 /* = 'V': form Q; */
68 /* = 'U': update a matrix X, by forming X*Q. */
69
70 /* UPLO (input) CHARACTER*1 */
71 /* = 'U': Upper triangle of A is stored; */
72 /* = 'L': Lower triangle of A is stored. */
73
74 /* N (input) INTEGER */
75 /* The order of the matrix A. N >= 0. */
76
77 /* KD (input) INTEGER */
78 /* The number of superdiagonals of the matrix A if UPLO = 'U', */
79 /* or the number of subdiagonals if UPLO = 'L'. KD >= 0. */
80
81 /* AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N) */
82 /* On entry, the upper or lower triangle of the symmetric band */
83 /* matrix A, stored in the first KD+1 rows of the array. The */
84 /* j-th column of A is stored in the j-th column of the array AB */
85 /* as follows: */
86 /* if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; */
87 /* if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). */
88 /* On exit, the diagonal elements of AB are overwritten by the */
89 /* diagonal elements of the tridiagonal matrix T; if KD > 0, the */
90 /* elements on the first superdiagonal (if UPLO = 'U') or the */
91 /* first subdiagonal (if UPLO = 'L') are overwritten by the */
92 /* off-diagonal elements of T; the rest of AB is overwritten by */
93 /* values generated during the reduction. */
94
95 /* LDAB (input) INTEGER */
96 /* The leading dimension of the array AB. LDAB >= KD+1. */
97
98 /* D (output) DOUBLE PRECISION array, dimension (N) */
99 /* The diagonal elements of the tridiagonal matrix T. */
100
101 /* E (output) DOUBLE PRECISION array, dimension (N-1) */
102 /* The off-diagonal elements of the tridiagonal matrix T: */
103 /* E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'. */
104
105 /* Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N) */
106 /* On entry, if VECT = 'U', then Q must contain an N-by-N */
107 /* matrix X; if VECT = 'N' or 'V', then Q need not be set. */
108
109 /* On exit: */
110 /* if VECT = 'V', Q contains the N-by-N orthogonal matrix Q; */
111 /* if VECT = 'U', Q contains the product X*Q; */
112 /* if VECT = 'N', the array Q is not referenced. */
113
114 /* LDQ (input) INTEGER */
115 /* The leading dimension of the array Q. */
116 /* LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'. */
117
118 /* WORK (workspace) DOUBLE PRECISION array, dimension (N) */
119
120 /* INFO (output) INTEGER */
121 /* = 0: successful exit */
122 /* < 0: if INFO = -i, the i-th argument had an illegal value */
123
124 /* Further Details */
125 /* =============== */
126
127 /* Modified by Linda Kaufman, Bell Labs. */
128
129 /* ===================================================================== */
130
131 /* .. Parameters .. */
132 /* .. */
133 /* .. Local Scalars .. */
134 /* .. */
135 /* .. External Subroutines .. */
136 /* .. */
137 /* .. Intrinsic Functions .. */
138 /* .. */
139 /* .. External Functions .. */
140 /* .. */
141 /* .. Executable Statements .. */
142
143 /* Test the input parameters */
144
145 /* Parameter adjustments */
146 ab_dim1 = *ldab;
147 ab_offset = 1 + ab_dim1;
148 ab -= ab_offset;
149 --d__;
150 --e;
151 q_dim1 = *ldq;
152 q_offset = 1 + q_dim1;
153 q -= q_offset;
154 --work;
155
156 /* Function Body */
157 initq = lsame_(vect, "V", (ftnlen)1, (ftnlen)1);
158 wantq = initq || lsame_(vect, "U", (ftnlen)1, (ftnlen)1);
159 upper = lsame_(uplo, "U", (ftnlen)1, (ftnlen)1);
160 kd1 = *kd + 1;
161 kdm1 = *kd - 1;
162 incx = *ldab - 1;
163 iqend = 1;
164
165 *info = 0;
166 if (! wantq && ! lsame_(vect, "N", (ftnlen)1, (ftnlen)1)) {
167 *info = -1;
168 } else if (! upper && ! lsame_(uplo, "L", (ftnlen)1, (ftnlen)1)) {
169 *info = -2;
170 } else if (*n < 0) {
171 *info = -3;
172 } else if (*kd < 0) {
173 *info = -4;
174 } else if (*ldab < kd1) {
175 *info = -6;
176 } else if (*ldq < max(1,*n) && wantq) {
177 *info = -10;
178 }
179 if (*info != 0) {
180 i__1 = -(*info);
181 xerbla_("DSBTRD", &i__1, (ftnlen)6);
182 return 0;
183 }
184
185 /* Quick return if possible */
186
187 if (*n == 0) {
188 return 0;
189 }
190
191 /* Initialize Q to the unit matrix, if needed */
192
193 if (initq) {
194 dlaset_("Full", n, n, &c_b9, &c_b10, &q[q_offset], ldq, (ftnlen)4);
195 }
196
197 /* Wherever possible, plane rotations are generated and applied in */
198 /* vector operations of length NR over the index set J1:J2:KD1. */
199
200 /* The cosines and sines of the plane rotations are stored in the */
201 /* arrays D and WORK. */
202
203 inca = kd1 * *ldab;
204 /* Computing MIN */
205 i__1 = *n - 1;
206 kdn = min(i__1,*kd);
207 if (upper) {
208
209 if (*kd > 1) {
210
211 /* Reduce to tridiagonal form, working with upper triangle */
212
213 nr = 0;
214 j1 = kdn + 2;
215 j2 = 1;
216
217 i__1 = *n - 2;
218 for (i__ = 1; i__ <= i__1; ++i__) {
219
220 /* Reduce i-th row of matrix to tridiagonal form */
221
222 for (k = kdn + 1; k >= 2; --k) {
223 j1 += kdn;
224 j2 += kdn;
225
226 if (nr > 0) {
227
228 /* generate plane rotations to annihilate nonzero */
229 /* elements which have been created outside the band */
230
231 dlargv_(&nr, &ab[(j1 - 1) * ab_dim1 + 1], &inca, &
232 work[j1], &kd1, &d__[j1], &kd1);
233
234 /* apply rotations from the right */
235
236
237 /* Dependent on the the number of diagonals either */
238 /* DLARTV or DROT is used */
239
240 if (nr >= (*kd << 1) - 1) {
241 i__2 = *kd - 1;
242 for (l = 1; l <= i__2; ++l) {
243 dlartv_(&nr, &ab[l + 1 + (j1 - 1) * ab_dim1],
244 &inca, &ab[l + j1 * ab_dim1], &inca, &
245 d__[j1], &work[j1], &kd1);
246 /* L10: */
247 }
248
249 } else {
250 jend = j1 + (nr - 1) * kd1;
251 i__2 = jend;
252 i__3 = kd1;
253 for (jinc = j1; i__3 < 0 ? jinc >= i__2 : jinc <=
254 i__2; jinc += i__3) {
255 drot_(&kdm1, &ab[(jinc - 1) * ab_dim1 + 2], &
256 c__1, &ab[jinc * ab_dim1 + 1], &c__1,
257 &d__[jinc], &work[jinc]);
258 /* L20: */
259 }
260 }
261 }
262
263
264 if (k > 2) {
265 if (k <= *n - i__ + 1) {
266
267 /* generate plane rotation to annihilate a(i,i+k-1) */
268 /* within the band */
269
270 dlartg_(&ab[*kd - k + 3 + (i__ + k - 2) * ab_dim1]
271 , &ab[*kd - k + 2 + (i__ + k - 1) *
272 ab_dim1], &d__[i__ + k - 1], &work[i__ +
273 k - 1], &temp);
274 ab[*kd - k + 3 + (i__ + k - 2) * ab_dim1] = temp;
275
276 /* apply rotation from the right */
277
278 i__3 = k - 3;
279 drot_(&i__3, &ab[*kd - k + 4 + (i__ + k - 2) *
280 ab_dim1], &c__1, &ab[*kd - k + 3 + (i__ +
281 k - 1) * ab_dim1], &c__1, &d__[i__ + k -
282 1], &work[i__ + k - 1]);
283 }
284 ++nr;
285 j1 = j1 - kdn - 1;
286 }
287
288 /* apply plane rotations from both sides to diagonal */
289 /* blocks */
290
291 if (nr > 0) {
292 dlar2v_(&nr, &ab[kd1 + (j1 - 1) * ab_dim1], &ab[kd1 +
293 j1 * ab_dim1], &ab[*kd + j1 * ab_dim1], &inca,
294 &d__[j1], &work[j1], &kd1);
295 }
296
297 /* apply plane rotations from the left */
298
299 if (nr > 0) {
300 if ((*kd << 1) - 1 < nr) {
301
302 /* Dependent on the the number of diagonals either */
303 /* DLARTV or DROT is used */
304
305 i__3 = *kd - 1;
306 for (l = 1; l <= i__3; ++l) {
307 if (j2 + l > *n) {
308 nrt = nr - 1;
309 } else {
310 nrt = nr;
311 }
312 if (nrt > 0) {
313 dlartv_(&nrt, &ab[*kd - l + (j1 + l) *
314 ab_dim1], &inca, &ab[*kd - l + 1
315 + (j1 + l) * ab_dim1], &inca, &
316 d__[j1], &work[j1], &kd1);
317 }
318 /* L30: */
319 }
320 } else {
321 j1end = j1 + kd1 * (nr - 2);
322 if (j1end >= j1) {
323 i__3 = j1end;
324 i__2 = kd1;
325 for (jin = j1; i__2 < 0 ? jin >= i__3 : jin <=
326 i__3; jin += i__2) {
327 i__4 = *kd - 1;
328 drot_(&i__4, &ab[*kd - 1 + (jin + 1) *
329 ab_dim1], &incx, &ab[*kd + (jin +
330 1) * ab_dim1], &incx, &d__[jin], &
331 work[jin]);
332 /* L40: */
333 }
334 }
335 /* Computing MIN */
336 i__2 = kdm1, i__3 = *n - j2;
337 lend = min(i__2,i__3);
338 last = j1end + kd1;
339 if (lend > 0) {
340 drot_(&lend, &ab[*kd - 1 + (last + 1) *
341 ab_dim1], &incx, &ab[*kd + (last + 1)
342 * ab_dim1], &incx, &d__[last], &work[
343 last]);
344 }
345 }
346 }
347
348 if (wantq) {
349
350 /* accumulate product of plane rotations in Q */
351
352 if (initq) {
353
354 /* take advantage of the fact that Q was */
355 /* initially the Identity matrix */
356
357 iqend = max(iqend,j2);
358 /* Computing MAX */
359 i__2 = 0, i__3 = k - 3;
360 i2 = max(i__2,i__3);
361 iqaend = i__ * *kd + 1;
362 if (k == 2) {
363 iqaend += *kd;
364 }
365 iqaend = min(iqaend,iqend);
366 i__2 = j2;
367 i__3 = kd1;
368 for (j = j1; i__3 < 0 ? j >= i__2 : j <= i__2; j
369 += i__3) {
370 ibl = i__ - i2 / kdm1;
371 ++i2;
372 /* Computing MAX */
373 i__4 = 1, i__5 = j - ibl;
374 iqb = max(i__4,i__5);
375 nq = iqaend + 1 - iqb;
376 /* Computing MIN */
377 i__4 = iqaend + *kd;
378 iqaend = min(i__4,iqend);
379 drot_(&nq, &q[iqb + (j - 1) * q_dim1], &c__1,
380 &q[iqb + j * q_dim1], &c__1, &d__[j],
381 &work[j]);
382 /* L50: */
383 }
384 } else {
385
386 i__3 = j2;
387 i__2 = kd1;
388 for (j = j1; i__2 < 0 ? j >= i__3 : j <= i__3; j
389 += i__2) {
390 drot_(n, &q[(j - 1) * q_dim1 + 1], &c__1, &q[
391 j * q_dim1 + 1], &c__1, &d__[j], &
392 work[j]);
393 /* L60: */
394 }
395 }
396
397 }
398
399 if (j2 + kdn > *n) {
400
401 /* adjust J2 to keep within the bounds of the matrix */
402
403 --nr;
404 j2 = j2 - kdn - 1;
405 }
406
407 i__2 = j2;
408 i__3 = kd1;
409 for (j = j1; i__3 < 0 ? j >= i__2 : j <= i__2; j += i__3)
410 {
411
412 /* create nonzero element a(j-1,j+kd) outside the band */
413 /* and store it in WORK */
414
415 work[j + *kd] = work[j] * ab[(j + *kd) * ab_dim1 + 1];
416 ab[(j + *kd) * ab_dim1 + 1] = d__[j] * ab[(j + *kd) *
417 ab_dim1 + 1];
418 /* L70: */
419 }
420 /* L80: */
421 }
422 /* L90: */
423 }
424 }
425
426 if (*kd > 0) {
427
428 /* copy off-diagonal elements to E */
429
430 i__1 = *n - 1;
431 for (i__ = 1; i__ <= i__1; ++i__) {
432 e[i__] = ab[*kd + (i__ + 1) * ab_dim1];
433 /* L100: */
434 }
435 } else {
436
437 /* set E to zero if original matrix was diagonal */
438
439 i__1 = *n - 1;
440 for (i__ = 1; i__ <= i__1; ++i__) {
441 e[i__] = 0.;
442 /* L110: */
443 }
444 }
445
446 /* copy diagonal elements to D */
447
448 i__1 = *n;
449 for (i__ = 1; i__ <= i__1; ++i__) {
450 d__[i__] = ab[kd1 + i__ * ab_dim1];
451 /* L120: */
452 }
453
454 } else {
455
456 if (*kd > 1) {
457
458 /* Reduce to tridiagonal form, working with lower triangle */
459
460 nr = 0;
461 j1 = kdn + 2;
462 j2 = 1;
463
464 i__1 = *n - 2;
465 for (i__ = 1; i__ <= i__1; ++i__) {
466
467 /* Reduce i-th column of matrix to tridiagonal form */
468
469 for (k = kdn + 1; k >= 2; --k) {
470 j1 += kdn;
471 j2 += kdn;
472
473 if (nr > 0) {
474
475 /* generate plane rotations to annihilate nonzero */
476 /* elements which have been created outside the band */
477
478 dlargv_(&nr, &ab[kd1 + (j1 - kd1) * ab_dim1], &inca, &
479 work[j1], &kd1, &d__[j1], &kd1);
480
481 /* apply plane rotations from one side */
482
483
484 /* Dependent on the the number of diagonals either */
485 /* DLARTV or DROT is used */
486
487 if (nr > (*kd << 1) - 1) {
488 i__3 = *kd - 1;
489 for (l = 1; l <= i__3; ++l) {
490 dlartv_(&nr, &ab[kd1 - l + (j1 - kd1 + l) *
491 ab_dim1], &inca, &ab[kd1 - l + 1 + (
492 j1 - kd1 + l) * ab_dim1], &inca, &d__[
493 j1], &work[j1], &kd1);
494 /* L130: */
495 }
496 } else {
497 jend = j1 + kd1 * (nr - 1);
498 i__3 = jend;
499 i__2 = kd1;
500 for (jinc = j1; i__2 < 0 ? jinc >= i__3 : jinc <=
501 i__3; jinc += i__2) {
502 drot_(&kdm1, &ab[*kd + (jinc - *kd) * ab_dim1]
503 , &incx, &ab[kd1 + (jinc - *kd) *
504 ab_dim1], &incx, &d__[jinc], &work[
505 jinc]);
506 /* L140: */
507 }
508 }
509
510 }
511
512 if (k > 2) {
513 if (k <= *n - i__ + 1) {
514
515 /* generate plane rotation to annihilate a(i+k-1,i) */
516 /* within the band */
517
518 dlartg_(&ab[k - 1 + i__ * ab_dim1], &ab[k + i__ *
519 ab_dim1], &d__[i__ + k - 1], &work[i__ +
520 k - 1], &temp);
521 ab[k - 1 + i__ * ab_dim1] = temp;
522
523 /* apply rotation from the left */
524
525 i__2 = k - 3;
526 i__3 = *ldab - 1;
527 i__4 = *ldab - 1;
528 drot_(&i__2, &ab[k - 2 + (i__ + 1) * ab_dim1], &
529 i__3, &ab[k - 1 + (i__ + 1) * ab_dim1], &
530 i__4, &d__[i__ + k - 1], &work[i__ + k -
531 1]);
532 }
533 ++nr;
534 j1 = j1 - kdn - 1;
535 }
536
537 /* apply plane rotations from both sides to diagonal */
538 /* blocks */
539
540 if (nr > 0) {
541 dlar2v_(&nr, &ab[(j1 - 1) * ab_dim1 + 1], &ab[j1 *
542 ab_dim1 + 1], &ab[(j1 - 1) * ab_dim1 + 2], &
543 inca, &d__[j1], &work[j1], &kd1);
544 }
545
546 /* apply plane rotations from the right */
547
548
549 /* Dependent on the the number of diagonals either */
550 /* DLARTV or DROT is used */
551
552 if (nr > 0) {
553 if (nr > (*kd << 1) - 1) {
554 i__2 = *kd - 1;
555 for (l = 1; l <= i__2; ++l) {
556 if (j2 + l > *n) {
557 nrt = nr - 1;
558 } else {
559 nrt = nr;
560 }
561 if (nrt > 0) {
562 dlartv_(&nrt, &ab[l + 2 + (j1 - 1) *
563 ab_dim1], &inca, &ab[l + 1 + j1 *
564 ab_dim1], &inca, &d__[j1], &work[
565 j1], &kd1);
566 }
567 /* L150: */
568 }
569 } else {
570 j1end = j1 + kd1 * (nr - 2);
571 if (j1end >= j1) {
572 i__2 = j1end;
573 i__3 = kd1;
574 for (j1inc = j1; i__3 < 0 ? j1inc >= i__2 :
575 j1inc <= i__2; j1inc += i__3) {
576 drot_(&kdm1, &ab[(j1inc - 1) * ab_dim1 +
577 3], &c__1, &ab[j1inc * ab_dim1 +
578 2], &c__1, &d__[j1inc], &work[
579 j1inc]);
580 /* L160: */
581 }
582 }
583 /* Computing MIN */
584 i__3 = kdm1, i__2 = *n - j2;
585 lend = min(i__3,i__2);
586 last = j1end + kd1;
587 if (lend > 0) {
588 drot_(&lend, &ab[(last - 1) * ab_dim1 + 3], &
589 c__1, &ab[last * ab_dim1 + 2], &c__1,
590 &d__[last], &work[last]);
591 }
592 }
593 }
594
595
596
597 if (wantq) {
598
599 /* accumulate product of plane rotations in Q */
600
601 if (initq) {
602
603 /* take advantage of the fact that Q was */
604 /* initially the Identity matrix */
605
606 iqend = max(iqend,j2);
607 /* Computing MAX */
608 i__3 = 0, i__2 = k - 3;
609 i2 = max(i__3,i__2);
610 iqaend = i__ * *kd + 1;
611 if (k == 2) {
612 iqaend += *kd;
613 }
614 iqaend = min(iqaend,iqend);
615 i__3 = j2;
616 i__2 = kd1;
617 for (j = j1; i__2 < 0 ? j >= i__3 : j <= i__3; j
618 += i__2) {
619 ibl = i__ - i2 / kdm1;
620 ++i2;
621 /* Computing MAX */
622 i__4 = 1, i__5 = j - ibl;
623 iqb = max(i__4,i__5);
624 nq = iqaend + 1 - iqb;
625 /* Computing MIN */
626 i__4 = iqaend + *kd;
627 iqaend = min(i__4,iqend);
628 drot_(&nq, &q[iqb + (j - 1) * q_dim1], &c__1,
629 &q[iqb + j * q_dim1], &c__1, &d__[j],
630 &work[j]);
631 /* L170: */
632 }
633 } else {
634
635 i__2 = j2;
636 i__3 = kd1;
637 for (j = j1; i__3 < 0 ? j >= i__2 : j <= i__2; j
638 += i__3) {
639 drot_(n, &q[(j - 1) * q_dim1 + 1], &c__1, &q[
640 j * q_dim1 + 1], &c__1, &d__[j], &
641 work[j]);
642 /* L180: */
643 }
644 }
645 }
646
647 if (j2 + kdn > *n) {
648
649 /* adjust J2 to keep within the bounds of the matrix */
650
651 --nr;
652 j2 = j2 - kdn - 1;
653 }
654
655 i__3 = j2;
656 i__2 = kd1;
657 for (j = j1; i__2 < 0 ? j >= i__3 : j <= i__3; j += i__2)
658 {
659
660 /* create nonzero element a(j+kd,j-1) outside the */
661 /* band and store it in WORK */
662
663 work[j + *kd] = work[j] * ab[kd1 + j * ab_dim1];
664 ab[kd1 + j * ab_dim1] = d__[j] * ab[kd1 + j * ab_dim1]
665 ;
666 /* L190: */
667 }
668 /* L200: */
669 }
670 /* L210: */
671 }
672 }
673
674 if (*kd > 0) {
675
676 /* copy off-diagonal elements to E */
677
678 i__1 = *n - 1;
679 for (i__ = 1; i__ <= i__1; ++i__) {
680 e[i__] = ab[i__ * ab_dim1 + 2];
681 /* L220: */
682 }
683 } else {
684
685 /* set E to zero if original matrix was diagonal */
686
687 i__1 = *n - 1;
688 for (i__ = 1; i__ <= i__1; ++i__) {
689 e[i__] = 0.;
690 /* L230: */
691 }
692 }
693
694 /* copy diagonal elements to D */
695
696 i__1 = *n;
697 for (i__ = 1; i__ <= i__1; ++i__) {
698 d__[i__] = ab[i__ * ab_dim1 + 1];
699 /* L240: */
700 }
701 }
702
703 return 0;
704
705 /* End of DSBTRD */
706
707 } /* dsbtrd_ */
708
709