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