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