1 /* ./src_f77/ssbgst.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 real c_b8 = 0.f;
11 static real c_b9 = 1.f;
12 static integer c__1 = 1;
13 static real c_b20 = -1.f;
14
ssbgst_(char * vect,char * uplo,integer * n,integer * ka,integer * kb,real * ab,integer * ldab,real * bb,integer * ldbb,real * x,integer * ldx,real * work,integer * info,ftnlen vect_len,ftnlen uplo_len)15 /* Subroutine */ int ssbgst_(char *vect, char *uplo, integer *n, integer *ka,
16 integer *kb, real *ab, integer *ldab, real *bb, integer *ldbb, real *
17 x, integer *ldx, real *work, integer *info, ftnlen vect_len, ftnlen
18 uplo_len)
19 {
20 /* System generated locals */
21 integer ab_dim1, ab_offset, bb_dim1, bb_offset, x_dim1, x_offset, i__1,
22 i__2, i__3, i__4;
23 real r__1;
24
25 /* Local variables */
26 static integer i__, j, k, l, m;
27 static real t;
28 static integer i0, i1, i2, j1, j2;
29 static real ra;
30 static integer nr, nx, ka1, kb1;
31 static real ra1;
32 static integer j1t, j2t;
33 static real bii;
34 static integer kbt, nrt, inca;
35 extern /* Subroutine */ int sger_(integer *, integer *, real *, real *,
36 integer *, real *, integer *, real *, integer *), srot_(integer *,
37 real *, integer *, real *, integer *, real *, real *);
38 extern logical lsame_(char *, char *, ftnlen, ftnlen);
39 extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
40 static logical upper, wantx;
41 extern /* Subroutine */ int slar2v_(integer *, real *, real *, real *,
42 integer *, real *, real *, integer *), xerbla_(char *, integer *,
43 ftnlen);
44 static logical update;
45 extern /* Subroutine */ int slaset_(char *, integer *, integer *, real *,
46 real *, real *, integer *, ftnlen), slartg_(real *, real *, real *
47 , real *, real *), slargv_(integer *, real *, integer *, real *,
48 integer *, real *, integer *), slartv_(integer *, real *, integer
49 *, real *, integer *, real *, real *, integer *);
50
51
52 /* -- LAPACK routine (version 3.0) -- */
53 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
54 /* Courant Institute, Argonne National Lab, and Rice University */
55 /* June 30, 1999 */
56
57 /* .. Scalar Arguments .. */
58 /* .. */
59 /* .. Array Arguments .. */
60 /* .. */
61
62 /* Purpose */
63 /* ======= */
64
65 /* SSBGST reduces a real symmetric-definite banded generalized */
66 /* eigenproblem A*x = lambda*B*x to standard form C*y = lambda*y, */
67 /* such that C has the same bandwidth as A. */
68
69 /* B must have been previously factorized as S**T*S by SPBSTF, using a */
70 /* split Cholesky factorization. A is overwritten by C = X**T*A*X, where */
71 /* X = S**(-1)*Q and Q is an orthogonal matrix chosen to preserve the */
72 /* bandwidth of A. */
73
74 /* Arguments */
75 /* ========= */
76
77 /* VECT (input) CHARACTER*1 */
78 /* = 'N': do not form the transformation matrix X; */
79 /* = 'V': form X. */
80
81 /* UPLO (input) CHARACTER*1 */
82 /* = 'U': Upper triangle of A is stored; */
83 /* = 'L': Lower triangle of A is stored. */
84
85 /* N (input) INTEGER */
86 /* The order of the matrices A and B. N >= 0. */
87
88 /* KA (input) INTEGER */
89 /* The number of superdiagonals of the matrix A if UPLO = 'U', */
90 /* or the number of subdiagonals if UPLO = 'L'. KA >= 0. */
91
92 /* KB (input) INTEGER */
93 /* The number of superdiagonals of the matrix B if UPLO = 'U', */
94 /* or the number of subdiagonals if UPLO = 'L'. KA >= KB >= 0. */
95
96 /* AB (input/output) REAL array, dimension (LDAB,N) */
97 /* On entry, the upper or lower triangle of the symmetric band */
98 /* matrix A, stored in the first ka+1 rows of the array. The */
99 /* j-th column of A is stored in the j-th column of the array AB */
100 /* as follows: */
101 /* if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j; */
102 /* if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+ka). */
103
104 /* On exit, the transformed matrix X**T*A*X, stored in the same */
105 /* format as A. */
106
107 /* LDAB (input) INTEGER */
108 /* The leading dimension of the array AB. LDAB >= KA+1. */
109
110 /* BB (input) REAL array, dimension (LDBB,N) */
111 /* The banded factor S from the split Cholesky factorization of */
112 /* B, as returned by SPBSTF, stored in the first KB+1 rows of */
113 /* the array. */
114
115 /* LDBB (input) INTEGER */
116 /* The leading dimension of the array BB. LDBB >= KB+1. */
117
118 /* X (output) REAL array, dimension (LDX,N) */
119 /* If VECT = 'V', the n-by-n matrix X. */
120 /* If VECT = 'N', the array X is not referenced. */
121
122 /* LDX (input) INTEGER */
123 /* The leading dimension of the array X. */
124 /* LDX >= max(1,N) if VECT = 'V'; LDX >= 1 otherwise. */
125
126 /* WORK (workspace) REAL array, dimension (2*N) */
127
128 /* INFO (output) INTEGER */
129 /* = 0: successful exit */
130 /* < 0: if INFO = -i, the i-th argument had an illegal value. */
131
132 /* ===================================================================== */
133
134 /* .. Parameters .. */
135 /* .. */
136 /* .. Local Scalars .. */
137 /* .. */
138 /* .. External Functions .. */
139 /* .. */
140 /* .. External Subroutines .. */
141 /* .. */
142 /* .. Intrinsic Functions .. */
143 /* .. */
144 /* .. Executable Statements .. */
145
146 /* Test the input parameters */
147
148 /* Parameter adjustments */
149 ab_dim1 = *ldab;
150 ab_offset = 1 + ab_dim1;
151 ab -= ab_offset;
152 bb_dim1 = *ldbb;
153 bb_offset = 1 + bb_dim1;
154 bb -= bb_offset;
155 x_dim1 = *ldx;
156 x_offset = 1 + x_dim1;
157 x -= x_offset;
158 --work;
159
160 /* Function Body */
161 wantx = lsame_(vect, "V", (ftnlen)1, (ftnlen)1);
162 upper = lsame_(uplo, "U", (ftnlen)1, (ftnlen)1);
163 ka1 = *ka + 1;
164 kb1 = *kb + 1;
165 *info = 0;
166 if (! wantx && ! 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 (*ka < 0) {
173 *info = -4;
174 } else if (*kb < 0) {
175 *info = -5;
176 } else if (*ldab < *ka + 1) {
177 *info = -7;
178 } else if (*ldbb < *kb + 1) {
179 *info = -9;
180 } else if (*ldx < 1 || wantx && *ldx < max(1,*n)) {
181 *info = -11;
182 }
183 if (*info != 0) {
184 i__1 = -(*info);
185 xerbla_("SSBGST", &i__1, (ftnlen)6);
186 return 0;
187 }
188
189 /* Quick return if possible */
190
191 if (*n == 0) {
192 return 0;
193 }
194
195 inca = *ldab * ka1;
196
197 /* Initialize X to the unit matrix, if needed */
198
199 if (wantx) {
200 slaset_("Full", n, n, &c_b8, &c_b9, &x[x_offset], ldx, (ftnlen)4);
201 }
202
203 /* Set M to the splitting point m. It must be the same value as is */
204 /* used in SPBSTF. The chosen value allows the arrays WORK and RWORK */
205 /* to be of dimension (N). */
206
207 m = (*n + *kb) / 2;
208
209 /* The routine works in two phases, corresponding to the two halves */
210 /* of the split Cholesky factorization of B as S**T*S where */
211
212 /* S = ( U ) */
213 /* ( M L ) */
214
215 /* with U upper triangular of order m, and L lower triangular of */
216 /* order n-m. S has the same bandwidth as B. */
217
218 /* S is treated as a product of elementary matrices: */
219
220 /* S = S(m)*S(m-1)*...*S(2)*S(1)*S(m+1)*S(m+2)*...*S(n-1)*S(n) */
221
222 /* where S(i) is determined by the i-th row of S. */
223
224 /* In phase 1, the index i takes the values n, n-1, ... , m+1; */
225 /* in phase 2, it takes the values 1, 2, ... , m. */
226
227 /* For each value of i, the current matrix A is updated by forming */
228 /* inv(S(i))**T*A*inv(S(i)). This creates a triangular bulge outside */
229 /* the band of A. The bulge is then pushed down toward the bottom of */
230 /* A in phase 1, and up toward the top of A in phase 2, by applying */
231 /* plane rotations. */
232
233 /* There are kb*(kb+1)/2 elements in the bulge, but at most 2*kb-1 */
234 /* of them are linearly independent, so annihilating a bulge requires */
235 /* only 2*kb-1 plane rotations. The rotations are divided into a 1st */
236 /* set of kb-1 rotations, and a 2nd set of kb rotations. */
237
238 /* Wherever possible, rotations are generated and applied in vector */
239 /* operations of length NR between the indices J1 and J2 (sometimes */
240 /* replaced by modified values NRT, J1T or J2T). */
241
242 /* The cosines and sines of the rotations are stored in the array */
243 /* WORK. The cosines of the 1st set of rotations are stored in */
244 /* elements n+2:n+m-kb-1 and the sines of the 1st set in elements */
245 /* 2:m-kb-1; the cosines of the 2nd set are stored in elements */
246 /* n+m-kb+1:2*n and the sines of the second set in elements m-kb+1:n. */
247
248 /* The bulges are not formed explicitly; nonzero elements outside the */
249 /* band are created only when they are required for generating new */
250 /* rotations; they are stored in the array WORK, in positions where */
251 /* they are later overwritten by the sines of the rotations which */
252 /* annihilate them. */
253
254 /* **************************** Phase 1 ***************************** */
255
256 /* The logical structure of this phase is: */
257
258 /* UPDATE = .TRUE. */
259 /* DO I = N, M + 1, -1 */
260 /* use S(i) to update A and create a new bulge */
261 /* apply rotations to push all bulges KA positions downward */
262 /* END DO */
263 /* UPDATE = .FALSE. */
264 /* DO I = M + KA + 1, N - 1 */
265 /* apply rotations to push all bulges KA positions downward */
266 /* END DO */
267
268 /* To avoid duplicating code, the two loops are merged. */
269
270 update = TRUE_;
271 i__ = *n + 1;
272 L10:
273 if (update) {
274 --i__;
275 /* Computing MIN */
276 i__1 = *kb, i__2 = i__ - 1;
277 kbt = min(i__1,i__2);
278 i0 = i__ - 1;
279 /* Computing MIN */
280 i__1 = *n, i__2 = i__ + *ka;
281 i1 = min(i__1,i__2);
282 i2 = i__ - kbt + ka1;
283 if (i__ < m + 1) {
284 update = FALSE_;
285 ++i__;
286 i0 = m;
287 if (*ka == 0) {
288 goto L480;
289 }
290 goto L10;
291 }
292 } else {
293 i__ += *ka;
294 if (i__ > *n - 1) {
295 goto L480;
296 }
297 }
298
299 if (upper) {
300
301 /* Transform A, working with the upper triangle */
302
303 if (update) {
304
305 /* Form inv(S(i))**T * A * inv(S(i)) */
306
307 bii = bb[kb1 + i__ * bb_dim1];
308 i__1 = i1;
309 for (j = i__; j <= i__1; ++j) {
310 ab[i__ - j + ka1 + j * ab_dim1] /= bii;
311 /* L20: */
312 }
313 /* Computing MAX */
314 i__1 = 1, i__2 = i__ - *ka;
315 i__3 = i__;
316 for (j = max(i__1,i__2); j <= i__3; ++j) {
317 ab[j - i__ + ka1 + i__ * ab_dim1] /= bii;
318 /* L30: */
319 }
320 i__3 = i__ - 1;
321 for (k = i__ - kbt; k <= i__3; ++k) {
322 i__1 = k;
323 for (j = i__ - kbt; j <= i__1; ++j) {
324 ab[j - k + ka1 + k * ab_dim1] = ab[j - k + ka1 + k *
325 ab_dim1] - bb[j - i__ + kb1 + i__ * bb_dim1] * ab[
326 k - i__ + ka1 + i__ * ab_dim1] - bb[k - i__ + kb1
327 + i__ * bb_dim1] * ab[j - i__ + ka1 + i__ *
328 ab_dim1] + ab[ka1 + i__ * ab_dim1] * bb[j - i__ +
329 kb1 + i__ * bb_dim1] * bb[k - i__ + kb1 + i__ *
330 bb_dim1];
331 /* L40: */
332 }
333 /* Computing MAX */
334 i__1 = 1, i__2 = i__ - *ka;
335 i__4 = i__ - kbt - 1;
336 for (j = max(i__1,i__2); j <= i__4; ++j) {
337 ab[j - k + ka1 + k * ab_dim1] -= bb[k - i__ + kb1 + i__ *
338 bb_dim1] * ab[j - i__ + ka1 + i__ * ab_dim1];
339 /* L50: */
340 }
341 /* L60: */
342 }
343 i__3 = i1;
344 for (j = i__; j <= i__3; ++j) {
345 /* Computing MAX */
346 i__4 = j - *ka, i__1 = i__ - kbt;
347 i__2 = i__ - 1;
348 for (k = max(i__4,i__1); k <= i__2; ++k) {
349 ab[k - j + ka1 + j * ab_dim1] -= bb[k - i__ + kb1 + i__ *
350 bb_dim1] * ab[i__ - j + ka1 + j * ab_dim1];
351 /* L70: */
352 }
353 /* L80: */
354 }
355
356 if (wantx) {
357
358 /* post-multiply X by inv(S(i)) */
359
360 i__3 = *n - m;
361 r__1 = 1.f / bii;
362 sscal_(&i__3, &r__1, &x[m + 1 + i__ * x_dim1], &c__1);
363 if (kbt > 0) {
364 i__3 = *n - m;
365 sger_(&i__3, &kbt, &c_b20, &x[m + 1 + i__ * x_dim1], &
366 c__1, &bb[kb1 - kbt + i__ * bb_dim1], &c__1, &x[m
367 + 1 + (i__ - kbt) * x_dim1], ldx);
368 }
369 }
370
371 /* store a(i,i1) in RA1 for use in next loop over K */
372
373 ra1 = ab[i__ - i1 + ka1 + i1 * ab_dim1];
374 }
375
376 /* Generate and apply vectors of rotations to chase all the */
377 /* existing bulges KA positions down toward the bottom of the */
378 /* band */
379
380 i__3 = *kb - 1;
381 for (k = 1; k <= i__3; ++k) {
382 if (update) {
383
384 /* Determine the rotations which would annihilate the bulge */
385 /* which has in theory just been created */
386
387 if (i__ - k + *ka < *n && i__ - k > 1) {
388
389 /* generate rotation to annihilate a(i,i-k+ka+1) */
390
391 slartg_(&ab[k + 1 + (i__ - k + *ka) * ab_dim1], &ra1, &
392 work[*n + i__ - k + *ka - m], &work[i__ - k + *ka
393 - m], &ra);
394
395 /* create nonzero element a(i-k,i-k+ka+1) outside the */
396 /* band and store it in WORK(i-k) */
397
398 t = -bb[kb1 - k + i__ * bb_dim1] * ra1;
399 work[i__ - k] = work[*n + i__ - k + *ka - m] * t - work[
400 i__ - k + *ka - m] * ab[(i__ - k + *ka) * ab_dim1
401 + 1];
402 ab[(i__ - k + *ka) * ab_dim1 + 1] = work[i__ - k + *ka -
403 m] * t + work[*n + i__ - k + *ka - m] * ab[(i__ -
404 k + *ka) * ab_dim1 + 1];
405 ra1 = ra;
406 }
407 }
408 /* Computing MAX */
409 i__2 = 1, i__4 = k - i0 + 2;
410 j2 = i__ - k - 1 + max(i__2,i__4) * ka1;
411 nr = (*n - j2 + *ka) / ka1;
412 j1 = j2 + (nr - 1) * ka1;
413 if (update) {
414 /* Computing MAX */
415 i__2 = j2, i__4 = i__ + (*ka << 1) - k + 1;
416 j2t = max(i__2,i__4);
417 } else {
418 j2t = j2;
419 }
420 nrt = (*n - j2t + *ka) / ka1;
421 i__2 = j1;
422 i__4 = ka1;
423 for (j = j2t; i__4 < 0 ? j >= i__2 : j <= i__2; j += i__4) {
424
425 /* create nonzero element a(j-ka,j+1) outside the band */
426 /* and store it in WORK(j-m) */
427
428 work[j - m] *= ab[(j + 1) * ab_dim1 + 1];
429 ab[(j + 1) * ab_dim1 + 1] = work[*n + j - m] * ab[(j + 1) *
430 ab_dim1 + 1];
431 /* L90: */
432 }
433
434 /* generate rotations in 1st set to annihilate elements which */
435 /* have been created outside the band */
436
437 if (nrt > 0) {
438 slargv_(&nrt, &ab[j2t * ab_dim1 + 1], &inca, &work[j2t - m], &
439 ka1, &work[*n + j2t - m], &ka1);
440 }
441 if (nr > 0) {
442
443 /* apply rotations in 1st set from the right */
444
445 i__4 = *ka - 1;
446 for (l = 1; l <= i__4; ++l) {
447 slartv_(&nr, &ab[ka1 - l + j2 * ab_dim1], &inca, &ab[*ka
448 - l + (j2 + 1) * ab_dim1], &inca, &work[*n + j2 -
449 m], &work[j2 - m], &ka1);
450 /* L100: */
451 }
452
453 /* apply rotations in 1st set from both sides to diagonal */
454 /* blocks */
455
456 slar2v_(&nr, &ab[ka1 + j2 * ab_dim1], &ab[ka1 + (j2 + 1) *
457 ab_dim1], &ab[*ka + (j2 + 1) * ab_dim1], &inca, &work[
458 *n + j2 - m], &work[j2 - m], &ka1);
459
460 }
461
462 /* start applying rotations in 1st set from the left */
463
464 i__4 = *kb - k + 1;
465 for (l = *ka - 1; l >= i__4; --l) {
466 nrt = (*n - j2 + l) / ka1;
467 if (nrt > 0) {
468 slartv_(&nrt, &ab[l + (j2 + ka1 - l) * ab_dim1], &inca, &
469 ab[l + 1 + (j2 + ka1 - l) * ab_dim1], &inca, &
470 work[*n + j2 - m], &work[j2 - m], &ka1);
471 }
472 /* L110: */
473 }
474
475 if (wantx) {
476
477 /* post-multiply X by product of rotations in 1st set */
478
479 i__4 = j1;
480 i__2 = ka1;
481 for (j = j2; i__2 < 0 ? j >= i__4 : j <= i__4; j += i__2) {
482 i__1 = *n - m;
483 srot_(&i__1, &x[m + 1 + j * x_dim1], &c__1, &x[m + 1 + (j
484 + 1) * x_dim1], &c__1, &work[*n + j - m], &work[j
485 - m]);
486 /* L120: */
487 }
488 }
489 /* L130: */
490 }
491
492 if (update) {
493 if (i2 <= *n && kbt > 0) {
494
495 /* create nonzero element a(i-kbt,i-kbt+ka+1) outside the */
496 /* band and store it in WORK(i-kbt) */
497
498 work[i__ - kbt] = -bb[kb1 - kbt + i__ * bb_dim1] * ra1;
499 }
500 }
501
502 for (k = *kb; k >= 1; --k) {
503 if (update) {
504 /* Computing MAX */
505 i__3 = 2, i__2 = k - i0 + 1;
506 j2 = i__ - k - 1 + max(i__3,i__2) * ka1;
507 } else {
508 /* Computing MAX */
509 i__3 = 1, i__2 = k - i0 + 1;
510 j2 = i__ - k - 1 + max(i__3,i__2) * ka1;
511 }
512
513 /* finish applying rotations in 2nd set from the left */
514
515 for (l = *kb - k; l >= 1; --l) {
516 nrt = (*n - j2 + *ka + l) / ka1;
517 if (nrt > 0) {
518 slartv_(&nrt, &ab[l + (j2 - l + 1) * ab_dim1], &inca, &ab[
519 l + 1 + (j2 - l + 1) * ab_dim1], &inca, &work[*n
520 + j2 - *ka], &work[j2 - *ka], &ka1);
521 }
522 /* L140: */
523 }
524 nr = (*n - j2 + *ka) / ka1;
525 j1 = j2 + (nr - 1) * ka1;
526 i__3 = j2;
527 i__2 = -ka1;
528 for (j = j1; i__2 < 0 ? j >= i__3 : j <= i__3; j += i__2) {
529 work[j] = work[j - *ka];
530 work[*n + j] = work[*n + j - *ka];
531 /* L150: */
532 }
533 i__2 = j1;
534 i__3 = ka1;
535 for (j = j2; i__3 < 0 ? j >= i__2 : j <= i__2; j += i__3) {
536
537 /* create nonzero element a(j-ka,j+1) outside the band */
538 /* and store it in WORK(j) */
539
540 work[j] *= ab[(j + 1) * ab_dim1 + 1];
541 ab[(j + 1) * ab_dim1 + 1] = work[*n + j] * ab[(j + 1) *
542 ab_dim1 + 1];
543 /* L160: */
544 }
545 if (update) {
546 if (i__ - k < *n - *ka && k <= kbt) {
547 work[i__ - k + *ka] = work[i__ - k];
548 }
549 }
550 /* L170: */
551 }
552
553 for (k = *kb; k >= 1; --k) {
554 /* Computing MAX */
555 i__3 = 1, i__2 = k - i0 + 1;
556 j2 = i__ - k - 1 + max(i__3,i__2) * ka1;
557 nr = (*n - j2 + *ka) / ka1;
558 j1 = j2 + (nr - 1) * ka1;
559 if (nr > 0) {
560
561 /* generate rotations in 2nd set to annihilate elements */
562 /* which have been created outside the band */
563
564 slargv_(&nr, &ab[j2 * ab_dim1 + 1], &inca, &work[j2], &ka1, &
565 work[*n + j2], &ka1);
566
567 /* apply rotations in 2nd set from the right */
568
569 i__3 = *ka - 1;
570 for (l = 1; l <= i__3; ++l) {
571 slartv_(&nr, &ab[ka1 - l + j2 * ab_dim1], &inca, &ab[*ka
572 - l + (j2 + 1) * ab_dim1], &inca, &work[*n + j2],
573 &work[j2], &ka1);
574 /* L180: */
575 }
576
577 /* apply rotations in 2nd set from both sides to diagonal */
578 /* blocks */
579
580 slar2v_(&nr, &ab[ka1 + j2 * ab_dim1], &ab[ka1 + (j2 + 1) *
581 ab_dim1], &ab[*ka + (j2 + 1) * ab_dim1], &inca, &work[
582 *n + j2], &work[j2], &ka1);
583
584 }
585
586 /* start applying rotations in 2nd set from the left */
587
588 i__3 = *kb - k + 1;
589 for (l = *ka - 1; l >= i__3; --l) {
590 nrt = (*n - j2 + l) / ka1;
591 if (nrt > 0) {
592 slartv_(&nrt, &ab[l + (j2 + ka1 - l) * ab_dim1], &inca, &
593 ab[l + 1 + (j2 + ka1 - l) * ab_dim1], &inca, &
594 work[*n + j2], &work[j2], &ka1);
595 }
596 /* L190: */
597 }
598
599 if (wantx) {
600
601 /* post-multiply X by product of rotations in 2nd set */
602
603 i__3 = j1;
604 i__2 = ka1;
605 for (j = j2; i__2 < 0 ? j >= i__3 : j <= i__3; j += i__2) {
606 i__4 = *n - m;
607 srot_(&i__4, &x[m + 1 + j * x_dim1], &c__1, &x[m + 1 + (j
608 + 1) * x_dim1], &c__1, &work[*n + j], &work[j]);
609 /* L200: */
610 }
611 }
612 /* L210: */
613 }
614
615 i__2 = *kb - 1;
616 for (k = 1; k <= i__2; ++k) {
617 /* Computing MAX */
618 i__3 = 1, i__4 = k - i0 + 2;
619 j2 = i__ - k - 1 + max(i__3,i__4) * ka1;
620
621 /* finish applying rotations in 1st set from the left */
622
623 for (l = *kb - k; l >= 1; --l) {
624 nrt = (*n - j2 + l) / ka1;
625 if (nrt > 0) {
626 slartv_(&nrt, &ab[l + (j2 + ka1 - l) * ab_dim1], &inca, &
627 ab[l + 1 + (j2 + ka1 - l) * ab_dim1], &inca, &
628 work[*n + j2 - m], &work[j2 - m], &ka1);
629 }
630 /* L220: */
631 }
632 /* L230: */
633 }
634
635 if (*kb > 1) {
636 i__2 = i__ - *kb + (*ka << 1) + 1;
637 for (j = *n - 1; j >= i__2; --j) {
638 work[*n + j - m] = work[*n + j - *ka - m];
639 work[j - m] = work[j - *ka - m];
640 /* L240: */
641 }
642 }
643
644 } else {
645
646 /* Transform A, working with the lower triangle */
647
648 if (update) {
649
650 /* Form inv(S(i))**T * A * inv(S(i)) */
651
652 bii = bb[i__ * bb_dim1 + 1];
653 i__2 = i1;
654 for (j = i__; j <= i__2; ++j) {
655 ab[j - i__ + 1 + i__ * ab_dim1] /= bii;
656 /* L250: */
657 }
658 /* Computing MAX */
659 i__2 = 1, i__3 = i__ - *ka;
660 i__4 = i__;
661 for (j = max(i__2,i__3); j <= i__4; ++j) {
662 ab[i__ - j + 1 + j * ab_dim1] /= bii;
663 /* L260: */
664 }
665 i__4 = i__ - 1;
666 for (k = i__ - kbt; k <= i__4; ++k) {
667 i__2 = k;
668 for (j = i__ - kbt; j <= i__2; ++j) {
669 ab[k - j + 1 + j * ab_dim1] = ab[k - j + 1 + j * ab_dim1]
670 - bb[i__ - j + 1 + j * bb_dim1] * ab[i__ - k + 1
671 + k * ab_dim1] - bb[i__ - k + 1 + k * bb_dim1] *
672 ab[i__ - j + 1 + j * ab_dim1] + ab[i__ * ab_dim1
673 + 1] * bb[i__ - j + 1 + j * bb_dim1] * bb[i__ - k
674 + 1 + k * bb_dim1];
675 /* L270: */
676 }
677 /* Computing MAX */
678 i__2 = 1, i__3 = i__ - *ka;
679 i__1 = i__ - kbt - 1;
680 for (j = max(i__2,i__3); j <= i__1; ++j) {
681 ab[k - j + 1 + j * ab_dim1] -= bb[i__ - k + 1 + k *
682 bb_dim1] * ab[i__ - j + 1 + j * ab_dim1];
683 /* L280: */
684 }
685 /* L290: */
686 }
687 i__4 = i1;
688 for (j = i__; j <= i__4; ++j) {
689 /* Computing MAX */
690 i__1 = j - *ka, i__2 = i__ - kbt;
691 i__3 = i__ - 1;
692 for (k = max(i__1,i__2); k <= i__3; ++k) {
693 ab[j - k + 1 + k * ab_dim1] -= bb[i__ - k + 1 + k *
694 bb_dim1] * ab[j - i__ + 1 + i__ * ab_dim1];
695 /* L300: */
696 }
697 /* L310: */
698 }
699
700 if (wantx) {
701
702 /* post-multiply X by inv(S(i)) */
703
704 i__4 = *n - m;
705 r__1 = 1.f / bii;
706 sscal_(&i__4, &r__1, &x[m + 1 + i__ * x_dim1], &c__1);
707 if (kbt > 0) {
708 i__4 = *n - m;
709 i__3 = *ldbb - 1;
710 sger_(&i__4, &kbt, &c_b20, &x[m + 1 + i__ * x_dim1], &
711 c__1, &bb[kbt + 1 + (i__ - kbt) * bb_dim1], &i__3,
712 &x[m + 1 + (i__ - kbt) * x_dim1], ldx);
713 }
714 }
715
716 /* store a(i1,i) in RA1 for use in next loop over K */
717
718 ra1 = ab[i1 - i__ + 1 + i__ * ab_dim1];
719 }
720
721 /* Generate and apply vectors of rotations to chase all the */
722 /* existing bulges KA positions down toward the bottom of the */
723 /* band */
724
725 i__4 = *kb - 1;
726 for (k = 1; k <= i__4; ++k) {
727 if (update) {
728
729 /* Determine the rotations which would annihilate the bulge */
730 /* which has in theory just been created */
731
732 if (i__ - k + *ka < *n && i__ - k > 1) {
733
734 /* generate rotation to annihilate a(i-k+ka+1,i) */
735
736 slartg_(&ab[ka1 - k + i__ * ab_dim1], &ra1, &work[*n +
737 i__ - k + *ka - m], &work[i__ - k + *ka - m], &ra)
738 ;
739
740 /* create nonzero element a(i-k+ka+1,i-k) outside the */
741 /* band and store it in WORK(i-k) */
742
743 t = -bb[k + 1 + (i__ - k) * bb_dim1] * ra1;
744 work[i__ - k] = work[*n + i__ - k + *ka - m] * t - work[
745 i__ - k + *ka - m] * ab[ka1 + (i__ - k) * ab_dim1]
746 ;
747 ab[ka1 + (i__ - k) * ab_dim1] = work[i__ - k + *ka - m] *
748 t + work[*n + i__ - k + *ka - m] * ab[ka1 + (i__
749 - k) * ab_dim1];
750 ra1 = ra;
751 }
752 }
753 /* Computing MAX */
754 i__3 = 1, i__1 = k - i0 + 2;
755 j2 = i__ - k - 1 + max(i__3,i__1) * ka1;
756 nr = (*n - j2 + *ka) / ka1;
757 j1 = j2 + (nr - 1) * ka1;
758 if (update) {
759 /* Computing MAX */
760 i__3 = j2, i__1 = i__ + (*ka << 1) - k + 1;
761 j2t = max(i__3,i__1);
762 } else {
763 j2t = j2;
764 }
765 nrt = (*n - j2t + *ka) / ka1;
766 i__3 = j1;
767 i__1 = ka1;
768 for (j = j2t; i__1 < 0 ? j >= i__3 : j <= i__3; j += i__1) {
769
770 /* create nonzero element a(j+1,j-ka) outside the band */
771 /* and store it in WORK(j-m) */
772
773 work[j - m] *= ab[ka1 + (j - *ka + 1) * ab_dim1];
774 ab[ka1 + (j - *ka + 1) * ab_dim1] = work[*n + j - m] * ab[ka1
775 + (j - *ka + 1) * ab_dim1];
776 /* L320: */
777 }
778
779 /* generate rotations in 1st set to annihilate elements which */
780 /* have been created outside the band */
781
782 if (nrt > 0) {
783 slargv_(&nrt, &ab[ka1 + (j2t - *ka) * ab_dim1], &inca, &work[
784 j2t - m], &ka1, &work[*n + j2t - m], &ka1);
785 }
786 if (nr > 0) {
787
788 /* apply rotations in 1st set from the left */
789
790 i__1 = *ka - 1;
791 for (l = 1; l <= i__1; ++l) {
792 slartv_(&nr, &ab[l + 1 + (j2 - l) * ab_dim1], &inca, &ab[
793 l + 2 + (j2 - l) * ab_dim1], &inca, &work[*n + j2
794 - m], &work[j2 - m], &ka1);
795 /* L330: */
796 }
797
798 /* apply rotations in 1st set from both sides to diagonal */
799 /* blocks */
800
801 slar2v_(&nr, &ab[j2 * ab_dim1 + 1], &ab[(j2 + 1) * ab_dim1 +
802 1], &ab[j2 * ab_dim1 + 2], &inca, &work[*n + j2 - m],
803 &work[j2 - m], &ka1);
804
805 }
806
807 /* start applying rotations in 1st set from the right */
808
809 i__1 = *kb - k + 1;
810 for (l = *ka - 1; l >= i__1; --l) {
811 nrt = (*n - j2 + l) / ka1;
812 if (nrt > 0) {
813 slartv_(&nrt, &ab[ka1 - l + 1 + j2 * ab_dim1], &inca, &ab[
814 ka1 - l + (j2 + 1) * ab_dim1], &inca, &work[*n +
815 j2 - m], &work[j2 - m], &ka1);
816 }
817 /* L340: */
818 }
819
820 if (wantx) {
821
822 /* post-multiply X by product of rotations in 1st set */
823
824 i__1 = j1;
825 i__3 = ka1;
826 for (j = j2; i__3 < 0 ? j >= i__1 : j <= i__1; j += i__3) {
827 i__2 = *n - m;
828 srot_(&i__2, &x[m + 1 + j * x_dim1], &c__1, &x[m + 1 + (j
829 + 1) * x_dim1], &c__1, &work[*n + j - m], &work[j
830 - m]);
831 /* L350: */
832 }
833 }
834 /* L360: */
835 }
836
837 if (update) {
838 if (i2 <= *n && kbt > 0) {
839
840 /* create nonzero element a(i-kbt+ka+1,i-kbt) outside the */
841 /* band and store it in WORK(i-kbt) */
842
843 work[i__ - kbt] = -bb[kbt + 1 + (i__ - kbt) * bb_dim1] * ra1;
844 }
845 }
846
847 for (k = *kb; k >= 1; --k) {
848 if (update) {
849 /* Computing MAX */
850 i__4 = 2, i__3 = k - i0 + 1;
851 j2 = i__ - k - 1 + max(i__4,i__3) * ka1;
852 } else {
853 /* Computing MAX */
854 i__4 = 1, i__3 = k - i0 + 1;
855 j2 = i__ - k - 1 + max(i__4,i__3) * ka1;
856 }
857
858 /* finish applying rotations in 2nd set from the right */
859
860 for (l = *kb - k; l >= 1; --l) {
861 nrt = (*n - j2 + *ka + l) / ka1;
862 if (nrt > 0) {
863 slartv_(&nrt, &ab[ka1 - l + 1 + (j2 - *ka) * ab_dim1], &
864 inca, &ab[ka1 - l + (j2 - *ka + 1) * ab_dim1], &
865 inca, &work[*n + j2 - *ka], &work[j2 - *ka], &ka1)
866 ;
867 }
868 /* L370: */
869 }
870 nr = (*n - j2 + *ka) / ka1;
871 j1 = j2 + (nr - 1) * ka1;
872 i__4 = j2;
873 i__3 = -ka1;
874 for (j = j1; i__3 < 0 ? j >= i__4 : j <= i__4; j += i__3) {
875 work[j] = work[j - *ka];
876 work[*n + j] = work[*n + j - *ka];
877 /* L380: */
878 }
879 i__3 = j1;
880 i__4 = ka1;
881 for (j = j2; i__4 < 0 ? j >= i__3 : j <= i__3; j += i__4) {
882
883 /* create nonzero element a(j+1,j-ka) outside the band */
884 /* and store it in WORK(j) */
885
886 work[j] *= ab[ka1 + (j - *ka + 1) * ab_dim1];
887 ab[ka1 + (j - *ka + 1) * ab_dim1] = work[*n + j] * ab[ka1 + (
888 j - *ka + 1) * ab_dim1];
889 /* L390: */
890 }
891 if (update) {
892 if (i__ - k < *n - *ka && k <= kbt) {
893 work[i__ - k + *ka] = work[i__ - k];
894 }
895 }
896 /* L400: */
897 }
898
899 for (k = *kb; k >= 1; --k) {
900 /* Computing MAX */
901 i__4 = 1, i__3 = k - i0 + 1;
902 j2 = i__ - k - 1 + max(i__4,i__3) * ka1;
903 nr = (*n - j2 + *ka) / ka1;
904 j1 = j2 + (nr - 1) * ka1;
905 if (nr > 0) {
906
907 /* generate rotations in 2nd set to annihilate elements */
908 /* which have been created outside the band */
909
910 slargv_(&nr, &ab[ka1 + (j2 - *ka) * ab_dim1], &inca, &work[j2]
911 , &ka1, &work[*n + j2], &ka1);
912
913 /* apply rotations in 2nd set from the left */
914
915 i__4 = *ka - 1;
916 for (l = 1; l <= i__4; ++l) {
917 slartv_(&nr, &ab[l + 1 + (j2 - l) * ab_dim1], &inca, &ab[
918 l + 2 + (j2 - l) * ab_dim1], &inca, &work[*n + j2]
919 , &work[j2], &ka1);
920 /* L410: */
921 }
922
923 /* apply rotations in 2nd set from both sides to diagonal */
924 /* blocks */
925
926 slar2v_(&nr, &ab[j2 * ab_dim1 + 1], &ab[(j2 + 1) * ab_dim1 +
927 1], &ab[j2 * ab_dim1 + 2], &inca, &work[*n + j2], &
928 work[j2], &ka1);
929
930 }
931
932 /* start applying rotations in 2nd set from the right */
933
934 i__4 = *kb - k + 1;
935 for (l = *ka - 1; l >= i__4; --l) {
936 nrt = (*n - j2 + l) / ka1;
937 if (nrt > 0) {
938 slartv_(&nrt, &ab[ka1 - l + 1 + j2 * ab_dim1], &inca, &ab[
939 ka1 - l + (j2 + 1) * ab_dim1], &inca, &work[*n +
940 j2], &work[j2], &ka1);
941 }
942 /* L420: */
943 }
944
945 if (wantx) {
946
947 /* post-multiply X by product of rotations in 2nd set */
948
949 i__4 = j1;
950 i__3 = ka1;
951 for (j = j2; i__3 < 0 ? j >= i__4 : j <= i__4; j += i__3) {
952 i__1 = *n - m;
953 srot_(&i__1, &x[m + 1 + j * x_dim1], &c__1, &x[m + 1 + (j
954 + 1) * x_dim1], &c__1, &work[*n + j], &work[j]);
955 /* L430: */
956 }
957 }
958 /* L440: */
959 }
960
961 i__3 = *kb - 1;
962 for (k = 1; k <= i__3; ++k) {
963 /* Computing MAX */
964 i__4 = 1, i__1 = k - i0 + 2;
965 j2 = i__ - k - 1 + max(i__4,i__1) * ka1;
966
967 /* finish applying rotations in 1st set from the right */
968
969 for (l = *kb - k; l >= 1; --l) {
970 nrt = (*n - j2 + l) / ka1;
971 if (nrt > 0) {
972 slartv_(&nrt, &ab[ka1 - l + 1 + j2 * ab_dim1], &inca, &ab[
973 ka1 - l + (j2 + 1) * ab_dim1], &inca, &work[*n +
974 j2 - m], &work[j2 - m], &ka1);
975 }
976 /* L450: */
977 }
978 /* L460: */
979 }
980
981 if (*kb > 1) {
982 i__3 = i__ - *kb + (*ka << 1) + 1;
983 for (j = *n - 1; j >= i__3; --j) {
984 work[*n + j - m] = work[*n + j - *ka - m];
985 work[j - m] = work[j - *ka - m];
986 /* L470: */
987 }
988 }
989
990 }
991
992 goto L10;
993
994 L480:
995
996 /* **************************** Phase 2 ***************************** */
997
998 /* The logical structure of this phase is: */
999
1000 /* UPDATE = .TRUE. */
1001 /* DO I = 1, M */
1002 /* use S(i) to update A and create a new bulge */
1003 /* apply rotations to push all bulges KA positions upward */
1004 /* END DO */
1005 /* UPDATE = .FALSE. */
1006 /* DO I = M - KA - 1, 2, -1 */
1007 /* apply rotations to push all bulges KA positions upward */
1008 /* END DO */
1009
1010 /* To avoid duplicating code, the two loops are merged. */
1011
1012 update = TRUE_;
1013 i__ = 0;
1014 L490:
1015 if (update) {
1016 ++i__;
1017 /* Computing MIN */
1018 i__3 = *kb, i__4 = m - i__;
1019 kbt = min(i__3,i__4);
1020 i0 = i__ + 1;
1021 /* Computing MAX */
1022 i__3 = 1, i__4 = i__ - *ka;
1023 i1 = max(i__3,i__4);
1024 i2 = i__ + kbt - ka1;
1025 if (i__ > m) {
1026 update = FALSE_;
1027 --i__;
1028 i0 = m + 1;
1029 if (*ka == 0) {
1030 return 0;
1031 }
1032 goto L490;
1033 }
1034 } else {
1035 i__ -= *ka;
1036 if (i__ < 2) {
1037 return 0;
1038 }
1039 }
1040
1041 if (i__ < m - kbt) {
1042 nx = m;
1043 } else {
1044 nx = *n;
1045 }
1046
1047 if (upper) {
1048
1049 /* Transform A, working with the upper triangle */
1050
1051 if (update) {
1052
1053 /* Form inv(S(i))**T * A * inv(S(i)) */
1054
1055 bii = bb[kb1 + i__ * bb_dim1];
1056 i__3 = i__;
1057 for (j = i1; j <= i__3; ++j) {
1058 ab[j - i__ + ka1 + i__ * ab_dim1] /= bii;
1059 /* L500: */
1060 }
1061 /* Computing MIN */
1062 i__4 = *n, i__1 = i__ + *ka;
1063 i__3 = min(i__4,i__1);
1064 for (j = i__; j <= i__3; ++j) {
1065 ab[i__ - j + ka1 + j * ab_dim1] /= bii;
1066 /* L510: */
1067 }
1068 i__3 = i__ + kbt;
1069 for (k = i__ + 1; k <= i__3; ++k) {
1070 i__4 = i__ + kbt;
1071 for (j = k; j <= i__4; ++j) {
1072 ab[k - j + ka1 + j * ab_dim1] = ab[k - j + ka1 + j *
1073 ab_dim1] - bb[i__ - j + kb1 + j * bb_dim1] * ab[
1074 i__ - k + ka1 + k * ab_dim1] - bb[i__ - k + kb1 +
1075 k * bb_dim1] * ab[i__ - j + ka1 + j * ab_dim1] +
1076 ab[ka1 + i__ * ab_dim1] * bb[i__ - j + kb1 + j *
1077 bb_dim1] * bb[i__ - k + kb1 + k * bb_dim1];
1078 /* L520: */
1079 }
1080 /* Computing MIN */
1081 i__1 = *n, i__2 = i__ + *ka;
1082 i__4 = min(i__1,i__2);
1083 for (j = i__ + kbt + 1; j <= i__4; ++j) {
1084 ab[k - j + ka1 + j * ab_dim1] -= bb[i__ - k + kb1 + k *
1085 bb_dim1] * ab[i__ - j + ka1 + j * ab_dim1];
1086 /* L530: */
1087 }
1088 /* L540: */
1089 }
1090 i__3 = i__;
1091 for (j = i1; j <= i__3; ++j) {
1092 /* Computing MIN */
1093 i__1 = j + *ka, i__2 = i__ + kbt;
1094 i__4 = min(i__1,i__2);
1095 for (k = i__ + 1; k <= i__4; ++k) {
1096 ab[j - k + ka1 + k * ab_dim1] -= bb[i__ - k + kb1 + k *
1097 bb_dim1] * ab[j - i__ + ka1 + i__ * ab_dim1];
1098 /* L550: */
1099 }
1100 /* L560: */
1101 }
1102
1103 if (wantx) {
1104
1105 /* post-multiply X by inv(S(i)) */
1106
1107 r__1 = 1.f / bii;
1108 sscal_(&nx, &r__1, &x[i__ * x_dim1 + 1], &c__1);
1109 if (kbt > 0) {
1110 i__3 = *ldbb - 1;
1111 sger_(&nx, &kbt, &c_b20, &x[i__ * x_dim1 + 1], &c__1, &bb[
1112 *kb + (i__ + 1) * bb_dim1], &i__3, &x[(i__ + 1) *
1113 x_dim1 + 1], ldx);
1114 }
1115 }
1116
1117 /* store a(i1,i) in RA1 for use in next loop over K */
1118
1119 ra1 = ab[i1 - i__ + ka1 + i__ * ab_dim1];
1120 }
1121
1122 /* Generate and apply vectors of rotations to chase all the */
1123 /* existing bulges KA positions up toward the top of the band */
1124
1125 i__3 = *kb - 1;
1126 for (k = 1; k <= i__3; ++k) {
1127 if (update) {
1128
1129 /* Determine the rotations which would annihilate the bulge */
1130 /* which has in theory just been created */
1131
1132 if (i__ + k - ka1 > 0 && i__ + k < m) {
1133
1134 /* generate rotation to annihilate a(i+k-ka-1,i) */
1135
1136 slartg_(&ab[k + 1 + i__ * ab_dim1], &ra1, &work[*n + i__
1137 + k - *ka], &work[i__ + k - *ka], &ra);
1138
1139 /* create nonzero element a(i+k-ka-1,i+k) outside the */
1140 /* band and store it in WORK(m-kb+i+k) */
1141
1142 t = -bb[kb1 - k + (i__ + k) * bb_dim1] * ra1;
1143 work[m - *kb + i__ + k] = work[*n + i__ + k - *ka] * t -
1144 work[i__ + k - *ka] * ab[(i__ + k) * ab_dim1 + 1];
1145 ab[(i__ + k) * ab_dim1 + 1] = work[i__ + k - *ka] * t +
1146 work[*n + i__ + k - *ka] * ab[(i__ + k) * ab_dim1
1147 + 1];
1148 ra1 = ra;
1149 }
1150 }
1151 /* Computing MAX */
1152 i__4 = 1, i__1 = k + i0 - m + 1;
1153 j2 = i__ + k + 1 - max(i__4,i__1) * ka1;
1154 nr = (j2 + *ka - 1) / ka1;
1155 j1 = j2 - (nr - 1) * ka1;
1156 if (update) {
1157 /* Computing MIN */
1158 i__4 = j2, i__1 = i__ - (*ka << 1) + k - 1;
1159 j2t = min(i__4,i__1);
1160 } else {
1161 j2t = j2;
1162 }
1163 nrt = (j2t + *ka - 1) / ka1;
1164 i__4 = j2t;
1165 i__1 = ka1;
1166 for (j = j1; i__1 < 0 ? j >= i__4 : j <= i__4; j += i__1) {
1167
1168 /* create nonzero element a(j-1,j+ka) outside the band */
1169 /* and store it in WORK(j) */
1170
1171 work[j] *= ab[(j + *ka - 1) * ab_dim1 + 1];
1172 ab[(j + *ka - 1) * ab_dim1 + 1] = work[*n + j] * ab[(j + *ka
1173 - 1) * ab_dim1 + 1];
1174 /* L570: */
1175 }
1176
1177 /* generate rotations in 1st set to annihilate elements which */
1178 /* have been created outside the band */
1179
1180 if (nrt > 0) {
1181 slargv_(&nrt, &ab[(j1 + *ka) * ab_dim1 + 1], &inca, &work[j1],
1182 &ka1, &work[*n + j1], &ka1);
1183 }
1184 if (nr > 0) {
1185
1186 /* apply rotations in 1st set from the left */
1187
1188 i__1 = *ka - 1;
1189 for (l = 1; l <= i__1; ++l) {
1190 slartv_(&nr, &ab[ka1 - l + (j1 + l) * ab_dim1], &inca, &
1191 ab[*ka - l + (j1 + l) * ab_dim1], &inca, &work[*n
1192 + j1], &work[j1], &ka1);
1193 /* L580: */
1194 }
1195
1196 /* apply rotations in 1st set from both sides to diagonal */
1197 /* blocks */
1198
1199 slar2v_(&nr, &ab[ka1 + j1 * ab_dim1], &ab[ka1 + (j1 - 1) *
1200 ab_dim1], &ab[*ka + j1 * ab_dim1], &inca, &work[*n +
1201 j1], &work[j1], &ka1);
1202
1203 }
1204
1205 /* start applying rotations in 1st set from the right */
1206
1207 i__1 = *kb - k + 1;
1208 for (l = *ka - 1; l >= i__1; --l) {
1209 nrt = (j2 + l - 1) / ka1;
1210 j1t = j2 - (nrt - 1) * ka1;
1211 if (nrt > 0) {
1212 slartv_(&nrt, &ab[l + j1t * ab_dim1], &inca, &ab[l + 1 + (
1213 j1t - 1) * ab_dim1], &inca, &work[*n + j1t], &
1214 work[j1t], &ka1);
1215 }
1216 /* L590: */
1217 }
1218
1219 if (wantx) {
1220
1221 /* post-multiply X by product of rotations in 1st set */
1222
1223 i__1 = j2;
1224 i__4 = ka1;
1225 for (j = j1; i__4 < 0 ? j >= i__1 : j <= i__1; j += i__4) {
1226 srot_(&nx, &x[j * x_dim1 + 1], &c__1, &x[(j - 1) * x_dim1
1227 + 1], &c__1, &work[*n + j], &work[j]);
1228 /* L600: */
1229 }
1230 }
1231 /* L610: */
1232 }
1233
1234 if (update) {
1235 if (i2 > 0 && kbt > 0) {
1236
1237 /* create nonzero element a(i+kbt-ka-1,i+kbt) outside the */
1238 /* band and store it in WORK(m-kb+i+kbt) */
1239
1240 work[m - *kb + i__ + kbt] = -bb[kb1 - kbt + (i__ + kbt) *
1241 bb_dim1] * ra1;
1242 }
1243 }
1244
1245 for (k = *kb; k >= 1; --k) {
1246 if (update) {
1247 /* Computing MAX */
1248 i__3 = 2, i__4 = k + i0 - m;
1249 j2 = i__ + k + 1 - max(i__3,i__4) * ka1;
1250 } else {
1251 /* Computing MAX */
1252 i__3 = 1, i__4 = k + i0 - m;
1253 j2 = i__ + k + 1 - max(i__3,i__4) * ka1;
1254 }
1255
1256 /* finish applying rotations in 2nd set from the right */
1257
1258 for (l = *kb - k; l >= 1; --l) {
1259 nrt = (j2 + *ka + l - 1) / ka1;
1260 j1t = j2 - (nrt - 1) * ka1;
1261 if (nrt > 0) {
1262 slartv_(&nrt, &ab[l + (j1t + *ka) * ab_dim1], &inca, &ab[
1263 l + 1 + (j1t + *ka - 1) * ab_dim1], &inca, &work[*
1264 n + m - *kb + j1t + *ka], &work[m - *kb + j1t + *
1265 ka], &ka1);
1266 }
1267 /* L620: */
1268 }
1269 nr = (j2 + *ka - 1) / ka1;
1270 j1 = j2 - (nr - 1) * ka1;
1271 i__3 = j2;
1272 i__4 = ka1;
1273 for (j = j1; i__4 < 0 ? j >= i__3 : j <= i__3; j += i__4) {
1274 work[m - *kb + j] = work[m - *kb + j + *ka];
1275 work[*n + m - *kb + j] = work[*n + m - *kb + j + *ka];
1276 /* L630: */
1277 }
1278 i__4 = j2;
1279 i__3 = ka1;
1280 for (j = j1; i__3 < 0 ? j >= i__4 : j <= i__4; j += i__3) {
1281
1282 /* create nonzero element a(j-1,j+ka) outside the band */
1283 /* and store it in WORK(m-kb+j) */
1284
1285 work[m - *kb + j] *= ab[(j + *ka - 1) * ab_dim1 + 1];
1286 ab[(j + *ka - 1) * ab_dim1 + 1] = work[*n + m - *kb + j] * ab[
1287 (j + *ka - 1) * ab_dim1 + 1];
1288 /* L640: */
1289 }
1290 if (update) {
1291 if (i__ + k > ka1 && k <= kbt) {
1292 work[m - *kb + i__ + k - *ka] = work[m - *kb + i__ + k];
1293 }
1294 }
1295 /* L650: */
1296 }
1297
1298 for (k = *kb; k >= 1; --k) {
1299 /* Computing MAX */
1300 i__3 = 1, i__4 = k + i0 - m;
1301 j2 = i__ + k + 1 - max(i__3,i__4) * ka1;
1302 nr = (j2 + *ka - 1) / ka1;
1303 j1 = j2 - (nr - 1) * ka1;
1304 if (nr > 0) {
1305
1306 /* generate rotations in 2nd set to annihilate elements */
1307 /* which have been created outside the band */
1308
1309 slargv_(&nr, &ab[(j1 + *ka) * ab_dim1 + 1], &inca, &work[m - *
1310 kb + j1], &ka1, &work[*n + m - *kb + j1], &ka1);
1311
1312 /* apply rotations in 2nd set from the left */
1313
1314 i__3 = *ka - 1;
1315 for (l = 1; l <= i__3; ++l) {
1316 slartv_(&nr, &ab[ka1 - l + (j1 + l) * ab_dim1], &inca, &
1317 ab[*ka - l + (j1 + l) * ab_dim1], &inca, &work[*n
1318 + m - *kb + j1], &work[m - *kb + j1], &ka1);
1319 /* L660: */
1320 }
1321
1322 /* apply rotations in 2nd set from both sides to diagonal */
1323 /* blocks */
1324
1325 slar2v_(&nr, &ab[ka1 + j1 * ab_dim1], &ab[ka1 + (j1 - 1) *
1326 ab_dim1], &ab[*ka + j1 * ab_dim1], &inca, &work[*n +
1327 m - *kb + j1], &work[m - *kb + j1], &ka1);
1328
1329 }
1330
1331 /* start applying rotations in 2nd set from the right */
1332
1333 i__3 = *kb - k + 1;
1334 for (l = *ka - 1; l >= i__3; --l) {
1335 nrt = (j2 + l - 1) / ka1;
1336 j1t = j2 - (nrt - 1) * ka1;
1337 if (nrt > 0) {
1338 slartv_(&nrt, &ab[l + j1t * ab_dim1], &inca, &ab[l + 1 + (
1339 j1t - 1) * ab_dim1], &inca, &work[*n + m - *kb +
1340 j1t], &work[m - *kb + j1t], &ka1);
1341 }
1342 /* L670: */
1343 }
1344
1345 if (wantx) {
1346
1347 /* post-multiply X by product of rotations in 2nd set */
1348
1349 i__3 = j2;
1350 i__4 = ka1;
1351 for (j = j1; i__4 < 0 ? j >= i__3 : j <= i__3; j += i__4) {
1352 srot_(&nx, &x[j * x_dim1 + 1], &c__1, &x[(j - 1) * x_dim1
1353 + 1], &c__1, &work[*n + m - *kb + j], &work[m - *
1354 kb + j]);
1355 /* L680: */
1356 }
1357 }
1358 /* L690: */
1359 }
1360
1361 i__4 = *kb - 1;
1362 for (k = 1; k <= i__4; ++k) {
1363 /* Computing MAX */
1364 i__3 = 1, i__1 = k + i0 - m + 1;
1365 j2 = i__ + k + 1 - max(i__3,i__1) * ka1;
1366
1367 /* finish applying rotations in 1st set from the right */
1368
1369 for (l = *kb - k; l >= 1; --l) {
1370 nrt = (j2 + l - 1) / ka1;
1371 j1t = j2 - (nrt - 1) * ka1;
1372 if (nrt > 0) {
1373 slartv_(&nrt, &ab[l + j1t * ab_dim1], &inca, &ab[l + 1 + (
1374 j1t - 1) * ab_dim1], &inca, &work[*n + j1t], &
1375 work[j1t], &ka1);
1376 }
1377 /* L700: */
1378 }
1379 /* L710: */
1380 }
1381
1382 if (*kb > 1) {
1383 /* Computing MIN */
1384 i__3 = i__ + *kb;
1385 i__4 = min(i__3,m) - (*ka << 1) - 1;
1386 for (j = 2; j <= i__4; ++j) {
1387 work[*n + j] = work[*n + j + *ka];
1388 work[j] = work[j + *ka];
1389 /* L720: */
1390 }
1391 }
1392
1393 } else {
1394
1395 /* Transform A, working with the lower triangle */
1396
1397 if (update) {
1398
1399 /* Form inv(S(i))**T * A * inv(S(i)) */
1400
1401 bii = bb[i__ * bb_dim1 + 1];
1402 i__4 = i__;
1403 for (j = i1; j <= i__4; ++j) {
1404 ab[i__ - j + 1 + j * ab_dim1] /= bii;
1405 /* L730: */
1406 }
1407 /* Computing MIN */
1408 i__3 = *n, i__1 = i__ + *ka;
1409 i__4 = min(i__3,i__1);
1410 for (j = i__; j <= i__4; ++j) {
1411 ab[j - i__ + 1 + i__ * ab_dim1] /= bii;
1412 /* L740: */
1413 }
1414 i__4 = i__ + kbt;
1415 for (k = i__ + 1; k <= i__4; ++k) {
1416 i__3 = i__ + kbt;
1417 for (j = k; j <= i__3; ++j) {
1418 ab[j - k + 1 + k * ab_dim1] = ab[j - k + 1 + k * ab_dim1]
1419 - bb[j - i__ + 1 + i__ * bb_dim1] * ab[k - i__ +
1420 1 + i__ * ab_dim1] - bb[k - i__ + 1 + i__ *
1421 bb_dim1] * ab[j - i__ + 1 + i__ * ab_dim1] + ab[
1422 i__ * ab_dim1 + 1] * bb[j - i__ + 1 + i__ *
1423 bb_dim1] * bb[k - i__ + 1 + i__ * bb_dim1];
1424 /* L750: */
1425 }
1426 /* Computing MIN */
1427 i__1 = *n, i__2 = i__ + *ka;
1428 i__3 = min(i__1,i__2);
1429 for (j = i__ + kbt + 1; j <= i__3; ++j) {
1430 ab[j - k + 1 + k * ab_dim1] -= bb[k - i__ + 1 + i__ *
1431 bb_dim1] * ab[j - i__ + 1 + i__ * ab_dim1];
1432 /* L760: */
1433 }
1434 /* L770: */
1435 }
1436 i__4 = i__;
1437 for (j = i1; j <= i__4; ++j) {
1438 /* Computing MIN */
1439 i__1 = j + *ka, i__2 = i__ + kbt;
1440 i__3 = min(i__1,i__2);
1441 for (k = i__ + 1; k <= i__3; ++k) {
1442 ab[k - j + 1 + j * ab_dim1] -= bb[k - i__ + 1 + i__ *
1443 bb_dim1] * ab[i__ - j + 1 + j * ab_dim1];
1444 /* L780: */
1445 }
1446 /* L790: */
1447 }
1448
1449 if (wantx) {
1450
1451 /* post-multiply X by inv(S(i)) */
1452
1453 r__1 = 1.f / bii;
1454 sscal_(&nx, &r__1, &x[i__ * x_dim1 + 1], &c__1);
1455 if (kbt > 0) {
1456 sger_(&nx, &kbt, &c_b20, &x[i__ * x_dim1 + 1], &c__1, &bb[
1457 i__ * bb_dim1 + 2], &c__1, &x[(i__ + 1) * x_dim1
1458 + 1], ldx);
1459 }
1460 }
1461
1462 /* store a(i,i1) in RA1 for use in next loop over K */
1463
1464 ra1 = ab[i__ - i1 + 1 + i1 * ab_dim1];
1465 }
1466
1467 /* Generate and apply vectors of rotations to chase all the */
1468 /* existing bulges KA positions up toward the top of the band */
1469
1470 i__4 = *kb - 1;
1471 for (k = 1; k <= i__4; ++k) {
1472 if (update) {
1473
1474 /* Determine the rotations which would annihilate the bulge */
1475 /* which has in theory just been created */
1476
1477 if (i__ + k - ka1 > 0 && i__ + k < m) {
1478
1479 /* generate rotation to annihilate a(i,i+k-ka-1) */
1480
1481 slartg_(&ab[ka1 - k + (i__ + k - *ka) * ab_dim1], &ra1, &
1482 work[*n + i__ + k - *ka], &work[i__ + k - *ka], &
1483 ra);
1484
1485 /* create nonzero element a(i+k,i+k-ka-1) outside the */
1486 /* band and store it in WORK(m-kb+i+k) */
1487
1488 t = -bb[k + 1 + i__ * bb_dim1] * ra1;
1489 work[m - *kb + i__ + k] = work[*n + i__ + k - *ka] * t -
1490 work[i__ + k - *ka] * ab[ka1 + (i__ + k - *ka) *
1491 ab_dim1];
1492 ab[ka1 + (i__ + k - *ka) * ab_dim1] = work[i__ + k - *ka]
1493 * t + work[*n + i__ + k - *ka] * ab[ka1 + (i__ +
1494 k - *ka) * ab_dim1];
1495 ra1 = ra;
1496 }
1497 }
1498 /* Computing MAX */
1499 i__3 = 1, i__1 = k + i0 - m + 1;
1500 j2 = i__ + k + 1 - max(i__3,i__1) * ka1;
1501 nr = (j2 + *ka - 1) / ka1;
1502 j1 = j2 - (nr - 1) * ka1;
1503 if (update) {
1504 /* Computing MIN */
1505 i__3 = j2, i__1 = i__ - (*ka << 1) + k - 1;
1506 j2t = min(i__3,i__1);
1507 } else {
1508 j2t = j2;
1509 }
1510 nrt = (j2t + *ka - 1) / ka1;
1511 i__3 = j2t;
1512 i__1 = ka1;
1513 for (j = j1; i__1 < 0 ? j >= i__3 : j <= i__3; j += i__1) {
1514
1515 /* create nonzero element a(j+ka,j-1) outside the band */
1516 /* and store it in WORK(j) */
1517
1518 work[j] *= ab[ka1 + (j - 1) * ab_dim1];
1519 ab[ka1 + (j - 1) * ab_dim1] = work[*n + j] * ab[ka1 + (j - 1)
1520 * ab_dim1];
1521 /* L800: */
1522 }
1523
1524 /* generate rotations in 1st set to annihilate elements which */
1525 /* have been created outside the band */
1526
1527 if (nrt > 0) {
1528 slargv_(&nrt, &ab[ka1 + j1 * ab_dim1], &inca, &work[j1], &ka1,
1529 &work[*n + j1], &ka1);
1530 }
1531 if (nr > 0) {
1532
1533 /* apply rotations in 1st set from the right */
1534
1535 i__1 = *ka - 1;
1536 for (l = 1; l <= i__1; ++l) {
1537 slartv_(&nr, &ab[l + 1 + j1 * ab_dim1], &inca, &ab[l + 2
1538 + (j1 - 1) * ab_dim1], &inca, &work[*n + j1], &
1539 work[j1], &ka1);
1540 /* L810: */
1541 }
1542
1543 /* apply rotations in 1st set from both sides to diagonal */
1544 /* blocks */
1545
1546 slar2v_(&nr, &ab[j1 * ab_dim1 + 1], &ab[(j1 - 1) * ab_dim1 +
1547 1], &ab[(j1 - 1) * ab_dim1 + 2], &inca, &work[*n + j1]
1548 , &work[j1], &ka1);
1549
1550 }
1551
1552 /* start applying rotations in 1st set from the left */
1553
1554 i__1 = *kb - k + 1;
1555 for (l = *ka - 1; l >= i__1; --l) {
1556 nrt = (j2 + l - 1) / ka1;
1557 j1t = j2 - (nrt - 1) * ka1;
1558 if (nrt > 0) {
1559 slartv_(&nrt, &ab[ka1 - l + 1 + (j1t - ka1 + l) * ab_dim1]
1560 , &inca, &ab[ka1 - l + (j1t - ka1 + l) * ab_dim1],
1561 &inca, &work[*n + j1t], &work[j1t], &ka1);
1562 }
1563 /* L820: */
1564 }
1565
1566 if (wantx) {
1567
1568 /* post-multiply X by product of rotations in 1st set */
1569
1570 i__1 = j2;
1571 i__3 = ka1;
1572 for (j = j1; i__3 < 0 ? j >= i__1 : j <= i__1; j += i__3) {
1573 srot_(&nx, &x[j * x_dim1 + 1], &c__1, &x[(j - 1) * x_dim1
1574 + 1], &c__1, &work[*n + j], &work[j]);
1575 /* L830: */
1576 }
1577 }
1578 /* L840: */
1579 }
1580
1581 if (update) {
1582 if (i2 > 0 && kbt > 0) {
1583
1584 /* create nonzero element a(i+kbt,i+kbt-ka-1) outside the */
1585 /* band and store it in WORK(m-kb+i+kbt) */
1586
1587 work[m - *kb + i__ + kbt] = -bb[kbt + 1 + i__ * bb_dim1] *
1588 ra1;
1589 }
1590 }
1591
1592 for (k = *kb; k >= 1; --k) {
1593 if (update) {
1594 /* Computing MAX */
1595 i__4 = 2, i__3 = k + i0 - m;
1596 j2 = i__ + k + 1 - max(i__4,i__3) * ka1;
1597 } else {
1598 /* Computing MAX */
1599 i__4 = 1, i__3 = k + i0 - m;
1600 j2 = i__ + k + 1 - max(i__4,i__3) * ka1;
1601 }
1602
1603 /* finish applying rotations in 2nd set from the left */
1604
1605 for (l = *kb - k; l >= 1; --l) {
1606 nrt = (j2 + *ka + l - 1) / ka1;
1607 j1t = j2 - (nrt - 1) * ka1;
1608 if (nrt > 0) {
1609 slartv_(&nrt, &ab[ka1 - l + 1 + (j1t + l - 1) * ab_dim1],
1610 &inca, &ab[ka1 - l + (j1t + l - 1) * ab_dim1], &
1611 inca, &work[*n + m - *kb + j1t + *ka], &work[m - *
1612 kb + j1t + *ka], &ka1);
1613 }
1614 /* L850: */
1615 }
1616 nr = (j2 + *ka - 1) / ka1;
1617 j1 = j2 - (nr - 1) * ka1;
1618 i__4 = j2;
1619 i__3 = ka1;
1620 for (j = j1; i__3 < 0 ? j >= i__4 : j <= i__4; j += i__3) {
1621 work[m - *kb + j] = work[m - *kb + j + *ka];
1622 work[*n + m - *kb + j] = work[*n + m - *kb + j + *ka];
1623 /* L860: */
1624 }
1625 i__3 = j2;
1626 i__4 = ka1;
1627 for (j = j1; i__4 < 0 ? j >= i__3 : j <= i__3; j += i__4) {
1628
1629 /* create nonzero element a(j+ka,j-1) outside the band */
1630 /* and store it in WORK(m-kb+j) */
1631
1632 work[m - *kb + j] *= ab[ka1 + (j - 1) * ab_dim1];
1633 ab[ka1 + (j - 1) * ab_dim1] = work[*n + m - *kb + j] * ab[ka1
1634 + (j - 1) * ab_dim1];
1635 /* L870: */
1636 }
1637 if (update) {
1638 if (i__ + k > ka1 && k <= kbt) {
1639 work[m - *kb + i__ + k - *ka] = work[m - *kb + i__ + k];
1640 }
1641 }
1642 /* L880: */
1643 }
1644
1645 for (k = *kb; k >= 1; --k) {
1646 /* Computing MAX */
1647 i__4 = 1, i__3 = k + i0 - m;
1648 j2 = i__ + k + 1 - max(i__4,i__3) * ka1;
1649 nr = (j2 + *ka - 1) / ka1;
1650 j1 = j2 - (nr - 1) * ka1;
1651 if (nr > 0) {
1652
1653 /* generate rotations in 2nd set to annihilate elements */
1654 /* which have been created outside the band */
1655
1656 slargv_(&nr, &ab[ka1 + j1 * ab_dim1], &inca, &work[m - *kb +
1657 j1], &ka1, &work[*n + m - *kb + j1], &ka1);
1658
1659 /* apply rotations in 2nd set from the right */
1660
1661 i__4 = *ka - 1;
1662 for (l = 1; l <= i__4; ++l) {
1663 slartv_(&nr, &ab[l + 1 + j1 * ab_dim1], &inca, &ab[l + 2
1664 + (j1 - 1) * ab_dim1], &inca, &work[*n + m - *kb
1665 + j1], &work[m - *kb + j1], &ka1);
1666 /* L890: */
1667 }
1668
1669 /* apply rotations in 2nd set from both sides to diagonal */
1670 /* blocks */
1671
1672 slar2v_(&nr, &ab[j1 * ab_dim1 + 1], &ab[(j1 - 1) * ab_dim1 +
1673 1], &ab[(j1 - 1) * ab_dim1 + 2], &inca, &work[*n + m
1674 - *kb + j1], &work[m - *kb + j1], &ka1);
1675
1676 }
1677
1678 /* start applying rotations in 2nd set from the left */
1679
1680 i__4 = *kb - k + 1;
1681 for (l = *ka - 1; l >= i__4; --l) {
1682 nrt = (j2 + l - 1) / ka1;
1683 j1t = j2 - (nrt - 1) * ka1;
1684 if (nrt > 0) {
1685 slartv_(&nrt, &ab[ka1 - l + 1 + (j1t - ka1 + l) * ab_dim1]
1686 , &inca, &ab[ka1 - l + (j1t - ka1 + l) * ab_dim1],
1687 &inca, &work[*n + m - *kb + j1t], &work[m - *kb
1688 + j1t], &ka1);
1689 }
1690 /* L900: */
1691 }
1692
1693 if (wantx) {
1694
1695 /* post-multiply X by product of rotations in 2nd set */
1696
1697 i__4 = j2;
1698 i__3 = ka1;
1699 for (j = j1; i__3 < 0 ? j >= i__4 : j <= i__4; j += i__3) {
1700 srot_(&nx, &x[j * x_dim1 + 1], &c__1, &x[(j - 1) * x_dim1
1701 + 1], &c__1, &work[*n + m - *kb + j], &work[m - *
1702 kb + j]);
1703 /* L910: */
1704 }
1705 }
1706 /* L920: */
1707 }
1708
1709 i__3 = *kb - 1;
1710 for (k = 1; k <= i__3; ++k) {
1711 /* Computing MAX */
1712 i__4 = 1, i__1 = k + i0 - m + 1;
1713 j2 = i__ + k + 1 - max(i__4,i__1) * ka1;
1714
1715 /* finish applying rotations in 1st set from the left */
1716
1717 for (l = *kb - k; l >= 1; --l) {
1718 nrt = (j2 + l - 1) / ka1;
1719 j1t = j2 - (nrt - 1) * ka1;
1720 if (nrt > 0) {
1721 slartv_(&nrt, &ab[ka1 - l + 1 + (j1t - ka1 + l) * ab_dim1]
1722 , &inca, &ab[ka1 - l + (j1t - ka1 + l) * ab_dim1],
1723 &inca, &work[*n + j1t], &work[j1t], &ka1);
1724 }
1725 /* L930: */
1726 }
1727 /* L940: */
1728 }
1729
1730 if (*kb > 1) {
1731 /* Computing MIN */
1732 i__4 = i__ + *kb;
1733 i__3 = min(i__4,m) - (*ka << 1) - 1;
1734 for (j = 2; j <= i__3; ++j) {
1735 work[*n + j] = work[*n + j + *ka];
1736 work[j] = work[j + *ka];
1737 /* L950: */
1738 }
1739 }
1740
1741 }
1742
1743 goto L490;
1744
1745 /* End of SSBGST */
1746
1747 } /* ssbgst_ */
1748
1749