1 SUBROUTINE DSBGST( VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X,
2 $ LDX, WORK, INFO )
3 *
4 * -- LAPACK routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER UPLO, VECT
11 INTEGER INFO, KA, KB, LDAB, LDBB, LDX, N
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
15 $ X( LDX, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DSBGST reduces a real symmetric-definite banded generalized
22 * eigenproblem A*x = lambda*B*x to standard form C*y = lambda*y,
23 * such that C has the same bandwidth as A.
24 *
25 * B must have been previously factorized as S**T*S by DPBSTF, using a
26 * split Cholesky factorization. A is overwritten by C = X**T*A*X, where
27 * X = S**(-1)*Q and Q is an orthogonal matrix chosen to preserve the
28 * bandwidth of A.
29 *
30 * Arguments
31 * =========
32 *
33 * VECT (input) CHARACTER*1
34 * = 'N': do not form the transformation matrix X;
35 * = 'V': form X.
36 *
37 * UPLO (input) CHARACTER*1
38 * = 'U': Upper triangle of A is stored;
39 * = 'L': Lower triangle of A is stored.
40 *
41 * N (input) INTEGER
42 * The order of the matrices A and B. N >= 0.
43 *
44 * KA (input) INTEGER
45 * The number of superdiagonals of the matrix A if UPLO = 'U',
46 * or the number of subdiagonals if UPLO = 'L'. KA >= 0.
47 *
48 * KB (input) INTEGER
49 * The number of superdiagonals of the matrix B if UPLO = 'U',
50 * or the number of subdiagonals if UPLO = 'L'. KA >= KB >= 0.
51 *
52 * AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
53 * On entry, the upper or lower triangle of the symmetric band
54 * matrix A, stored in the first ka+1 rows of the array. The
55 * j-th column of A is stored in the j-th column of the array AB
56 * as follows:
57 * if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
58 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+ka).
59 *
60 * On exit, the transformed matrix X**T*A*X, stored in the same
61 * format as A.
62 *
63 * LDAB (input) INTEGER
64 * The leading dimension of the array AB. LDAB >= KA+1.
65 *
66 * BB (input) DOUBLE PRECISION array, dimension (LDBB,N)
67 * The banded factor S from the split Cholesky factorization of
68 * B, as returned by DPBSTF, stored in the first KB+1 rows of
69 * the array.
70 *
71 * LDBB (input) INTEGER
72 * The leading dimension of the array BB. LDBB >= KB+1.
73 *
74 * X (output) DOUBLE PRECISION array, dimension (LDX,N)
75 * If VECT = 'V', the n-by-n matrix X.
76 * If VECT = 'N', the array X is not referenced.
77 *
78 * LDX (input) INTEGER
79 * The leading dimension of the array X.
80 * LDX >= max(1,N) if VECT = 'V'; LDX >= 1 otherwise.
81 *
82 * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
83 *
84 * INFO (output) INTEGER
85 * = 0: successful exit
86 * < 0: if INFO = -i, the i-th argument had an illegal value.
87 *
88 * =====================================================================
89 *
90 * .. Parameters ..
91 DOUBLE PRECISION ZERO, ONE
92 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
93 * ..
94 * .. Local Scalars ..
95 LOGICAL UPDATE, UPPER, WANTX
96 INTEGER I, I0, I1, I2, INCA, J, J1, J1T, J2, J2T, K,
97 $ KA1, KB1, KBT, L, M, NR, NRT, NX
98 DOUBLE PRECISION BII, RA, RA1, T
99 * ..
100 * .. External Functions ..
101 LOGICAL LSAME
102 EXTERNAL LSAME
103 * ..
104 * .. External Subroutines ..
105 EXTERNAL DGER, DLAR2V, DLARGV, DLARTG, DLARTV, DLASET,
106 $ DROT, DSCAL, XERBLA
107 * ..
108 * .. Intrinsic Functions ..
109 INTRINSIC MAX, MIN
110 * ..
111 * .. Executable Statements ..
112 *
113 * Test the input parameters
114 *
115 WANTX = LSAME( VECT, 'V' )
116 UPPER = LSAME( UPLO, 'U' )
117 KA1 = KA + 1
118 KB1 = KB + 1
119 INFO = 0
120 IF( .NOT.WANTX .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
121 INFO = -1
122 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
123 INFO = -2
124 ELSE IF( N.LT.0 ) THEN
125 INFO = -3
126 ELSE IF( KA.LT.0 ) THEN
127 INFO = -4
128 ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN
129 INFO = -5
130 ELSE IF( LDAB.LT.KA+1 ) THEN
131 INFO = -7
132 ELSE IF( LDBB.LT.KB+1 ) THEN
133 INFO = -9
134 ELSE IF( LDX.LT.1 .OR. WANTX .AND. LDX.LT.MAX( 1, N ) ) THEN
135 INFO = -11
136 END IF
137 IF( INFO.NE.0 ) THEN
138 CALL XERBLA( 'DSBGST', -INFO )
139 RETURN
140 END IF
141 *
142 * Quick return if possible
143 *
144 IF( N.EQ.0 )
145 $ RETURN
146 *
147 INCA = LDAB*KA1
148 *
149 * Initialize X to the unit matrix, if needed
150 *
151 IF( WANTX )
152 $ CALL DLASET( 'Full', N, N, ZERO, ONE, X, LDX )
153 *
154 * Set M to the splitting point m. It must be the same value as is
155 * used in DPBSTF. The chosen value allows the arrays WORK and RWORK
156 * to be of dimension (N).
157 *
158 M = ( N+KB ) / 2
159 *
160 * The routine works in two phases, corresponding to the two halves
161 * of the split Cholesky factorization of B as S**T*S where
162 *
163 * S = ( U )
164 * ( M L )
165 *
166 * with U upper triangular of order m, and L lower triangular of
167 * order n-m. S has the same bandwidth as B.
168 *
169 * S is treated as a product of elementary matrices:
170 *
171 * S = S(m)*S(m-1)*...*S(2)*S(1)*S(m+1)*S(m+2)*...*S(n-1)*S(n)
172 *
173 * where S(i) is determined by the i-th row of S.
174 *
175 * In phase 1, the index i takes the values n, n-1, ... , m+1;
176 * in phase 2, it takes the values 1, 2, ... , m.
177 *
178 * For each value of i, the current matrix A is updated by forming
179 * inv(S(i))**T*A*inv(S(i)). This creates a triangular bulge outside
180 * the band of A. The bulge is then pushed down toward the bottom of
181 * A in phase 1, and up toward the top of A in phase 2, by applying
182 * plane rotations.
183 *
184 * There are kb*(kb+1)/2 elements in the bulge, but at most 2*kb-1
185 * of them are linearly independent, so annihilating a bulge requires
186 * only 2*kb-1 plane rotations. The rotations are divided into a 1st
187 * set of kb-1 rotations, and a 2nd set of kb rotations.
188 *
189 * Wherever possible, rotations are generated and applied in vector
190 * operations of length NR between the indices J1 and J2 (sometimes
191 * replaced by modified values NRT, J1T or J2T).
192 *
193 * The cosines and sines of the rotations are stored in the array
194 * WORK. The cosines of the 1st set of rotations are stored in
195 * elements n+2:n+m-kb-1 and the sines of the 1st set in elements
196 * 2:m-kb-1; the cosines of the 2nd set are stored in elements
197 * n+m-kb+1:2*n and the sines of the second set in elements m-kb+1:n.
198 *
199 * The bulges are not formed explicitly; nonzero elements outside the
200 * band are created only when they are required for generating new
201 * rotations; they are stored in the array WORK, in positions where
202 * they are later overwritten by the sines of the rotations which
203 * annihilate them.
204 *
205 * **************************** Phase 1 *****************************
206 *
207 * The logical structure of this phase is:
208 *
209 * UPDATE = .TRUE.
210 * DO I = N, M + 1, -1
211 * use S(i) to update A and create a new bulge
212 * apply rotations to push all bulges KA positions downward
213 * END DO
214 * UPDATE = .FALSE.
215 * DO I = M + KA + 1, N - 1
216 * apply rotations to push all bulges KA positions downward
217 * END DO
218 *
219 * To avoid duplicating code, the two loops are merged.
220 *
221 UPDATE = .TRUE.
222 I = N + 1
223 10 CONTINUE
224 IF( UPDATE ) THEN
225 I = I - 1
226 KBT = MIN( KB, I-1 )
227 I0 = I - 1
228 I1 = MIN( N, I+KA )
229 I2 = I - KBT + KA1
230 IF( I.LT.M+1 ) THEN
231 UPDATE = .FALSE.
232 I = I + 1
233 I0 = M
234 IF( KA.EQ.0 )
235 $ GO TO 480
236 GO TO 10
237 END IF
238 ELSE
239 I = I + KA
240 IF( I.GT.N-1 )
241 $ GO TO 480
242 END IF
243 *
244 IF( UPPER ) THEN
245 *
246 * Transform A, working with the upper triangle
247 *
248 IF( UPDATE ) THEN
249 *
250 * Form inv(S(i))**T * A * inv(S(i))
251 *
252 BII = BB( KB1, I )
253 DO 20 J = I, I1
254 AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
255 20 CONTINUE
256 DO 30 J = MAX( 1, I-KA ), I
257 AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
258 30 CONTINUE
259 DO 60 K = I - KBT, I - 1
260 DO 40 J = I - KBT, K
261 AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
262 $ BB( J-I+KB1, I )*AB( K-I+KA1, I ) -
263 $ BB( K-I+KB1, I )*AB( J-I+KA1, I ) +
264 $ AB( KA1, I )*BB( J-I+KB1, I )*
265 $ BB( K-I+KB1, I )
266 40 CONTINUE
267 DO 50 J = MAX( 1, I-KA ), I - KBT - 1
268 AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
269 $ BB( K-I+KB1, I )*AB( J-I+KA1, I )
270 50 CONTINUE
271 60 CONTINUE
272 DO 80 J = I, I1
273 DO 70 K = MAX( J-KA, I-KBT ), I - 1
274 AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
275 $ BB( K-I+KB1, I )*AB( I-J+KA1, J )
276 70 CONTINUE
277 80 CONTINUE
278 *
279 IF( WANTX ) THEN
280 *
281 * post-multiply X by inv(S(i))
282 *
283 CALL DSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
284 IF( KBT.GT.0 )
285 $ CALL DGER( N-M, KBT, -ONE, X( M+1, I ), 1,
286 $ BB( KB1-KBT, I ), 1, X( M+1, I-KBT ), LDX )
287 END IF
288 *
289 * store a(i,i1) in RA1 for use in next loop over K
290 *
291 RA1 = AB( I-I1+KA1, I1 )
292 END IF
293 *
294 * Generate and apply vectors of rotations to chase all the
295 * existing bulges KA positions down toward the bottom of the
296 * band
297 *
298 DO 130 K = 1, KB - 1
299 IF( UPDATE ) THEN
300 *
301 * Determine the rotations which would annihilate the bulge
302 * which has in theory just been created
303 *
304 IF( I-K+KA.LT.N .AND. I-K.GT.1 ) THEN
305 *
306 * generate rotation to annihilate a(i,i-k+ka+1)
307 *
308 CALL DLARTG( AB( K+1, I-K+KA ), RA1,
309 $ WORK( N+I-K+KA-M ), WORK( I-K+KA-M ),
310 $ RA )
311 *
312 * create nonzero element a(i-k,i-k+ka+1) outside the
313 * band and store it in WORK(i-k)
314 *
315 T = -BB( KB1-K, I )*RA1
316 WORK( I-K ) = WORK( N+I-K+KA-M )*T -
317 $ WORK( I-K+KA-M )*AB( 1, I-K+KA )
318 AB( 1, I-K+KA ) = WORK( I-K+KA-M )*T +
319 $ WORK( N+I-K+KA-M )*AB( 1, I-K+KA )
320 RA1 = RA
321 END IF
322 END IF
323 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
324 NR = ( N-J2+KA ) / KA1
325 J1 = J2 + ( NR-1 )*KA1
326 IF( UPDATE ) THEN
327 J2T = MAX( J2, I+2*KA-K+1 )
328 ELSE
329 J2T = J2
330 END IF
331 NRT = ( N-J2T+KA ) / KA1
332 DO 90 J = J2T, J1, KA1
333 *
334 * create nonzero element a(j-ka,j+1) outside the band
335 * and store it in WORK(j-m)
336 *
337 WORK( J-M ) = WORK( J-M )*AB( 1, J+1 )
338 AB( 1, J+1 ) = WORK( N+J-M )*AB( 1, J+1 )
339 90 CONTINUE
340 *
341 * generate rotations in 1st set to annihilate elements which
342 * have been created outside the band
343 *
344 IF( NRT.GT.0 )
345 $ CALL DLARGV( NRT, AB( 1, J2T ), INCA, WORK( J2T-M ), KA1,
346 $ WORK( N+J2T-M ), KA1 )
347 IF( NR.GT.0 ) THEN
348 *
349 * apply rotations in 1st set from the right
350 *
351 DO 100 L = 1, KA - 1
352 CALL DLARTV( NR, AB( KA1-L, J2 ), INCA,
353 $ AB( KA-L, J2+1 ), INCA, WORK( N+J2-M ),
354 $ WORK( J2-M ), KA1 )
355 100 CONTINUE
356 *
357 * apply rotations in 1st set from both sides to diagonal
358 * blocks
359 *
360 CALL DLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
361 $ AB( KA, J2+1 ), INCA, WORK( N+J2-M ),
362 $ WORK( J2-M ), KA1 )
363 *
364 END IF
365 *
366 * start applying rotations in 1st set from the left
367 *
368 DO 110 L = KA - 1, KB - K + 1, -1
369 NRT = ( N-J2+L ) / KA1
370 IF( NRT.GT.0 )
371 $ CALL DLARTV( NRT, AB( L, J2+KA1-L ), INCA,
372 $ AB( L+1, J2+KA1-L ), INCA,
373 $ WORK( N+J2-M ), WORK( J2-M ), KA1 )
374 110 CONTINUE
375 *
376 IF( WANTX ) THEN
377 *
378 * post-multiply X by product of rotations in 1st set
379 *
380 DO 120 J = J2, J1, KA1
381 CALL DROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
382 $ WORK( N+J-M ), WORK( J-M ) )
383 120 CONTINUE
384 END IF
385 130 CONTINUE
386 *
387 IF( UPDATE ) THEN
388 IF( I2.LE.N .AND. KBT.GT.0 ) THEN
389 *
390 * create nonzero element a(i-kbt,i-kbt+ka+1) outside the
391 * band and store it in WORK(i-kbt)
392 *
393 WORK( I-KBT ) = -BB( KB1-KBT, I )*RA1
394 END IF
395 END IF
396 *
397 DO 170 K = KB, 1, -1
398 IF( UPDATE ) THEN
399 J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
400 ELSE
401 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
402 END IF
403 *
404 * finish applying rotations in 2nd set from the left
405 *
406 DO 140 L = KB - K, 1, -1
407 NRT = ( N-J2+KA+L ) / KA1
408 IF( NRT.GT.0 )
409 $ CALL DLARTV( NRT, AB( L, J2-L+1 ), INCA,
410 $ AB( L+1, J2-L+1 ), INCA, WORK( N+J2-KA ),
411 $ WORK( J2-KA ), KA1 )
412 140 CONTINUE
413 NR = ( N-J2+KA ) / KA1
414 J1 = J2 + ( NR-1 )*KA1
415 DO 150 J = J1, J2, -KA1
416 WORK( J ) = WORK( J-KA )
417 WORK( N+J ) = WORK( N+J-KA )
418 150 CONTINUE
419 DO 160 J = J2, J1, KA1
420 *
421 * create nonzero element a(j-ka,j+1) outside the band
422 * and store it in WORK(j)
423 *
424 WORK( J ) = WORK( J )*AB( 1, J+1 )
425 AB( 1, J+1 ) = WORK( N+J )*AB( 1, J+1 )
426 160 CONTINUE
427 IF( UPDATE ) THEN
428 IF( I-K.LT.N-KA .AND. K.LE.KBT )
429 $ WORK( I-K+KA ) = WORK( I-K )
430 END IF
431 170 CONTINUE
432 *
433 DO 210 K = KB, 1, -1
434 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
435 NR = ( N-J2+KA ) / KA1
436 J1 = J2 + ( NR-1 )*KA1
437 IF( NR.GT.0 ) THEN
438 *
439 * generate rotations in 2nd set to annihilate elements
440 * which have been created outside the band
441 *
442 CALL DLARGV( NR, AB( 1, J2 ), INCA, WORK( J2 ), KA1,
443 $ WORK( N+J2 ), KA1 )
444 *
445 * apply rotations in 2nd set from the right
446 *
447 DO 180 L = 1, KA - 1
448 CALL DLARTV( NR, AB( KA1-L, J2 ), INCA,
449 $ AB( KA-L, J2+1 ), INCA, WORK( N+J2 ),
450 $ WORK( J2 ), KA1 )
451 180 CONTINUE
452 *
453 * apply rotations in 2nd set from both sides to diagonal
454 * blocks
455 *
456 CALL DLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
457 $ AB( KA, J2+1 ), INCA, WORK( N+J2 ),
458 $ WORK( J2 ), KA1 )
459 *
460 END IF
461 *
462 * start applying rotations in 2nd set from the left
463 *
464 DO 190 L = KA - 1, KB - K + 1, -1
465 NRT = ( N-J2+L ) / KA1
466 IF( NRT.GT.0 )
467 $ CALL DLARTV( NRT, AB( L, J2+KA1-L ), INCA,
468 $ AB( L+1, J2+KA1-L ), INCA, WORK( N+J2 ),
469 $ WORK( J2 ), KA1 )
470 190 CONTINUE
471 *
472 IF( WANTX ) THEN
473 *
474 * post-multiply X by product of rotations in 2nd set
475 *
476 DO 200 J = J2, J1, KA1
477 CALL DROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
478 $ WORK( N+J ), WORK( J ) )
479 200 CONTINUE
480 END IF
481 210 CONTINUE
482 *
483 DO 230 K = 1, KB - 1
484 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
485 *
486 * finish applying rotations in 1st set from the left
487 *
488 DO 220 L = KB - K, 1, -1
489 NRT = ( N-J2+L ) / KA1
490 IF( NRT.GT.0 )
491 $ CALL DLARTV( NRT, AB( L, J2+KA1-L ), INCA,
492 $ AB( L+1, J2+KA1-L ), INCA,
493 $ WORK( N+J2-M ), WORK( J2-M ), KA1 )
494 220 CONTINUE
495 230 CONTINUE
496 *
497 IF( KB.GT.1 ) THEN
498 DO 240 J = N - 1, I - KB + 2*KA + 1, -1
499 WORK( N+J-M ) = WORK( N+J-KA-M )
500 WORK( J-M ) = WORK( J-KA-M )
501 240 CONTINUE
502 END IF
503 *
504 ELSE
505 *
506 * Transform A, working with the lower triangle
507 *
508 IF( UPDATE ) THEN
509 *
510 * Form inv(S(i))**T * A * inv(S(i))
511 *
512 BII = BB( 1, I )
513 DO 250 J = I, I1
514 AB( J-I+1, I ) = AB( J-I+1, I ) / BII
515 250 CONTINUE
516 DO 260 J = MAX( 1, I-KA ), I
517 AB( I-J+1, J ) = AB( I-J+1, J ) / BII
518 260 CONTINUE
519 DO 290 K = I - KBT, I - 1
520 DO 270 J = I - KBT, K
521 AB( K-J+1, J ) = AB( K-J+1, J ) -
522 $ BB( I-J+1, J )*AB( I-K+1, K ) -
523 $ BB( I-K+1, K )*AB( I-J+1, J ) +
524 $ AB( 1, I )*BB( I-J+1, J )*
525 $ BB( I-K+1, K )
526 270 CONTINUE
527 DO 280 J = MAX( 1, I-KA ), I - KBT - 1
528 AB( K-J+1, J ) = AB( K-J+1, J ) -
529 $ BB( I-K+1, K )*AB( I-J+1, J )
530 280 CONTINUE
531 290 CONTINUE
532 DO 310 J = I, I1
533 DO 300 K = MAX( J-KA, I-KBT ), I - 1
534 AB( J-K+1, K ) = AB( J-K+1, K ) -
535 $ BB( I-K+1, K )*AB( J-I+1, I )
536 300 CONTINUE
537 310 CONTINUE
538 *
539 IF( WANTX ) THEN
540 *
541 * post-multiply X by inv(S(i))
542 *
543 CALL DSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
544 IF( KBT.GT.0 )
545 $ CALL DGER( N-M, KBT, -ONE, X( M+1, I ), 1,
546 $ BB( KBT+1, I-KBT ), LDBB-1,
547 $ X( M+1, I-KBT ), LDX )
548 END IF
549 *
550 * store a(i1,i) in RA1 for use in next loop over K
551 *
552 RA1 = AB( I1-I+1, I )
553 END IF
554 *
555 * Generate and apply vectors of rotations to chase all the
556 * existing bulges KA positions down toward the bottom of the
557 * band
558 *
559 DO 360 K = 1, KB - 1
560 IF( UPDATE ) THEN
561 *
562 * Determine the rotations which would annihilate the bulge
563 * which has in theory just been created
564 *
565 IF( I-K+KA.LT.N .AND. I-K.GT.1 ) THEN
566 *
567 * generate rotation to annihilate a(i-k+ka+1,i)
568 *
569 CALL DLARTG( AB( KA1-K, I ), RA1, WORK( N+I-K+KA-M ),
570 $ WORK( I-K+KA-M ), RA )
571 *
572 * create nonzero element a(i-k+ka+1,i-k) outside the
573 * band and store it in WORK(i-k)
574 *
575 T = -BB( K+1, I-K )*RA1
576 WORK( I-K ) = WORK( N+I-K+KA-M )*T -
577 $ WORK( I-K+KA-M )*AB( KA1, I-K )
578 AB( KA1, I-K ) = WORK( I-K+KA-M )*T +
579 $ WORK( N+I-K+KA-M )*AB( KA1, I-K )
580 RA1 = RA
581 END IF
582 END IF
583 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
584 NR = ( N-J2+KA ) / KA1
585 J1 = J2 + ( NR-1 )*KA1
586 IF( UPDATE ) THEN
587 J2T = MAX( J2, I+2*KA-K+1 )
588 ELSE
589 J2T = J2
590 END IF
591 NRT = ( N-J2T+KA ) / KA1
592 DO 320 J = J2T, J1, KA1
593 *
594 * create nonzero element a(j+1,j-ka) outside the band
595 * and store it in WORK(j-m)
596 *
597 WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 )
598 AB( KA1, J-KA+1 ) = WORK( N+J-M )*AB( KA1, J-KA+1 )
599 320 CONTINUE
600 *
601 * generate rotations in 1st set to annihilate elements which
602 * have been created outside the band
603 *
604 IF( NRT.GT.0 )
605 $ CALL DLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ),
606 $ KA1, WORK( N+J2T-M ), KA1 )
607 IF( NR.GT.0 ) THEN
608 *
609 * apply rotations in 1st set from the left
610 *
611 DO 330 L = 1, KA - 1
612 CALL DLARTV( NR, AB( L+1, J2-L ), INCA,
613 $ AB( L+2, J2-L ), INCA, WORK( N+J2-M ),
614 $ WORK( J2-M ), KA1 )
615 330 CONTINUE
616 *
617 * apply rotations in 1st set from both sides to diagonal
618 * blocks
619 *
620 CALL DLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
621 $ INCA, WORK( N+J2-M ), WORK( J2-M ), KA1 )
622 *
623 END IF
624 *
625 * start applying rotations in 1st set from the right
626 *
627 DO 340 L = KA - 1, KB - K + 1, -1
628 NRT = ( N-J2+L ) / KA1
629 IF( NRT.GT.0 )
630 $ CALL DLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
631 $ AB( KA1-L, J2+1 ), INCA, WORK( N+J2-M ),
632 $ WORK( J2-M ), KA1 )
633 340 CONTINUE
634 *
635 IF( WANTX ) THEN
636 *
637 * post-multiply X by product of rotations in 1st set
638 *
639 DO 350 J = J2, J1, KA1
640 CALL DROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
641 $ WORK( N+J-M ), WORK( J-M ) )
642 350 CONTINUE
643 END IF
644 360 CONTINUE
645 *
646 IF( UPDATE ) THEN
647 IF( I2.LE.N .AND. KBT.GT.0 ) THEN
648 *
649 * create nonzero element a(i-kbt+ka+1,i-kbt) outside the
650 * band and store it in WORK(i-kbt)
651 *
652 WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1
653 END IF
654 END IF
655 *
656 DO 400 K = KB, 1, -1
657 IF( UPDATE ) THEN
658 J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
659 ELSE
660 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
661 END IF
662 *
663 * finish applying rotations in 2nd set from the right
664 *
665 DO 370 L = KB - K, 1, -1
666 NRT = ( N-J2+KA+L ) / KA1
667 IF( NRT.GT.0 )
668 $ CALL DLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA,
669 $ AB( KA1-L, J2-KA+1 ), INCA,
670 $ WORK( N+J2-KA ), WORK( J2-KA ), KA1 )
671 370 CONTINUE
672 NR = ( N-J2+KA ) / KA1
673 J1 = J2 + ( NR-1 )*KA1
674 DO 380 J = J1, J2, -KA1
675 WORK( J ) = WORK( J-KA )
676 WORK( N+J ) = WORK( N+J-KA )
677 380 CONTINUE
678 DO 390 J = J2, J1, KA1
679 *
680 * create nonzero element a(j+1,j-ka) outside the band
681 * and store it in WORK(j)
682 *
683 WORK( J ) = WORK( J )*AB( KA1, J-KA+1 )
684 AB( KA1, J-KA+1 ) = WORK( N+J )*AB( KA1, J-KA+1 )
685 390 CONTINUE
686 IF( UPDATE ) THEN
687 IF( I-K.LT.N-KA .AND. K.LE.KBT )
688 $ WORK( I-K+KA ) = WORK( I-K )
689 END IF
690 400 CONTINUE
691 *
692 DO 440 K = KB, 1, -1
693 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
694 NR = ( N-J2+KA ) / KA1
695 J1 = J2 + ( NR-1 )*KA1
696 IF( NR.GT.0 ) THEN
697 *
698 * generate rotations in 2nd set to annihilate elements
699 * which have been created outside the band
700 *
701 CALL DLARGV( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1,
702 $ WORK( N+J2 ), KA1 )
703 *
704 * apply rotations in 2nd set from the left
705 *
706 DO 410 L = 1, KA - 1
707 CALL DLARTV( NR, AB( L+1, J2-L ), INCA,
708 $ AB( L+2, J2-L ), INCA, WORK( N+J2 ),
709 $ WORK( J2 ), KA1 )
710 410 CONTINUE
711 *
712 * apply rotations in 2nd set from both sides to diagonal
713 * blocks
714 *
715 CALL DLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
716 $ INCA, WORK( N+J2 ), WORK( J2 ), KA1 )
717 *
718 END IF
719 *
720 * start applying rotations in 2nd set from the right
721 *
722 DO 420 L = KA - 1, KB - K + 1, -1
723 NRT = ( N-J2+L ) / KA1
724 IF( NRT.GT.0 )
725 $ CALL DLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
726 $ AB( KA1-L, J2+1 ), INCA, WORK( N+J2 ),
727 $ WORK( J2 ), KA1 )
728 420 CONTINUE
729 *
730 IF( WANTX ) THEN
731 *
732 * post-multiply X by product of rotations in 2nd set
733 *
734 DO 430 J = J2, J1, KA1
735 CALL DROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
736 $ WORK( N+J ), WORK( J ) )
737 430 CONTINUE
738 END IF
739 440 CONTINUE
740 *
741 DO 460 K = 1, KB - 1
742 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
743 *
744 * finish applying rotations in 1st set from the right
745 *
746 DO 450 L = KB - K, 1, -1
747 NRT = ( N-J2+L ) / KA1
748 IF( NRT.GT.0 )
749 $ CALL DLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
750 $ AB( KA1-L, J2+1 ), INCA, WORK( N+J2-M ),
751 $ WORK( J2-M ), KA1 )
752 450 CONTINUE
753 460 CONTINUE
754 *
755 IF( KB.GT.1 ) THEN
756 DO 470 J = N - 1, I - KB + 2*KA + 1, -1
757 WORK( N+J-M ) = WORK( N+J-KA-M )
758 WORK( J-M ) = WORK( J-KA-M )
759 470 CONTINUE
760 END IF
761 *
762 END IF
763 *
764 GO TO 10
765 *
766 480 CONTINUE
767 *
768 * **************************** Phase 2 *****************************
769 *
770 * The logical structure of this phase is:
771 *
772 * UPDATE = .TRUE.
773 * DO I = 1, M
774 * use S(i) to update A and create a new bulge
775 * apply rotations to push all bulges KA positions upward
776 * END DO
777 * UPDATE = .FALSE.
778 * DO I = M - KA - 1, 2, -1
779 * apply rotations to push all bulges KA positions upward
780 * END DO
781 *
782 * To avoid duplicating code, the two loops are merged.
783 *
784 UPDATE = .TRUE.
785 I = 0
786 490 CONTINUE
787 IF( UPDATE ) THEN
788 I = I + 1
789 KBT = MIN( KB, M-I )
790 I0 = I + 1
791 I1 = MAX( 1, I-KA )
792 I2 = I + KBT - KA1
793 IF( I.GT.M ) THEN
794 UPDATE = .FALSE.
795 I = I - 1
796 I0 = M + 1
797 IF( KA.EQ.0 )
798 $ RETURN
799 GO TO 490
800 END IF
801 ELSE
802 I = I - KA
803 IF( I.LT.2 )
804 $ RETURN
805 END IF
806 *
807 IF( I.LT.M-KBT ) THEN
808 NX = M
809 ELSE
810 NX = N
811 END IF
812 *
813 IF( UPPER ) THEN
814 *
815 * Transform A, working with the upper triangle
816 *
817 IF( UPDATE ) THEN
818 *
819 * Form inv(S(i))**T * A * inv(S(i))
820 *
821 BII = BB( KB1, I )
822 DO 500 J = I1, I
823 AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
824 500 CONTINUE
825 DO 510 J = I, MIN( N, I+KA )
826 AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
827 510 CONTINUE
828 DO 540 K = I + 1, I + KBT
829 DO 520 J = K, I + KBT
830 AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
831 $ BB( I-J+KB1, J )*AB( I-K+KA1, K ) -
832 $ BB( I-K+KB1, K )*AB( I-J+KA1, J ) +
833 $ AB( KA1, I )*BB( I-J+KB1, J )*
834 $ BB( I-K+KB1, K )
835 520 CONTINUE
836 DO 530 J = I + KBT + 1, MIN( N, I+KA )
837 AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
838 $ BB( I-K+KB1, K )*AB( I-J+KA1, J )
839 530 CONTINUE
840 540 CONTINUE
841 DO 560 J = I1, I
842 DO 550 K = I + 1, MIN( J+KA, I+KBT )
843 AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
844 $ BB( I-K+KB1, K )*AB( J-I+KA1, I )
845 550 CONTINUE
846 560 CONTINUE
847 *
848 IF( WANTX ) THEN
849 *
850 * post-multiply X by inv(S(i))
851 *
852 CALL DSCAL( NX, ONE / BII, X( 1, I ), 1 )
853 IF( KBT.GT.0 )
854 $ CALL DGER( NX, KBT, -ONE, X( 1, I ), 1, BB( KB, I+1 ),
855 $ LDBB-1, X( 1, I+1 ), LDX )
856 END IF
857 *
858 * store a(i1,i) in RA1 for use in next loop over K
859 *
860 RA1 = AB( I1-I+KA1, I )
861 END IF
862 *
863 * Generate and apply vectors of rotations to chase all the
864 * existing bulges KA positions up toward the top of the band
865 *
866 DO 610 K = 1, KB - 1
867 IF( UPDATE ) THEN
868 *
869 * Determine the rotations which would annihilate the bulge
870 * which has in theory just been created
871 *
872 IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
873 *
874 * generate rotation to annihilate a(i+k-ka-1,i)
875 *
876 CALL DLARTG( AB( K+1, I ), RA1, WORK( N+I+K-KA ),
877 $ WORK( I+K-KA ), RA )
878 *
879 * create nonzero element a(i+k-ka-1,i+k) outside the
880 * band and store it in WORK(m-kb+i+k)
881 *
882 T = -BB( KB1-K, I+K )*RA1
883 WORK( M-KB+I+K ) = WORK( N+I+K-KA )*T -
884 $ WORK( I+K-KA )*AB( 1, I+K )
885 AB( 1, I+K ) = WORK( I+K-KA )*T +
886 $ WORK( N+I+K-KA )*AB( 1, I+K )
887 RA1 = RA
888 END IF
889 END IF
890 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
891 NR = ( J2+KA-1 ) / KA1
892 J1 = J2 - ( NR-1 )*KA1
893 IF( UPDATE ) THEN
894 J2T = MIN( J2, I-2*KA+K-1 )
895 ELSE
896 J2T = J2
897 END IF
898 NRT = ( J2T+KA-1 ) / KA1
899 DO 570 J = J1, J2T, KA1
900 *
901 * create nonzero element a(j-1,j+ka) outside the band
902 * and store it in WORK(j)
903 *
904 WORK( J ) = WORK( J )*AB( 1, J+KA-1 )
905 AB( 1, J+KA-1 ) = WORK( N+J )*AB( 1, J+KA-1 )
906 570 CONTINUE
907 *
908 * generate rotations in 1st set to annihilate elements which
909 * have been created outside the band
910 *
911 IF( NRT.GT.0 )
912 $ CALL DLARGV( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1,
913 $ WORK( N+J1 ), KA1 )
914 IF( NR.GT.0 ) THEN
915 *
916 * apply rotations in 1st set from the left
917 *
918 DO 580 L = 1, KA - 1
919 CALL DLARTV( NR, AB( KA1-L, J1+L ), INCA,
920 $ AB( KA-L, J1+L ), INCA, WORK( N+J1 ),
921 $ WORK( J1 ), KA1 )
922 580 CONTINUE
923 *
924 * apply rotations in 1st set from both sides to diagonal
925 * blocks
926 *
927 CALL DLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
928 $ AB( KA, J1 ), INCA, WORK( N+J1 ),
929 $ WORK( J1 ), KA1 )
930 *
931 END IF
932 *
933 * start applying rotations in 1st set from the right
934 *
935 DO 590 L = KA - 1, KB - K + 1, -1
936 NRT = ( J2+L-1 ) / KA1
937 J1T = J2 - ( NRT-1 )*KA1
938 IF( NRT.GT.0 )
939 $ CALL DLARTV( NRT, AB( L, J1T ), INCA,
940 $ AB( L+1, J1T-1 ), INCA, WORK( N+J1T ),
941 $ WORK( J1T ), KA1 )
942 590 CONTINUE
943 *
944 IF( WANTX ) THEN
945 *
946 * post-multiply X by product of rotations in 1st set
947 *
948 DO 600 J = J1, J2, KA1
949 CALL DROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
950 $ WORK( N+J ), WORK( J ) )
951 600 CONTINUE
952 END IF
953 610 CONTINUE
954 *
955 IF( UPDATE ) THEN
956 IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
957 *
958 * create nonzero element a(i+kbt-ka-1,i+kbt) outside the
959 * band and store it in WORK(m-kb+i+kbt)
960 *
961 WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1
962 END IF
963 END IF
964 *
965 DO 650 K = KB, 1, -1
966 IF( UPDATE ) THEN
967 J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
968 ELSE
969 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
970 END IF
971 *
972 * finish applying rotations in 2nd set from the right
973 *
974 DO 620 L = KB - K, 1, -1
975 NRT = ( J2+KA+L-1 ) / KA1
976 J1T = J2 - ( NRT-1 )*KA1
977 IF( NRT.GT.0 )
978 $ CALL DLARTV( NRT, AB( L, J1T+KA ), INCA,
979 $ AB( L+1, J1T+KA-1 ), INCA,
980 $ WORK( N+M-KB+J1T+KA ),
981 $ WORK( M-KB+J1T+KA ), KA1 )
982 620 CONTINUE
983 NR = ( J2+KA-1 ) / KA1
984 J1 = J2 - ( NR-1 )*KA1
985 DO 630 J = J1, J2, KA1
986 WORK( M-KB+J ) = WORK( M-KB+J+KA )
987 WORK( N+M-KB+J ) = WORK( N+M-KB+J+KA )
988 630 CONTINUE
989 DO 640 J = J1, J2, KA1
990 *
991 * create nonzero element a(j-1,j+ka) outside the band
992 * and store it in WORK(m-kb+j)
993 *
994 WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 )
995 AB( 1, J+KA-1 ) = WORK( N+M-KB+J )*AB( 1, J+KA-1 )
996 640 CONTINUE
997 IF( UPDATE ) THEN
998 IF( I+K.GT.KA1 .AND. K.LE.KBT )
999 $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
1000 END IF
1001 650 CONTINUE
1002 *
1003 DO 690 K = KB, 1, -1
1004 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1005 NR = ( J2+KA-1 ) / KA1
1006 J1 = J2 - ( NR-1 )*KA1
1007 IF( NR.GT.0 ) THEN
1008 *
1009 * generate rotations in 2nd set to annihilate elements
1010 * which have been created outside the band
1011 *
1012 CALL DLARGV( NR, AB( 1, J1+KA ), INCA, WORK( M-KB+J1 ),
1013 $ KA1, WORK( N+M-KB+J1 ), KA1 )
1014 *
1015 * apply rotations in 2nd set from the left
1016 *
1017 DO 660 L = 1, KA - 1
1018 CALL DLARTV( NR, AB( KA1-L, J1+L ), INCA,
1019 $ AB( KA-L, J1+L ), INCA,
1020 $ WORK( N+M-KB+J1 ), WORK( M-KB+J1 ), KA1 )
1021 660 CONTINUE
1022 *
1023 * apply rotations in 2nd set from both sides to diagonal
1024 * blocks
1025 *
1026 CALL DLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
1027 $ AB( KA, J1 ), INCA, WORK( N+M-KB+J1 ),
1028 $ WORK( M-KB+J1 ), KA1 )
1029 *
1030 END IF
1031 *
1032 * start applying rotations in 2nd set from the right
1033 *
1034 DO 670 L = KA - 1, KB - K + 1, -1
1035 NRT = ( J2+L-1 ) / KA1
1036 J1T = J2 - ( NRT-1 )*KA1
1037 IF( NRT.GT.0 )
1038 $ CALL DLARTV( NRT, AB( L, J1T ), INCA,
1039 $ AB( L+1, J1T-1 ), INCA,
1040 $ WORK( N+M-KB+J1T ), WORK( M-KB+J1T ),
1041 $ KA1 )
1042 670 CONTINUE
1043 *
1044 IF( WANTX ) THEN
1045 *
1046 * post-multiply X by product of rotations in 2nd set
1047 *
1048 DO 680 J = J1, J2, KA1
1049 CALL DROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1050 $ WORK( N+M-KB+J ), WORK( M-KB+J ) )
1051 680 CONTINUE
1052 END IF
1053 690 CONTINUE
1054 *
1055 DO 710 K = 1, KB - 1
1056 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1057 *
1058 * finish applying rotations in 1st set from the right
1059 *
1060 DO 700 L = KB - K, 1, -1
1061 NRT = ( J2+L-1 ) / KA1
1062 J1T = J2 - ( NRT-1 )*KA1
1063 IF( NRT.GT.0 )
1064 $ CALL DLARTV( NRT, AB( L, J1T ), INCA,
1065 $ AB( L+1, J1T-1 ), INCA, WORK( N+J1T ),
1066 $ WORK( J1T ), KA1 )
1067 700 CONTINUE
1068 710 CONTINUE
1069 *
1070 IF( KB.GT.1 ) THEN
1071 DO 720 J = 2, MIN( I+KB, M ) - 2*KA - 1
1072 WORK( N+J ) = WORK( N+J+KA )
1073 WORK( J ) = WORK( J+KA )
1074 720 CONTINUE
1075 END IF
1076 *
1077 ELSE
1078 *
1079 * Transform A, working with the lower triangle
1080 *
1081 IF( UPDATE ) THEN
1082 *
1083 * Form inv(S(i))**T * A * inv(S(i))
1084 *
1085 BII = BB( 1, I )
1086 DO 730 J = I1, I
1087 AB( I-J+1, J ) = AB( I-J+1, J ) / BII
1088 730 CONTINUE
1089 DO 740 J = I, MIN( N, I+KA )
1090 AB( J-I+1, I ) = AB( J-I+1, I ) / BII
1091 740 CONTINUE
1092 DO 770 K = I + 1, I + KBT
1093 DO 750 J = K, I + KBT
1094 AB( J-K+1, K ) = AB( J-K+1, K ) -
1095 $ BB( J-I+1, I )*AB( K-I+1, I ) -
1096 $ BB( K-I+1, I )*AB( J-I+1, I ) +
1097 $ AB( 1, I )*BB( J-I+1, I )*
1098 $ BB( K-I+1, I )
1099 750 CONTINUE
1100 DO 760 J = I + KBT + 1, MIN( N, I+KA )
1101 AB( J-K+1, K ) = AB( J-K+1, K ) -
1102 $ BB( K-I+1, I )*AB( J-I+1, I )
1103 760 CONTINUE
1104 770 CONTINUE
1105 DO 790 J = I1, I
1106 DO 780 K = I + 1, MIN( J+KA, I+KBT )
1107 AB( K-J+1, J ) = AB( K-J+1, J ) -
1108 $ BB( K-I+1, I )*AB( I-J+1, J )
1109 780 CONTINUE
1110 790 CONTINUE
1111 *
1112 IF( WANTX ) THEN
1113 *
1114 * post-multiply X by inv(S(i))
1115 *
1116 CALL DSCAL( NX, ONE / BII, X( 1, I ), 1 )
1117 IF( KBT.GT.0 )
1118 $ CALL DGER( NX, KBT, -ONE, X( 1, I ), 1, BB( 2, I ), 1,
1119 $ X( 1, I+1 ), LDX )
1120 END IF
1121 *
1122 * store a(i,i1) in RA1 for use in next loop over K
1123 *
1124 RA1 = AB( I-I1+1, I1 )
1125 END IF
1126 *
1127 * Generate and apply vectors of rotations to chase all the
1128 * existing bulges KA positions up toward the top of the band
1129 *
1130 DO 840 K = 1, KB - 1
1131 IF( UPDATE ) THEN
1132 *
1133 * Determine the rotations which would annihilate the bulge
1134 * which has in theory just been created
1135 *
1136 IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
1137 *
1138 * generate rotation to annihilate a(i,i+k-ka-1)
1139 *
1140 CALL DLARTG( AB( KA1-K, I+K-KA ), RA1,
1141 $ WORK( N+I+K-KA ), WORK( I+K-KA ), RA )
1142 *
1143 * create nonzero element a(i+k,i+k-ka-1) outside the
1144 * band and store it in WORK(m-kb+i+k)
1145 *
1146 T = -BB( K+1, I )*RA1
1147 WORK( M-KB+I+K ) = WORK( N+I+K-KA )*T -
1148 $ WORK( I+K-KA )*AB( KA1, I+K-KA )
1149 AB( KA1, I+K-KA ) = WORK( I+K-KA )*T +
1150 $ WORK( N+I+K-KA )*AB( KA1, I+K-KA )
1151 RA1 = RA
1152 END IF
1153 END IF
1154 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1155 NR = ( J2+KA-1 ) / KA1
1156 J1 = J2 - ( NR-1 )*KA1
1157 IF( UPDATE ) THEN
1158 J2T = MIN( J2, I-2*KA+K-1 )
1159 ELSE
1160 J2T = J2
1161 END IF
1162 NRT = ( J2T+KA-1 ) / KA1
1163 DO 800 J = J1, J2T, KA1
1164 *
1165 * create nonzero element a(j+ka,j-1) outside the band
1166 * and store it in WORK(j)
1167 *
1168 WORK( J ) = WORK( J )*AB( KA1, J-1 )
1169 AB( KA1, J-1 ) = WORK( N+J )*AB( KA1, J-1 )
1170 800 CONTINUE
1171 *
1172 * generate rotations in 1st set to annihilate elements which
1173 * have been created outside the band
1174 *
1175 IF( NRT.GT.0 )
1176 $ CALL DLARGV( NRT, AB( KA1, J1 ), INCA, WORK( J1 ), KA1,
1177 $ WORK( N+J1 ), KA1 )
1178 IF( NR.GT.0 ) THEN
1179 *
1180 * apply rotations in 1st set from the right
1181 *
1182 DO 810 L = 1, KA - 1
1183 CALL DLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
1184 $ INCA, WORK( N+J1 ), WORK( J1 ), KA1 )
1185 810 CONTINUE
1186 *
1187 * apply rotations in 1st set from both sides to diagonal
1188 * blocks
1189 *
1190 CALL DLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
1191 $ AB( 2, J1-1 ), INCA, WORK( N+J1 ),
1192 $ WORK( J1 ), KA1 )
1193 *
1194 END IF
1195 *
1196 * start applying rotations in 1st set from the left
1197 *
1198 DO 820 L = KA - 1, KB - K + 1, -1
1199 NRT = ( J2+L-1 ) / KA1
1200 J1T = J2 - ( NRT-1 )*KA1
1201 IF( NRT.GT.0 )
1202 $ CALL DLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1203 $ AB( KA1-L, J1T-KA1+L ), INCA,
1204 $ WORK( N+J1T ), WORK( J1T ), KA1 )
1205 820 CONTINUE
1206 *
1207 IF( WANTX ) THEN
1208 *
1209 * post-multiply X by product of rotations in 1st set
1210 *
1211 DO 830 J = J1, J2, KA1
1212 CALL DROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1213 $ WORK( N+J ), WORK( J ) )
1214 830 CONTINUE
1215 END IF
1216 840 CONTINUE
1217 *
1218 IF( UPDATE ) THEN
1219 IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
1220 *
1221 * create nonzero element a(i+kbt,i+kbt-ka-1) outside the
1222 * band and store it in WORK(m-kb+i+kbt)
1223 *
1224 WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1
1225 END IF
1226 END IF
1227 *
1228 DO 880 K = KB, 1, -1
1229 IF( UPDATE ) THEN
1230 J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
1231 ELSE
1232 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1233 END IF
1234 *
1235 * finish applying rotations in 2nd set from the left
1236 *
1237 DO 850 L = KB - K, 1, -1
1238 NRT = ( J2+KA+L-1 ) / KA1
1239 J1T = J2 - ( NRT-1 )*KA1
1240 IF( NRT.GT.0 )
1241 $ CALL DLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA,
1242 $ AB( KA1-L, J1T+L-1 ), INCA,
1243 $ WORK( N+M-KB+J1T+KA ),
1244 $ WORK( M-KB+J1T+KA ), KA1 )
1245 850 CONTINUE
1246 NR = ( J2+KA-1 ) / KA1
1247 J1 = J2 - ( NR-1 )*KA1
1248 DO 860 J = J1, J2, KA1
1249 WORK( M-KB+J ) = WORK( M-KB+J+KA )
1250 WORK( N+M-KB+J ) = WORK( N+M-KB+J+KA )
1251 860 CONTINUE
1252 DO 870 J = J1, J2, KA1
1253 *
1254 * create nonzero element a(j+ka,j-1) outside the band
1255 * and store it in WORK(m-kb+j)
1256 *
1257 WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 )
1258 AB( KA1, J-1 ) = WORK( N+M-KB+J )*AB( KA1, J-1 )
1259 870 CONTINUE
1260 IF( UPDATE ) THEN
1261 IF( I+K.GT.KA1 .AND. K.LE.KBT )
1262 $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
1263 END IF
1264 880 CONTINUE
1265 *
1266 DO 920 K = KB, 1, -1
1267 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1268 NR = ( J2+KA-1 ) / KA1
1269 J1 = J2 - ( NR-1 )*KA1
1270 IF( NR.GT.0 ) THEN
1271 *
1272 * generate rotations in 2nd set to annihilate elements
1273 * which have been created outside the band
1274 *
1275 CALL DLARGV( NR, AB( KA1, J1 ), INCA, WORK( M-KB+J1 ),
1276 $ KA1, WORK( N+M-KB+J1 ), KA1 )
1277 *
1278 * apply rotations in 2nd set from the right
1279 *
1280 DO 890 L = 1, KA - 1
1281 CALL DLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
1282 $ INCA, WORK( N+M-KB+J1 ), WORK( M-KB+J1 ),
1283 $ KA1 )
1284 890 CONTINUE
1285 *
1286 * apply rotations in 2nd set from both sides to diagonal
1287 * blocks
1288 *
1289 CALL DLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
1290 $ AB( 2, J1-1 ), INCA, WORK( N+M-KB+J1 ),
1291 $ WORK( M-KB+J1 ), KA1 )
1292 *
1293 END IF
1294 *
1295 * start applying rotations in 2nd set from the left
1296 *
1297 DO 900 L = KA - 1, KB - K + 1, -1
1298 NRT = ( J2+L-1 ) / KA1
1299 J1T = J2 - ( NRT-1 )*KA1
1300 IF( NRT.GT.0 )
1301 $ CALL DLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1302 $ AB( KA1-L, J1T-KA1+L ), INCA,
1303 $ WORK( N+M-KB+J1T ), WORK( M-KB+J1T ),
1304 $ KA1 )
1305 900 CONTINUE
1306 *
1307 IF( WANTX ) THEN
1308 *
1309 * post-multiply X by product of rotations in 2nd set
1310 *
1311 DO 910 J = J1, J2, KA1
1312 CALL DROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1313 $ WORK( N+M-KB+J ), WORK( M-KB+J ) )
1314 910 CONTINUE
1315 END IF
1316 920 CONTINUE
1317 *
1318 DO 940 K = 1, KB - 1
1319 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1320 *
1321 * finish applying rotations in 1st set from the left
1322 *
1323 DO 930 L = KB - K, 1, -1
1324 NRT = ( J2+L-1 ) / KA1
1325 J1T = J2 - ( NRT-1 )*KA1
1326 IF( NRT.GT.0 )
1327 $ CALL DLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1328 $ AB( KA1-L, J1T-KA1+L ), INCA,
1329 $ WORK( N+J1T ), WORK( J1T ), KA1 )
1330 930 CONTINUE
1331 940 CONTINUE
1332 *
1333 IF( KB.GT.1 ) THEN
1334 DO 950 J = 2, MIN( I+KB, M ) - 2*KA - 1
1335 WORK( N+J ) = WORK( N+J+KA )
1336 WORK( J ) = WORK( J+KA )
1337 950 CONTINUE
1338 END IF
1339 *
1340 END IF
1341 *
1342 GO TO 490
1343 *
1344 * End of DSBGST
1345 *
1346 END
2 $ LDX, WORK, INFO )
3 *
4 * -- LAPACK routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER UPLO, VECT
11 INTEGER INFO, KA, KB, LDAB, LDBB, LDX, N
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
15 $ X( LDX, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DSBGST reduces a real symmetric-definite banded generalized
22 * eigenproblem A*x = lambda*B*x to standard form C*y = lambda*y,
23 * such that C has the same bandwidth as A.
24 *
25 * B must have been previously factorized as S**T*S by DPBSTF, using a
26 * split Cholesky factorization. A is overwritten by C = X**T*A*X, where
27 * X = S**(-1)*Q and Q is an orthogonal matrix chosen to preserve the
28 * bandwidth of A.
29 *
30 * Arguments
31 * =========
32 *
33 * VECT (input) CHARACTER*1
34 * = 'N': do not form the transformation matrix X;
35 * = 'V': form X.
36 *
37 * UPLO (input) CHARACTER*1
38 * = 'U': Upper triangle of A is stored;
39 * = 'L': Lower triangle of A is stored.
40 *
41 * N (input) INTEGER
42 * The order of the matrices A and B. N >= 0.
43 *
44 * KA (input) INTEGER
45 * The number of superdiagonals of the matrix A if UPLO = 'U',
46 * or the number of subdiagonals if UPLO = 'L'. KA >= 0.
47 *
48 * KB (input) INTEGER
49 * The number of superdiagonals of the matrix B if UPLO = 'U',
50 * or the number of subdiagonals if UPLO = 'L'. KA >= KB >= 0.
51 *
52 * AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
53 * On entry, the upper or lower triangle of the symmetric band
54 * matrix A, stored in the first ka+1 rows of the array. The
55 * j-th column of A is stored in the j-th column of the array AB
56 * as follows:
57 * if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
58 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+ka).
59 *
60 * On exit, the transformed matrix X**T*A*X, stored in the same
61 * format as A.
62 *
63 * LDAB (input) INTEGER
64 * The leading dimension of the array AB. LDAB >= KA+1.
65 *
66 * BB (input) DOUBLE PRECISION array, dimension (LDBB,N)
67 * The banded factor S from the split Cholesky factorization of
68 * B, as returned by DPBSTF, stored in the first KB+1 rows of
69 * the array.
70 *
71 * LDBB (input) INTEGER
72 * The leading dimension of the array BB. LDBB >= KB+1.
73 *
74 * X (output) DOUBLE PRECISION array, dimension (LDX,N)
75 * If VECT = 'V', the n-by-n matrix X.
76 * If VECT = 'N', the array X is not referenced.
77 *
78 * LDX (input) INTEGER
79 * The leading dimension of the array X.
80 * LDX >= max(1,N) if VECT = 'V'; LDX >= 1 otherwise.
81 *
82 * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
83 *
84 * INFO (output) INTEGER
85 * = 0: successful exit
86 * < 0: if INFO = -i, the i-th argument had an illegal value.
87 *
88 * =====================================================================
89 *
90 * .. Parameters ..
91 DOUBLE PRECISION ZERO, ONE
92 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
93 * ..
94 * .. Local Scalars ..
95 LOGICAL UPDATE, UPPER, WANTX
96 INTEGER I, I0, I1, I2, INCA, J, J1, J1T, J2, J2T, K,
97 $ KA1, KB1, KBT, L, M, NR, NRT, NX
98 DOUBLE PRECISION BII, RA, RA1, T
99 * ..
100 * .. External Functions ..
101 LOGICAL LSAME
102 EXTERNAL LSAME
103 * ..
104 * .. External Subroutines ..
105 EXTERNAL DGER, DLAR2V, DLARGV, DLARTG, DLARTV, DLASET,
106 $ DROT, DSCAL, XERBLA
107 * ..
108 * .. Intrinsic Functions ..
109 INTRINSIC MAX, MIN
110 * ..
111 * .. Executable Statements ..
112 *
113 * Test the input parameters
114 *
115 WANTX = LSAME( VECT, 'V' )
116 UPPER = LSAME( UPLO, 'U' )
117 KA1 = KA + 1
118 KB1 = KB + 1
119 INFO = 0
120 IF( .NOT.WANTX .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
121 INFO = -1
122 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
123 INFO = -2
124 ELSE IF( N.LT.0 ) THEN
125 INFO = -3
126 ELSE IF( KA.LT.0 ) THEN
127 INFO = -4
128 ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN
129 INFO = -5
130 ELSE IF( LDAB.LT.KA+1 ) THEN
131 INFO = -7
132 ELSE IF( LDBB.LT.KB+1 ) THEN
133 INFO = -9
134 ELSE IF( LDX.LT.1 .OR. WANTX .AND. LDX.LT.MAX( 1, N ) ) THEN
135 INFO = -11
136 END IF
137 IF( INFO.NE.0 ) THEN
138 CALL XERBLA( 'DSBGST', -INFO )
139 RETURN
140 END IF
141 *
142 * Quick return if possible
143 *
144 IF( N.EQ.0 )
145 $ RETURN
146 *
147 INCA = LDAB*KA1
148 *
149 * Initialize X to the unit matrix, if needed
150 *
151 IF( WANTX )
152 $ CALL DLASET( 'Full', N, N, ZERO, ONE, X, LDX )
153 *
154 * Set M to the splitting point m. It must be the same value as is
155 * used in DPBSTF. The chosen value allows the arrays WORK and RWORK
156 * to be of dimension (N).
157 *
158 M = ( N+KB ) / 2
159 *
160 * The routine works in two phases, corresponding to the two halves
161 * of the split Cholesky factorization of B as S**T*S where
162 *
163 * S = ( U )
164 * ( M L )
165 *
166 * with U upper triangular of order m, and L lower triangular of
167 * order n-m. S has the same bandwidth as B.
168 *
169 * S is treated as a product of elementary matrices:
170 *
171 * S = S(m)*S(m-1)*...*S(2)*S(1)*S(m+1)*S(m+2)*...*S(n-1)*S(n)
172 *
173 * where S(i) is determined by the i-th row of S.
174 *
175 * In phase 1, the index i takes the values n, n-1, ... , m+1;
176 * in phase 2, it takes the values 1, 2, ... , m.
177 *
178 * For each value of i, the current matrix A is updated by forming
179 * inv(S(i))**T*A*inv(S(i)). This creates a triangular bulge outside
180 * the band of A. The bulge is then pushed down toward the bottom of
181 * A in phase 1, and up toward the top of A in phase 2, by applying
182 * plane rotations.
183 *
184 * There are kb*(kb+1)/2 elements in the bulge, but at most 2*kb-1
185 * of them are linearly independent, so annihilating a bulge requires
186 * only 2*kb-1 plane rotations. The rotations are divided into a 1st
187 * set of kb-1 rotations, and a 2nd set of kb rotations.
188 *
189 * Wherever possible, rotations are generated and applied in vector
190 * operations of length NR between the indices J1 and J2 (sometimes
191 * replaced by modified values NRT, J1T or J2T).
192 *
193 * The cosines and sines of the rotations are stored in the array
194 * WORK. The cosines of the 1st set of rotations are stored in
195 * elements n+2:n+m-kb-1 and the sines of the 1st set in elements
196 * 2:m-kb-1; the cosines of the 2nd set are stored in elements
197 * n+m-kb+1:2*n and the sines of the second set in elements m-kb+1:n.
198 *
199 * The bulges are not formed explicitly; nonzero elements outside the
200 * band are created only when they are required for generating new
201 * rotations; they are stored in the array WORK, in positions where
202 * they are later overwritten by the sines of the rotations which
203 * annihilate them.
204 *
205 * **************************** Phase 1 *****************************
206 *
207 * The logical structure of this phase is:
208 *
209 * UPDATE = .TRUE.
210 * DO I = N, M + 1, -1
211 * use S(i) to update A and create a new bulge
212 * apply rotations to push all bulges KA positions downward
213 * END DO
214 * UPDATE = .FALSE.
215 * DO I = M + KA + 1, N - 1
216 * apply rotations to push all bulges KA positions downward
217 * END DO
218 *
219 * To avoid duplicating code, the two loops are merged.
220 *
221 UPDATE = .TRUE.
222 I = N + 1
223 10 CONTINUE
224 IF( UPDATE ) THEN
225 I = I - 1
226 KBT = MIN( KB, I-1 )
227 I0 = I - 1
228 I1 = MIN( N, I+KA )
229 I2 = I - KBT + KA1
230 IF( I.LT.M+1 ) THEN
231 UPDATE = .FALSE.
232 I = I + 1
233 I0 = M
234 IF( KA.EQ.0 )
235 $ GO TO 480
236 GO TO 10
237 END IF
238 ELSE
239 I = I + KA
240 IF( I.GT.N-1 )
241 $ GO TO 480
242 END IF
243 *
244 IF( UPPER ) THEN
245 *
246 * Transform A, working with the upper triangle
247 *
248 IF( UPDATE ) THEN
249 *
250 * Form inv(S(i))**T * A * inv(S(i))
251 *
252 BII = BB( KB1, I )
253 DO 20 J = I, I1
254 AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
255 20 CONTINUE
256 DO 30 J = MAX( 1, I-KA ), I
257 AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
258 30 CONTINUE
259 DO 60 K = I - KBT, I - 1
260 DO 40 J = I - KBT, K
261 AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
262 $ BB( J-I+KB1, I )*AB( K-I+KA1, I ) -
263 $ BB( K-I+KB1, I )*AB( J-I+KA1, I ) +
264 $ AB( KA1, I )*BB( J-I+KB1, I )*
265 $ BB( K-I+KB1, I )
266 40 CONTINUE
267 DO 50 J = MAX( 1, I-KA ), I - KBT - 1
268 AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
269 $ BB( K-I+KB1, I )*AB( J-I+KA1, I )
270 50 CONTINUE
271 60 CONTINUE
272 DO 80 J = I, I1
273 DO 70 K = MAX( J-KA, I-KBT ), I - 1
274 AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
275 $ BB( K-I+KB1, I )*AB( I-J+KA1, J )
276 70 CONTINUE
277 80 CONTINUE
278 *
279 IF( WANTX ) THEN
280 *
281 * post-multiply X by inv(S(i))
282 *
283 CALL DSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
284 IF( KBT.GT.0 )
285 $ CALL DGER( N-M, KBT, -ONE, X( M+1, I ), 1,
286 $ BB( KB1-KBT, I ), 1, X( M+1, I-KBT ), LDX )
287 END IF
288 *
289 * store a(i,i1) in RA1 for use in next loop over K
290 *
291 RA1 = AB( I-I1+KA1, I1 )
292 END IF
293 *
294 * Generate and apply vectors of rotations to chase all the
295 * existing bulges KA positions down toward the bottom of the
296 * band
297 *
298 DO 130 K = 1, KB - 1
299 IF( UPDATE ) THEN
300 *
301 * Determine the rotations which would annihilate the bulge
302 * which has in theory just been created
303 *
304 IF( I-K+KA.LT.N .AND. I-K.GT.1 ) THEN
305 *
306 * generate rotation to annihilate a(i,i-k+ka+1)
307 *
308 CALL DLARTG( AB( K+1, I-K+KA ), RA1,
309 $ WORK( N+I-K+KA-M ), WORK( I-K+KA-M ),
310 $ RA )
311 *
312 * create nonzero element a(i-k,i-k+ka+1) outside the
313 * band and store it in WORK(i-k)
314 *
315 T = -BB( KB1-K, I )*RA1
316 WORK( I-K ) = WORK( N+I-K+KA-M )*T -
317 $ WORK( I-K+KA-M )*AB( 1, I-K+KA )
318 AB( 1, I-K+KA ) = WORK( I-K+KA-M )*T +
319 $ WORK( N+I-K+KA-M )*AB( 1, I-K+KA )
320 RA1 = RA
321 END IF
322 END IF
323 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
324 NR = ( N-J2+KA ) / KA1
325 J1 = J2 + ( NR-1 )*KA1
326 IF( UPDATE ) THEN
327 J2T = MAX( J2, I+2*KA-K+1 )
328 ELSE
329 J2T = J2
330 END IF
331 NRT = ( N-J2T+KA ) / KA1
332 DO 90 J = J2T, J1, KA1
333 *
334 * create nonzero element a(j-ka,j+1) outside the band
335 * and store it in WORK(j-m)
336 *
337 WORK( J-M ) = WORK( J-M )*AB( 1, J+1 )
338 AB( 1, J+1 ) = WORK( N+J-M )*AB( 1, J+1 )
339 90 CONTINUE
340 *
341 * generate rotations in 1st set to annihilate elements which
342 * have been created outside the band
343 *
344 IF( NRT.GT.0 )
345 $ CALL DLARGV( NRT, AB( 1, J2T ), INCA, WORK( J2T-M ), KA1,
346 $ WORK( N+J2T-M ), KA1 )
347 IF( NR.GT.0 ) THEN
348 *
349 * apply rotations in 1st set from the right
350 *
351 DO 100 L = 1, KA - 1
352 CALL DLARTV( NR, AB( KA1-L, J2 ), INCA,
353 $ AB( KA-L, J2+1 ), INCA, WORK( N+J2-M ),
354 $ WORK( J2-M ), KA1 )
355 100 CONTINUE
356 *
357 * apply rotations in 1st set from both sides to diagonal
358 * blocks
359 *
360 CALL DLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
361 $ AB( KA, J2+1 ), INCA, WORK( N+J2-M ),
362 $ WORK( J2-M ), KA1 )
363 *
364 END IF
365 *
366 * start applying rotations in 1st set from the left
367 *
368 DO 110 L = KA - 1, KB - K + 1, -1
369 NRT = ( N-J2+L ) / KA1
370 IF( NRT.GT.0 )
371 $ CALL DLARTV( NRT, AB( L, J2+KA1-L ), INCA,
372 $ AB( L+1, J2+KA1-L ), INCA,
373 $ WORK( N+J2-M ), WORK( J2-M ), KA1 )
374 110 CONTINUE
375 *
376 IF( WANTX ) THEN
377 *
378 * post-multiply X by product of rotations in 1st set
379 *
380 DO 120 J = J2, J1, KA1
381 CALL DROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
382 $ WORK( N+J-M ), WORK( J-M ) )
383 120 CONTINUE
384 END IF
385 130 CONTINUE
386 *
387 IF( UPDATE ) THEN
388 IF( I2.LE.N .AND. KBT.GT.0 ) THEN
389 *
390 * create nonzero element a(i-kbt,i-kbt+ka+1) outside the
391 * band and store it in WORK(i-kbt)
392 *
393 WORK( I-KBT ) = -BB( KB1-KBT, I )*RA1
394 END IF
395 END IF
396 *
397 DO 170 K = KB, 1, -1
398 IF( UPDATE ) THEN
399 J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
400 ELSE
401 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
402 END IF
403 *
404 * finish applying rotations in 2nd set from the left
405 *
406 DO 140 L = KB - K, 1, -1
407 NRT = ( N-J2+KA+L ) / KA1
408 IF( NRT.GT.0 )
409 $ CALL DLARTV( NRT, AB( L, J2-L+1 ), INCA,
410 $ AB( L+1, J2-L+1 ), INCA, WORK( N+J2-KA ),
411 $ WORK( J2-KA ), KA1 )
412 140 CONTINUE
413 NR = ( N-J2+KA ) / KA1
414 J1 = J2 + ( NR-1 )*KA1
415 DO 150 J = J1, J2, -KA1
416 WORK( J ) = WORK( J-KA )
417 WORK( N+J ) = WORK( N+J-KA )
418 150 CONTINUE
419 DO 160 J = J2, J1, KA1
420 *
421 * create nonzero element a(j-ka,j+1) outside the band
422 * and store it in WORK(j)
423 *
424 WORK( J ) = WORK( J )*AB( 1, J+1 )
425 AB( 1, J+1 ) = WORK( N+J )*AB( 1, J+1 )
426 160 CONTINUE
427 IF( UPDATE ) THEN
428 IF( I-K.LT.N-KA .AND. K.LE.KBT )
429 $ WORK( I-K+KA ) = WORK( I-K )
430 END IF
431 170 CONTINUE
432 *
433 DO 210 K = KB, 1, -1
434 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
435 NR = ( N-J2+KA ) / KA1
436 J1 = J2 + ( NR-1 )*KA1
437 IF( NR.GT.0 ) THEN
438 *
439 * generate rotations in 2nd set to annihilate elements
440 * which have been created outside the band
441 *
442 CALL DLARGV( NR, AB( 1, J2 ), INCA, WORK( J2 ), KA1,
443 $ WORK( N+J2 ), KA1 )
444 *
445 * apply rotations in 2nd set from the right
446 *
447 DO 180 L = 1, KA - 1
448 CALL DLARTV( NR, AB( KA1-L, J2 ), INCA,
449 $ AB( KA-L, J2+1 ), INCA, WORK( N+J2 ),
450 $ WORK( J2 ), KA1 )
451 180 CONTINUE
452 *
453 * apply rotations in 2nd set from both sides to diagonal
454 * blocks
455 *
456 CALL DLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
457 $ AB( KA, J2+1 ), INCA, WORK( N+J2 ),
458 $ WORK( J2 ), KA1 )
459 *
460 END IF
461 *
462 * start applying rotations in 2nd set from the left
463 *
464 DO 190 L = KA - 1, KB - K + 1, -1
465 NRT = ( N-J2+L ) / KA1
466 IF( NRT.GT.0 )
467 $ CALL DLARTV( NRT, AB( L, J2+KA1-L ), INCA,
468 $ AB( L+1, J2+KA1-L ), INCA, WORK( N+J2 ),
469 $ WORK( J2 ), KA1 )
470 190 CONTINUE
471 *
472 IF( WANTX ) THEN
473 *
474 * post-multiply X by product of rotations in 2nd set
475 *
476 DO 200 J = J2, J1, KA1
477 CALL DROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
478 $ WORK( N+J ), WORK( J ) )
479 200 CONTINUE
480 END IF
481 210 CONTINUE
482 *
483 DO 230 K = 1, KB - 1
484 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
485 *
486 * finish applying rotations in 1st set from the left
487 *
488 DO 220 L = KB - K, 1, -1
489 NRT = ( N-J2+L ) / KA1
490 IF( NRT.GT.0 )
491 $ CALL DLARTV( NRT, AB( L, J2+KA1-L ), INCA,
492 $ AB( L+1, J2+KA1-L ), INCA,
493 $ WORK( N+J2-M ), WORK( J2-M ), KA1 )
494 220 CONTINUE
495 230 CONTINUE
496 *
497 IF( KB.GT.1 ) THEN
498 DO 240 J = N - 1, I - KB + 2*KA + 1, -1
499 WORK( N+J-M ) = WORK( N+J-KA-M )
500 WORK( J-M ) = WORK( J-KA-M )
501 240 CONTINUE
502 END IF
503 *
504 ELSE
505 *
506 * Transform A, working with the lower triangle
507 *
508 IF( UPDATE ) THEN
509 *
510 * Form inv(S(i))**T * A * inv(S(i))
511 *
512 BII = BB( 1, I )
513 DO 250 J = I, I1
514 AB( J-I+1, I ) = AB( J-I+1, I ) / BII
515 250 CONTINUE
516 DO 260 J = MAX( 1, I-KA ), I
517 AB( I-J+1, J ) = AB( I-J+1, J ) / BII
518 260 CONTINUE
519 DO 290 K = I - KBT, I - 1
520 DO 270 J = I - KBT, K
521 AB( K-J+1, J ) = AB( K-J+1, J ) -
522 $ BB( I-J+1, J )*AB( I-K+1, K ) -
523 $ BB( I-K+1, K )*AB( I-J+1, J ) +
524 $ AB( 1, I )*BB( I-J+1, J )*
525 $ BB( I-K+1, K )
526 270 CONTINUE
527 DO 280 J = MAX( 1, I-KA ), I - KBT - 1
528 AB( K-J+1, J ) = AB( K-J+1, J ) -
529 $ BB( I-K+1, K )*AB( I-J+1, J )
530 280 CONTINUE
531 290 CONTINUE
532 DO 310 J = I, I1
533 DO 300 K = MAX( J-KA, I-KBT ), I - 1
534 AB( J-K+1, K ) = AB( J-K+1, K ) -
535 $ BB( I-K+1, K )*AB( J-I+1, I )
536 300 CONTINUE
537 310 CONTINUE
538 *
539 IF( WANTX ) THEN
540 *
541 * post-multiply X by inv(S(i))
542 *
543 CALL DSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
544 IF( KBT.GT.0 )
545 $ CALL DGER( N-M, KBT, -ONE, X( M+1, I ), 1,
546 $ BB( KBT+1, I-KBT ), LDBB-1,
547 $ X( M+1, I-KBT ), LDX )
548 END IF
549 *
550 * store a(i1,i) in RA1 for use in next loop over K
551 *
552 RA1 = AB( I1-I+1, I )
553 END IF
554 *
555 * Generate and apply vectors of rotations to chase all the
556 * existing bulges KA positions down toward the bottom of the
557 * band
558 *
559 DO 360 K = 1, KB - 1
560 IF( UPDATE ) THEN
561 *
562 * Determine the rotations which would annihilate the bulge
563 * which has in theory just been created
564 *
565 IF( I-K+KA.LT.N .AND. I-K.GT.1 ) THEN
566 *
567 * generate rotation to annihilate a(i-k+ka+1,i)
568 *
569 CALL DLARTG( AB( KA1-K, I ), RA1, WORK( N+I-K+KA-M ),
570 $ WORK( I-K+KA-M ), RA )
571 *
572 * create nonzero element a(i-k+ka+1,i-k) outside the
573 * band and store it in WORK(i-k)
574 *
575 T = -BB( K+1, I-K )*RA1
576 WORK( I-K ) = WORK( N+I-K+KA-M )*T -
577 $ WORK( I-K+KA-M )*AB( KA1, I-K )
578 AB( KA1, I-K ) = WORK( I-K+KA-M )*T +
579 $ WORK( N+I-K+KA-M )*AB( KA1, I-K )
580 RA1 = RA
581 END IF
582 END IF
583 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
584 NR = ( N-J2+KA ) / KA1
585 J1 = J2 + ( NR-1 )*KA1
586 IF( UPDATE ) THEN
587 J2T = MAX( J2, I+2*KA-K+1 )
588 ELSE
589 J2T = J2
590 END IF
591 NRT = ( N-J2T+KA ) / KA1
592 DO 320 J = J2T, J1, KA1
593 *
594 * create nonzero element a(j+1,j-ka) outside the band
595 * and store it in WORK(j-m)
596 *
597 WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 )
598 AB( KA1, J-KA+1 ) = WORK( N+J-M )*AB( KA1, J-KA+1 )
599 320 CONTINUE
600 *
601 * generate rotations in 1st set to annihilate elements which
602 * have been created outside the band
603 *
604 IF( NRT.GT.0 )
605 $ CALL DLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ),
606 $ KA1, WORK( N+J2T-M ), KA1 )
607 IF( NR.GT.0 ) THEN
608 *
609 * apply rotations in 1st set from the left
610 *
611 DO 330 L = 1, KA - 1
612 CALL DLARTV( NR, AB( L+1, J2-L ), INCA,
613 $ AB( L+2, J2-L ), INCA, WORK( N+J2-M ),
614 $ WORK( J2-M ), KA1 )
615 330 CONTINUE
616 *
617 * apply rotations in 1st set from both sides to diagonal
618 * blocks
619 *
620 CALL DLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
621 $ INCA, WORK( N+J2-M ), WORK( J2-M ), KA1 )
622 *
623 END IF
624 *
625 * start applying rotations in 1st set from the right
626 *
627 DO 340 L = KA - 1, KB - K + 1, -1
628 NRT = ( N-J2+L ) / KA1
629 IF( NRT.GT.0 )
630 $ CALL DLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
631 $ AB( KA1-L, J2+1 ), INCA, WORK( N+J2-M ),
632 $ WORK( J2-M ), KA1 )
633 340 CONTINUE
634 *
635 IF( WANTX ) THEN
636 *
637 * post-multiply X by product of rotations in 1st set
638 *
639 DO 350 J = J2, J1, KA1
640 CALL DROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
641 $ WORK( N+J-M ), WORK( J-M ) )
642 350 CONTINUE
643 END IF
644 360 CONTINUE
645 *
646 IF( UPDATE ) THEN
647 IF( I2.LE.N .AND. KBT.GT.0 ) THEN
648 *
649 * create nonzero element a(i-kbt+ka+1,i-kbt) outside the
650 * band and store it in WORK(i-kbt)
651 *
652 WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1
653 END IF
654 END IF
655 *
656 DO 400 K = KB, 1, -1
657 IF( UPDATE ) THEN
658 J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
659 ELSE
660 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
661 END IF
662 *
663 * finish applying rotations in 2nd set from the right
664 *
665 DO 370 L = KB - K, 1, -1
666 NRT = ( N-J2+KA+L ) / KA1
667 IF( NRT.GT.0 )
668 $ CALL DLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA,
669 $ AB( KA1-L, J2-KA+1 ), INCA,
670 $ WORK( N+J2-KA ), WORK( J2-KA ), KA1 )
671 370 CONTINUE
672 NR = ( N-J2+KA ) / KA1
673 J1 = J2 + ( NR-1 )*KA1
674 DO 380 J = J1, J2, -KA1
675 WORK( J ) = WORK( J-KA )
676 WORK( N+J ) = WORK( N+J-KA )
677 380 CONTINUE
678 DO 390 J = J2, J1, KA1
679 *
680 * create nonzero element a(j+1,j-ka) outside the band
681 * and store it in WORK(j)
682 *
683 WORK( J ) = WORK( J )*AB( KA1, J-KA+1 )
684 AB( KA1, J-KA+1 ) = WORK( N+J )*AB( KA1, J-KA+1 )
685 390 CONTINUE
686 IF( UPDATE ) THEN
687 IF( I-K.LT.N-KA .AND. K.LE.KBT )
688 $ WORK( I-K+KA ) = WORK( I-K )
689 END IF
690 400 CONTINUE
691 *
692 DO 440 K = KB, 1, -1
693 J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
694 NR = ( N-J2+KA ) / KA1
695 J1 = J2 + ( NR-1 )*KA1
696 IF( NR.GT.0 ) THEN
697 *
698 * generate rotations in 2nd set to annihilate elements
699 * which have been created outside the band
700 *
701 CALL DLARGV( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1,
702 $ WORK( N+J2 ), KA1 )
703 *
704 * apply rotations in 2nd set from the left
705 *
706 DO 410 L = 1, KA - 1
707 CALL DLARTV( NR, AB( L+1, J2-L ), INCA,
708 $ AB( L+2, J2-L ), INCA, WORK( N+J2 ),
709 $ WORK( J2 ), KA1 )
710 410 CONTINUE
711 *
712 * apply rotations in 2nd set from both sides to diagonal
713 * blocks
714 *
715 CALL DLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
716 $ INCA, WORK( N+J2 ), WORK( J2 ), KA1 )
717 *
718 END IF
719 *
720 * start applying rotations in 2nd set from the right
721 *
722 DO 420 L = KA - 1, KB - K + 1, -1
723 NRT = ( N-J2+L ) / KA1
724 IF( NRT.GT.0 )
725 $ CALL DLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
726 $ AB( KA1-L, J2+1 ), INCA, WORK( N+J2 ),
727 $ WORK( J2 ), KA1 )
728 420 CONTINUE
729 *
730 IF( WANTX ) THEN
731 *
732 * post-multiply X by product of rotations in 2nd set
733 *
734 DO 430 J = J2, J1, KA1
735 CALL DROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
736 $ WORK( N+J ), WORK( J ) )
737 430 CONTINUE
738 END IF
739 440 CONTINUE
740 *
741 DO 460 K = 1, KB - 1
742 J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
743 *
744 * finish applying rotations in 1st set from the right
745 *
746 DO 450 L = KB - K, 1, -1
747 NRT = ( N-J2+L ) / KA1
748 IF( NRT.GT.0 )
749 $ CALL DLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
750 $ AB( KA1-L, J2+1 ), INCA, WORK( N+J2-M ),
751 $ WORK( J2-M ), KA1 )
752 450 CONTINUE
753 460 CONTINUE
754 *
755 IF( KB.GT.1 ) THEN
756 DO 470 J = N - 1, I - KB + 2*KA + 1, -1
757 WORK( N+J-M ) = WORK( N+J-KA-M )
758 WORK( J-M ) = WORK( J-KA-M )
759 470 CONTINUE
760 END IF
761 *
762 END IF
763 *
764 GO TO 10
765 *
766 480 CONTINUE
767 *
768 * **************************** Phase 2 *****************************
769 *
770 * The logical structure of this phase is:
771 *
772 * UPDATE = .TRUE.
773 * DO I = 1, M
774 * use S(i) to update A and create a new bulge
775 * apply rotations to push all bulges KA positions upward
776 * END DO
777 * UPDATE = .FALSE.
778 * DO I = M - KA - 1, 2, -1
779 * apply rotations to push all bulges KA positions upward
780 * END DO
781 *
782 * To avoid duplicating code, the two loops are merged.
783 *
784 UPDATE = .TRUE.
785 I = 0
786 490 CONTINUE
787 IF( UPDATE ) THEN
788 I = I + 1
789 KBT = MIN( KB, M-I )
790 I0 = I + 1
791 I1 = MAX( 1, I-KA )
792 I2 = I + KBT - KA1
793 IF( I.GT.M ) THEN
794 UPDATE = .FALSE.
795 I = I - 1
796 I0 = M + 1
797 IF( KA.EQ.0 )
798 $ RETURN
799 GO TO 490
800 END IF
801 ELSE
802 I = I - KA
803 IF( I.LT.2 )
804 $ RETURN
805 END IF
806 *
807 IF( I.LT.M-KBT ) THEN
808 NX = M
809 ELSE
810 NX = N
811 END IF
812 *
813 IF( UPPER ) THEN
814 *
815 * Transform A, working with the upper triangle
816 *
817 IF( UPDATE ) THEN
818 *
819 * Form inv(S(i))**T * A * inv(S(i))
820 *
821 BII = BB( KB1, I )
822 DO 500 J = I1, I
823 AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
824 500 CONTINUE
825 DO 510 J = I, MIN( N, I+KA )
826 AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
827 510 CONTINUE
828 DO 540 K = I + 1, I + KBT
829 DO 520 J = K, I + KBT
830 AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
831 $ BB( I-J+KB1, J )*AB( I-K+KA1, K ) -
832 $ BB( I-K+KB1, K )*AB( I-J+KA1, J ) +
833 $ AB( KA1, I )*BB( I-J+KB1, J )*
834 $ BB( I-K+KB1, K )
835 520 CONTINUE
836 DO 530 J = I + KBT + 1, MIN( N, I+KA )
837 AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
838 $ BB( I-K+KB1, K )*AB( I-J+KA1, J )
839 530 CONTINUE
840 540 CONTINUE
841 DO 560 J = I1, I
842 DO 550 K = I + 1, MIN( J+KA, I+KBT )
843 AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
844 $ BB( I-K+KB1, K )*AB( J-I+KA1, I )
845 550 CONTINUE
846 560 CONTINUE
847 *
848 IF( WANTX ) THEN
849 *
850 * post-multiply X by inv(S(i))
851 *
852 CALL DSCAL( NX, ONE / BII, X( 1, I ), 1 )
853 IF( KBT.GT.0 )
854 $ CALL DGER( NX, KBT, -ONE, X( 1, I ), 1, BB( KB, I+1 ),
855 $ LDBB-1, X( 1, I+1 ), LDX )
856 END IF
857 *
858 * store a(i1,i) in RA1 for use in next loop over K
859 *
860 RA1 = AB( I1-I+KA1, I )
861 END IF
862 *
863 * Generate and apply vectors of rotations to chase all the
864 * existing bulges KA positions up toward the top of the band
865 *
866 DO 610 K = 1, KB - 1
867 IF( UPDATE ) THEN
868 *
869 * Determine the rotations which would annihilate the bulge
870 * which has in theory just been created
871 *
872 IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
873 *
874 * generate rotation to annihilate a(i+k-ka-1,i)
875 *
876 CALL DLARTG( AB( K+1, I ), RA1, WORK( N+I+K-KA ),
877 $ WORK( I+K-KA ), RA )
878 *
879 * create nonzero element a(i+k-ka-1,i+k) outside the
880 * band and store it in WORK(m-kb+i+k)
881 *
882 T = -BB( KB1-K, I+K )*RA1
883 WORK( M-KB+I+K ) = WORK( N+I+K-KA )*T -
884 $ WORK( I+K-KA )*AB( 1, I+K )
885 AB( 1, I+K ) = WORK( I+K-KA )*T +
886 $ WORK( N+I+K-KA )*AB( 1, I+K )
887 RA1 = RA
888 END IF
889 END IF
890 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
891 NR = ( J2+KA-1 ) / KA1
892 J1 = J2 - ( NR-1 )*KA1
893 IF( UPDATE ) THEN
894 J2T = MIN( J2, I-2*KA+K-1 )
895 ELSE
896 J2T = J2
897 END IF
898 NRT = ( J2T+KA-1 ) / KA1
899 DO 570 J = J1, J2T, KA1
900 *
901 * create nonzero element a(j-1,j+ka) outside the band
902 * and store it in WORK(j)
903 *
904 WORK( J ) = WORK( J )*AB( 1, J+KA-1 )
905 AB( 1, J+KA-1 ) = WORK( N+J )*AB( 1, J+KA-1 )
906 570 CONTINUE
907 *
908 * generate rotations in 1st set to annihilate elements which
909 * have been created outside the band
910 *
911 IF( NRT.GT.0 )
912 $ CALL DLARGV( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1,
913 $ WORK( N+J1 ), KA1 )
914 IF( NR.GT.0 ) THEN
915 *
916 * apply rotations in 1st set from the left
917 *
918 DO 580 L = 1, KA - 1
919 CALL DLARTV( NR, AB( KA1-L, J1+L ), INCA,
920 $ AB( KA-L, J1+L ), INCA, WORK( N+J1 ),
921 $ WORK( J1 ), KA1 )
922 580 CONTINUE
923 *
924 * apply rotations in 1st set from both sides to diagonal
925 * blocks
926 *
927 CALL DLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
928 $ AB( KA, J1 ), INCA, WORK( N+J1 ),
929 $ WORK( J1 ), KA1 )
930 *
931 END IF
932 *
933 * start applying rotations in 1st set from the right
934 *
935 DO 590 L = KA - 1, KB - K + 1, -1
936 NRT = ( J2+L-1 ) / KA1
937 J1T = J2 - ( NRT-1 )*KA1
938 IF( NRT.GT.0 )
939 $ CALL DLARTV( NRT, AB( L, J1T ), INCA,
940 $ AB( L+1, J1T-1 ), INCA, WORK( N+J1T ),
941 $ WORK( J1T ), KA1 )
942 590 CONTINUE
943 *
944 IF( WANTX ) THEN
945 *
946 * post-multiply X by product of rotations in 1st set
947 *
948 DO 600 J = J1, J2, KA1
949 CALL DROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
950 $ WORK( N+J ), WORK( J ) )
951 600 CONTINUE
952 END IF
953 610 CONTINUE
954 *
955 IF( UPDATE ) THEN
956 IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
957 *
958 * create nonzero element a(i+kbt-ka-1,i+kbt) outside the
959 * band and store it in WORK(m-kb+i+kbt)
960 *
961 WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1
962 END IF
963 END IF
964 *
965 DO 650 K = KB, 1, -1
966 IF( UPDATE ) THEN
967 J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
968 ELSE
969 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
970 END IF
971 *
972 * finish applying rotations in 2nd set from the right
973 *
974 DO 620 L = KB - K, 1, -1
975 NRT = ( J2+KA+L-1 ) / KA1
976 J1T = J2 - ( NRT-1 )*KA1
977 IF( NRT.GT.0 )
978 $ CALL DLARTV( NRT, AB( L, J1T+KA ), INCA,
979 $ AB( L+1, J1T+KA-1 ), INCA,
980 $ WORK( N+M-KB+J1T+KA ),
981 $ WORK( M-KB+J1T+KA ), KA1 )
982 620 CONTINUE
983 NR = ( J2+KA-1 ) / KA1
984 J1 = J2 - ( NR-1 )*KA1
985 DO 630 J = J1, J2, KA1
986 WORK( M-KB+J ) = WORK( M-KB+J+KA )
987 WORK( N+M-KB+J ) = WORK( N+M-KB+J+KA )
988 630 CONTINUE
989 DO 640 J = J1, J2, KA1
990 *
991 * create nonzero element a(j-1,j+ka) outside the band
992 * and store it in WORK(m-kb+j)
993 *
994 WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 )
995 AB( 1, J+KA-1 ) = WORK( N+M-KB+J )*AB( 1, J+KA-1 )
996 640 CONTINUE
997 IF( UPDATE ) THEN
998 IF( I+K.GT.KA1 .AND. K.LE.KBT )
999 $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
1000 END IF
1001 650 CONTINUE
1002 *
1003 DO 690 K = KB, 1, -1
1004 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1005 NR = ( J2+KA-1 ) / KA1
1006 J1 = J2 - ( NR-1 )*KA1
1007 IF( NR.GT.0 ) THEN
1008 *
1009 * generate rotations in 2nd set to annihilate elements
1010 * which have been created outside the band
1011 *
1012 CALL DLARGV( NR, AB( 1, J1+KA ), INCA, WORK( M-KB+J1 ),
1013 $ KA1, WORK( N+M-KB+J1 ), KA1 )
1014 *
1015 * apply rotations in 2nd set from the left
1016 *
1017 DO 660 L = 1, KA - 1
1018 CALL DLARTV( NR, AB( KA1-L, J1+L ), INCA,
1019 $ AB( KA-L, J1+L ), INCA,
1020 $ WORK( N+M-KB+J1 ), WORK( M-KB+J1 ), KA1 )
1021 660 CONTINUE
1022 *
1023 * apply rotations in 2nd set from both sides to diagonal
1024 * blocks
1025 *
1026 CALL DLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
1027 $ AB( KA, J1 ), INCA, WORK( N+M-KB+J1 ),
1028 $ WORK( M-KB+J1 ), KA1 )
1029 *
1030 END IF
1031 *
1032 * start applying rotations in 2nd set from the right
1033 *
1034 DO 670 L = KA - 1, KB - K + 1, -1
1035 NRT = ( J2+L-1 ) / KA1
1036 J1T = J2 - ( NRT-1 )*KA1
1037 IF( NRT.GT.0 )
1038 $ CALL DLARTV( NRT, AB( L, J1T ), INCA,
1039 $ AB( L+1, J1T-1 ), INCA,
1040 $ WORK( N+M-KB+J1T ), WORK( M-KB+J1T ),
1041 $ KA1 )
1042 670 CONTINUE
1043 *
1044 IF( WANTX ) THEN
1045 *
1046 * post-multiply X by product of rotations in 2nd set
1047 *
1048 DO 680 J = J1, J2, KA1
1049 CALL DROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1050 $ WORK( N+M-KB+J ), WORK( M-KB+J ) )
1051 680 CONTINUE
1052 END IF
1053 690 CONTINUE
1054 *
1055 DO 710 K = 1, KB - 1
1056 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1057 *
1058 * finish applying rotations in 1st set from the right
1059 *
1060 DO 700 L = KB - K, 1, -1
1061 NRT = ( J2+L-1 ) / KA1
1062 J1T = J2 - ( NRT-1 )*KA1
1063 IF( NRT.GT.0 )
1064 $ CALL DLARTV( NRT, AB( L, J1T ), INCA,
1065 $ AB( L+1, J1T-1 ), INCA, WORK( N+J1T ),
1066 $ WORK( J1T ), KA1 )
1067 700 CONTINUE
1068 710 CONTINUE
1069 *
1070 IF( KB.GT.1 ) THEN
1071 DO 720 J = 2, MIN( I+KB, M ) - 2*KA - 1
1072 WORK( N+J ) = WORK( N+J+KA )
1073 WORK( J ) = WORK( J+KA )
1074 720 CONTINUE
1075 END IF
1076 *
1077 ELSE
1078 *
1079 * Transform A, working with the lower triangle
1080 *
1081 IF( UPDATE ) THEN
1082 *
1083 * Form inv(S(i))**T * A * inv(S(i))
1084 *
1085 BII = BB( 1, I )
1086 DO 730 J = I1, I
1087 AB( I-J+1, J ) = AB( I-J+1, J ) / BII
1088 730 CONTINUE
1089 DO 740 J = I, MIN( N, I+KA )
1090 AB( J-I+1, I ) = AB( J-I+1, I ) / BII
1091 740 CONTINUE
1092 DO 770 K = I + 1, I + KBT
1093 DO 750 J = K, I + KBT
1094 AB( J-K+1, K ) = AB( J-K+1, K ) -
1095 $ BB( J-I+1, I )*AB( K-I+1, I ) -
1096 $ BB( K-I+1, I )*AB( J-I+1, I ) +
1097 $ AB( 1, I )*BB( J-I+1, I )*
1098 $ BB( K-I+1, I )
1099 750 CONTINUE
1100 DO 760 J = I + KBT + 1, MIN( N, I+KA )
1101 AB( J-K+1, K ) = AB( J-K+1, K ) -
1102 $ BB( K-I+1, I )*AB( J-I+1, I )
1103 760 CONTINUE
1104 770 CONTINUE
1105 DO 790 J = I1, I
1106 DO 780 K = I + 1, MIN( J+KA, I+KBT )
1107 AB( K-J+1, J ) = AB( K-J+1, J ) -
1108 $ BB( K-I+1, I )*AB( I-J+1, J )
1109 780 CONTINUE
1110 790 CONTINUE
1111 *
1112 IF( WANTX ) THEN
1113 *
1114 * post-multiply X by inv(S(i))
1115 *
1116 CALL DSCAL( NX, ONE / BII, X( 1, I ), 1 )
1117 IF( KBT.GT.0 )
1118 $ CALL DGER( NX, KBT, -ONE, X( 1, I ), 1, BB( 2, I ), 1,
1119 $ X( 1, I+1 ), LDX )
1120 END IF
1121 *
1122 * store a(i,i1) in RA1 for use in next loop over K
1123 *
1124 RA1 = AB( I-I1+1, I1 )
1125 END IF
1126 *
1127 * Generate and apply vectors of rotations to chase all the
1128 * existing bulges KA positions up toward the top of the band
1129 *
1130 DO 840 K = 1, KB - 1
1131 IF( UPDATE ) THEN
1132 *
1133 * Determine the rotations which would annihilate the bulge
1134 * which has in theory just been created
1135 *
1136 IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
1137 *
1138 * generate rotation to annihilate a(i,i+k-ka-1)
1139 *
1140 CALL DLARTG( AB( KA1-K, I+K-KA ), RA1,
1141 $ WORK( N+I+K-KA ), WORK( I+K-KA ), RA )
1142 *
1143 * create nonzero element a(i+k,i+k-ka-1) outside the
1144 * band and store it in WORK(m-kb+i+k)
1145 *
1146 T = -BB( K+1, I )*RA1
1147 WORK( M-KB+I+K ) = WORK( N+I+K-KA )*T -
1148 $ WORK( I+K-KA )*AB( KA1, I+K-KA )
1149 AB( KA1, I+K-KA ) = WORK( I+K-KA )*T +
1150 $ WORK( N+I+K-KA )*AB( KA1, I+K-KA )
1151 RA1 = RA
1152 END IF
1153 END IF
1154 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1155 NR = ( J2+KA-1 ) / KA1
1156 J1 = J2 - ( NR-1 )*KA1
1157 IF( UPDATE ) THEN
1158 J2T = MIN( J2, I-2*KA+K-1 )
1159 ELSE
1160 J2T = J2
1161 END IF
1162 NRT = ( J2T+KA-1 ) / KA1
1163 DO 800 J = J1, J2T, KA1
1164 *
1165 * create nonzero element a(j+ka,j-1) outside the band
1166 * and store it in WORK(j)
1167 *
1168 WORK( J ) = WORK( J )*AB( KA1, J-1 )
1169 AB( KA1, J-1 ) = WORK( N+J )*AB( KA1, J-1 )
1170 800 CONTINUE
1171 *
1172 * generate rotations in 1st set to annihilate elements which
1173 * have been created outside the band
1174 *
1175 IF( NRT.GT.0 )
1176 $ CALL DLARGV( NRT, AB( KA1, J1 ), INCA, WORK( J1 ), KA1,
1177 $ WORK( N+J1 ), KA1 )
1178 IF( NR.GT.0 ) THEN
1179 *
1180 * apply rotations in 1st set from the right
1181 *
1182 DO 810 L = 1, KA - 1
1183 CALL DLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
1184 $ INCA, WORK( N+J1 ), WORK( J1 ), KA1 )
1185 810 CONTINUE
1186 *
1187 * apply rotations in 1st set from both sides to diagonal
1188 * blocks
1189 *
1190 CALL DLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
1191 $ AB( 2, J1-1 ), INCA, WORK( N+J1 ),
1192 $ WORK( J1 ), KA1 )
1193 *
1194 END IF
1195 *
1196 * start applying rotations in 1st set from the left
1197 *
1198 DO 820 L = KA - 1, KB - K + 1, -1
1199 NRT = ( J2+L-1 ) / KA1
1200 J1T = J2 - ( NRT-1 )*KA1
1201 IF( NRT.GT.0 )
1202 $ CALL DLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1203 $ AB( KA1-L, J1T-KA1+L ), INCA,
1204 $ WORK( N+J1T ), WORK( J1T ), KA1 )
1205 820 CONTINUE
1206 *
1207 IF( WANTX ) THEN
1208 *
1209 * post-multiply X by product of rotations in 1st set
1210 *
1211 DO 830 J = J1, J2, KA1
1212 CALL DROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1213 $ WORK( N+J ), WORK( J ) )
1214 830 CONTINUE
1215 END IF
1216 840 CONTINUE
1217 *
1218 IF( UPDATE ) THEN
1219 IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
1220 *
1221 * create nonzero element a(i+kbt,i+kbt-ka-1) outside the
1222 * band and store it in WORK(m-kb+i+kbt)
1223 *
1224 WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1
1225 END IF
1226 END IF
1227 *
1228 DO 880 K = KB, 1, -1
1229 IF( UPDATE ) THEN
1230 J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
1231 ELSE
1232 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1233 END IF
1234 *
1235 * finish applying rotations in 2nd set from the left
1236 *
1237 DO 850 L = KB - K, 1, -1
1238 NRT = ( J2+KA+L-1 ) / KA1
1239 J1T = J2 - ( NRT-1 )*KA1
1240 IF( NRT.GT.0 )
1241 $ CALL DLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA,
1242 $ AB( KA1-L, J1T+L-1 ), INCA,
1243 $ WORK( N+M-KB+J1T+KA ),
1244 $ WORK( M-KB+J1T+KA ), KA1 )
1245 850 CONTINUE
1246 NR = ( J2+KA-1 ) / KA1
1247 J1 = J2 - ( NR-1 )*KA1
1248 DO 860 J = J1, J2, KA1
1249 WORK( M-KB+J ) = WORK( M-KB+J+KA )
1250 WORK( N+M-KB+J ) = WORK( N+M-KB+J+KA )
1251 860 CONTINUE
1252 DO 870 J = J1, J2, KA1
1253 *
1254 * create nonzero element a(j+ka,j-1) outside the band
1255 * and store it in WORK(m-kb+j)
1256 *
1257 WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 )
1258 AB( KA1, J-1 ) = WORK( N+M-KB+J )*AB( KA1, J-1 )
1259 870 CONTINUE
1260 IF( UPDATE ) THEN
1261 IF( I+K.GT.KA1 .AND. K.LE.KBT )
1262 $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
1263 END IF
1264 880 CONTINUE
1265 *
1266 DO 920 K = KB, 1, -1
1267 J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1268 NR = ( J2+KA-1 ) / KA1
1269 J1 = J2 - ( NR-1 )*KA1
1270 IF( NR.GT.0 ) THEN
1271 *
1272 * generate rotations in 2nd set to annihilate elements
1273 * which have been created outside the band
1274 *
1275 CALL DLARGV( NR, AB( KA1, J1 ), INCA, WORK( M-KB+J1 ),
1276 $ KA1, WORK( N+M-KB+J1 ), KA1 )
1277 *
1278 * apply rotations in 2nd set from the right
1279 *
1280 DO 890 L = 1, KA - 1
1281 CALL DLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
1282 $ INCA, WORK( N+M-KB+J1 ), WORK( M-KB+J1 ),
1283 $ KA1 )
1284 890 CONTINUE
1285 *
1286 * apply rotations in 2nd set from both sides to diagonal
1287 * blocks
1288 *
1289 CALL DLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
1290 $ AB( 2, J1-1 ), INCA, WORK( N+M-KB+J1 ),
1291 $ WORK( M-KB+J1 ), KA1 )
1292 *
1293 END IF
1294 *
1295 * start applying rotations in 2nd set from the left
1296 *
1297 DO 900 L = KA - 1, KB - K + 1, -1
1298 NRT = ( J2+L-1 ) / KA1
1299 J1T = J2 - ( NRT-1 )*KA1
1300 IF( NRT.GT.0 )
1301 $ CALL DLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1302 $ AB( KA1-L, J1T-KA1+L ), INCA,
1303 $ WORK( N+M-KB+J1T ), WORK( M-KB+J1T ),
1304 $ KA1 )
1305 900 CONTINUE
1306 *
1307 IF( WANTX ) THEN
1308 *
1309 * post-multiply X by product of rotations in 2nd set
1310 *
1311 DO 910 J = J1, J2, KA1
1312 CALL DROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1313 $ WORK( N+M-KB+J ), WORK( M-KB+J ) )
1314 910 CONTINUE
1315 END IF
1316 920 CONTINUE
1317 *
1318 DO 940 K = 1, KB - 1
1319 J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1320 *
1321 * finish applying rotations in 1st set from the left
1322 *
1323 DO 930 L = KB - K, 1, -1
1324 NRT = ( J2+L-1 ) / KA1
1325 J1T = J2 - ( NRT-1 )*KA1
1326 IF( NRT.GT.0 )
1327 $ CALL DLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1328 $ AB( KA1-L, J1T-KA1+L ), INCA,
1329 $ WORK( N+J1T ), WORK( J1T ), KA1 )
1330 930 CONTINUE
1331 940 CONTINUE
1332 *
1333 IF( KB.GT.1 ) THEN
1334 DO 950 J = 2, MIN( I+KB, M ) - 2*KA - 1
1335 WORK( N+J ) = WORK( N+J+KA )
1336 WORK( J ) = WORK( J+KA )
1337 950 CONTINUE
1338 END IF
1339 *
1340 END IF
1341 *
1342 GO TO 490
1343 *
1344 * End of DSBGST
1345 *
1346 END