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