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