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