LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
chbgst.f
Go to the documentation of this file.
1 *> \brief \b CHBGST
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CHBGST + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chbgst.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chbgst.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chbgst.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CHBGST( 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 * REAL RWORK( * )
30 * COMPLEX AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
31 * $ X( LDX, * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> CHBGST 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 CPBSTF, 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 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 array, dimension (LDBB,N)
110 *> The banded factor S from the split Cholesky factorization of
111 *> B, as returned by CPBSTF, 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 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 array, dimension (N)
138 *> \endverbatim
139 *>
140 *> \param[out] RWORK
141 *> \verbatim
142 *> RWORK is REAL 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 complexOTHERcomputational
163 *
164 * =====================================================================
165  SUBROUTINE chbgst( 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  REAL rwork( * )
179  COMPLEX ab( ldab, * ), bb( ldbb, * ), work( * ),
180  $ x( ldx, * )
181 * ..
182 *
183 * =====================================================================
184 *
185 * .. Parameters ..
186  COMPLEX czero, cone
187  REAL one
188  parameter( czero = ( 0.0e+0, 0.0e+0 ),
189  $ cone = ( 1.0e+0, 0.0e+0 ), one = 1.0e+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  REAL bii
196  COMPLEX ra, ra1, t
197 * ..
198 * .. External Functions ..
199  LOGICAL lsame
200  EXTERNAL lsame
201 * ..
202 * .. External Subroutines ..
203  EXTERNAL cgerc, cgeru, clacgv, clar2v, clargv, clartg,
204  $ clartv, claset, crot, csscal, xerbla
205 * ..
206 * .. Intrinsic Functions ..
207  INTRINSIC conjg, max, min, real
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( 'CHBGST', -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 claset( '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 CPBSTF. 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 = REAL( BB( KB1, I ) )
349  ab( ka1, i ) = ( REAL( 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  $ conjg( ab( k-i+ka1, i ) ) -
361  $ conjg( bb( k-i+kb1, i ) )*
362  $ ab( j-i+ka1, i ) +
363  $ REAL( AB( KA1, I ) )*
364  $ bb( j-i+kb1, i )*
365  $ conjg( 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  $ conjg( 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 csscal( n-m, one / bii, x( m+1, i ), 1 )
385  IF( kbt.GT.0 )
386  $ CALL cgerc( 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 clartg( 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  $ conjg( 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 clargv( 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 clartv( 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 clar2v( 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 clacgv( 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 clartv( 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 crot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
485  $ rwork( j-m ), conjg( 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 clartv( 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 clargv( 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 clartv( 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 clar2v( nr, ab( ka1, j2 ), ab( ka1, j2+1 ),
560  $ ab( ka, j2+1 ), inca, rwork( j2 ),
561  $ work( j2 ), ka1 )
562 *
563  CALL clacgv( 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 clartv( 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 crot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
582  $ rwork( j ), conjg( 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 clartv( 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 = REAL( BB( 1, I ) )
617  ab( 1, i ) = ( REAL( 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 )*conjg( ab( i-k+1,
628  $ k ) ) - conjg( bb( i-k+1, k ) )*
629  $ ab( i-j+1, j ) + REAL( AB( 1, I ) )*
630  $ bb( i-j+1, j )*conjg( 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  $ conjg( 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 csscal( n-m, one / bii, x( m+1, i ), 1 )
651  IF( kbt.GT.0 )
652  $ CALL cgeru( 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 clartg( 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  $ conjg( work( i-k+ka-m ) )*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 clargv( nrt, ab( ka1, j2t-ka ), inca, work( j2t-m ),
713  $ ka1, rwork( j2t-m ), ka1 )
714  IF( nr.GT.0 ) THEN
715 *
716 * apply rotations in 1st set from the left
717 *
718  DO 330 l = 1, ka - 1
719  CALL clartv( nr, ab( l+1, j2-l ), inca,
720  $ ab( l+2, j2-l ), inca, rwork( j2-m ),
721  $ work( j2-m ), ka1 )
722  330 continue
723 *
724 * apply rotations in 1st set from both sides to diagonal
725 * blocks
726 *
727  CALL clar2v( nr, ab( 1, j2 ), ab( 1, j2+1 ), ab( 2, j2 ),
728  $ inca, rwork( j2-m ), work( j2-m ), ka1 )
729 *
730  CALL clacgv( nr, work( j2-m ), ka1 )
731  END IF
732 *
733 * start applying rotations in 1st set from the right
734 *
735  DO 340 l = ka - 1, kb - k + 1, -1
736  nrt = ( n-j2+l ) / ka1
737  IF( nrt.GT.0 )
738  $ CALL clartv( nrt, ab( ka1-l+1, j2 ), inca,
739  $ ab( ka1-l, j2+1 ), inca, rwork( j2-m ),
740  $ work( j2-m ), ka1 )
741  340 continue
742 *
743  IF( wantx ) THEN
744 *
745 * post-multiply X by product of rotations in 1st set
746 *
747  DO 350 j = j2, j1, ka1
748  CALL crot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
749  $ rwork( j-m ), work( j-m ) )
750  350 continue
751  END IF
752  360 continue
753 *
754  IF( update ) THEN
755  IF( i2.LE.n .AND. kbt.GT.0 ) THEN
756 *
757 * create nonzero element a(i-kbt+ka+1,i-kbt) outside the
758 * band and store it in WORK(i-kbt)
759 *
760  work( i-kbt ) = -bb( kbt+1, i-kbt )*ra1
761  END IF
762  END IF
763 *
764  DO 400 k = kb, 1, -1
765  IF( update ) THEN
766  j2 = i - k - 1 + max( 2, k-i0+1 )*ka1
767  ELSE
768  j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
769  END IF
770 *
771 * finish applying rotations in 2nd set from the right
772 *
773  DO 370 l = kb - k, 1, -1
774  nrt = ( n-j2+ka+l ) / ka1
775  IF( nrt.GT.0 )
776  $ CALL clartv( nrt, ab( ka1-l+1, j2-ka ), inca,
777  $ ab( ka1-l, j2-ka+1 ), inca,
778  $ rwork( j2-ka ), work( j2-ka ), ka1 )
779  370 continue
780  nr = ( n-j2+ka ) / ka1
781  j1 = j2 + ( nr-1 )*ka1
782  DO 380 j = j1, j2, -ka1
783  work( j ) = work( j-ka )
784  rwork( j ) = rwork( j-ka )
785  380 continue
786  DO 390 j = j2, j1, ka1
787 *
788 * create nonzero element a(j+1,j-ka) outside the band
789 * and store it in WORK(j)
790 *
791  work( j ) = work( j )*ab( ka1, j-ka+1 )
792  ab( ka1, j-ka+1 ) = rwork( j )*ab( ka1, j-ka+1 )
793  390 continue
794  IF( update ) THEN
795  IF( i-k.LT.n-ka .AND. k.LE.kbt )
796  $ work( i-k+ka ) = work( i-k )
797  END IF
798  400 continue
799 *
800  DO 440 k = kb, 1, -1
801  j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
802  nr = ( n-j2+ka ) / ka1
803  j1 = j2 + ( nr-1 )*ka1
804  IF( nr.GT.0 ) THEN
805 *
806 * generate rotations in 2nd set to annihilate elements
807 * which have been created outside the band
808 *
809  CALL clargv( nr, ab( ka1, j2-ka ), inca, work( j2 ), ka1,
810  $ rwork( j2 ), ka1 )
811 *
812 * apply rotations in 2nd set from the left
813 *
814  DO 410 l = 1, ka - 1
815  CALL clartv( nr, ab( l+1, j2-l ), inca,
816  $ ab( l+2, j2-l ), inca, rwork( j2 ),
817  $ work( j2 ), ka1 )
818  410 continue
819 *
820 * apply rotations in 2nd set from both sides to diagonal
821 * blocks
822 *
823  CALL clar2v( nr, ab( 1, j2 ), ab( 1, j2+1 ), ab( 2, j2 ),
824  $ inca, rwork( j2 ), work( j2 ), ka1 )
825 *
826  CALL clacgv( nr, work( j2 ), ka1 )
827  END IF
828 *
829 * start applying rotations in 2nd set from the right
830 *
831  DO 420 l = ka - 1, kb - k + 1, -1
832  nrt = ( n-j2+l ) / ka1
833  IF( nrt.GT.0 )
834  $ CALL clartv( nrt, ab( ka1-l+1, j2 ), inca,
835  $ ab( ka1-l, j2+1 ), inca, rwork( j2 ),
836  $ work( j2 ), ka1 )
837  420 continue
838 *
839  IF( wantx ) THEN
840 *
841 * post-multiply X by product of rotations in 2nd set
842 *
843  DO 430 j = j2, j1, ka1
844  CALL crot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
845  $ rwork( j ), work( j ) )
846  430 continue
847  END IF
848  440 continue
849 *
850  DO 460 k = 1, kb - 1
851  j2 = i - k - 1 + max( 1, k-i0+2 )*ka1
852 *
853 * finish applying rotations in 1st set from the right
854 *
855  DO 450 l = kb - k, 1, -1
856  nrt = ( n-j2+l ) / ka1
857  IF( nrt.GT.0 )
858  $ CALL clartv( nrt, ab( ka1-l+1, j2 ), inca,
859  $ ab( ka1-l, j2+1 ), inca, rwork( j2-m ),
860  $ work( j2-m ), ka1 )
861  450 continue
862  460 continue
863 *
864  IF( kb.GT.1 ) THEN
865  DO 470 j = n - 1, j2 + ka, -1
866  rwork( j-m ) = rwork( j-ka-m )
867  work( j-m ) = work( j-ka-m )
868  470 continue
869  END IF
870 *
871  END IF
872 *
873  go to 10
874 *
875  480 continue
876 *
877 * **************************** Phase 2 *****************************
878 *
879 * The logical structure of this phase is:
880 *
881 * UPDATE = .TRUE.
882 * DO I = 1, M
883 * use S(i) to update A and create a new bulge
884 * apply rotations to push all bulges KA positions upward
885 * END DO
886 * UPDATE = .FALSE.
887 * DO I = M - KA - 1, 2, -1
888 * apply rotations to push all bulges KA positions upward
889 * END DO
890 *
891 * To avoid duplicating code, the two loops are merged.
892 *
893  update = .true.
894  i = 0
895  490 continue
896  IF( update ) THEN
897  i = i + 1
898  kbt = min( kb, m-i )
899  i0 = i + 1
900  i1 = max( 1, i-ka )
901  i2 = i + kbt - ka1
902  IF( i.GT.m ) THEN
903  update = .false.
904  i = i - 1
905  i0 = m + 1
906  IF( ka.EQ.0 )
907  $ return
908  go to 490
909  END IF
910  ELSE
911  i = i - ka
912  IF( i.LT.2 )
913  $ return
914  END IF
915 *
916  IF( i.LT.m-kbt ) THEN
917  nx = m
918  ELSE
919  nx = n
920  END IF
921 *
922  IF( upper ) THEN
923 *
924 * Transform A, working with the upper triangle
925 *
926  IF( update ) THEN
927 *
928 * Form inv(S(i))**H * A * inv(S(i))
929 *
930  bii = REAL( BB( KB1, I ) )
931  ab( ka1, i ) = ( REAL( AB( KA1, I ) ) / bii ) / bii
932  DO 500 j = i1, i - 1
933  ab( j-i+ka1, i ) = ab( j-i+ka1, i ) / bii
934  500 continue
935  DO 510 j = i + 1, min( n, i+ka )
936  ab( i-j+ka1, j ) = ab( i-j+ka1, j ) / bii
937  510 continue
938  DO 540 k = i + 1, i + kbt
939  DO 520 j = k, i + kbt
940  ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -
941  $ bb( i-j+kb1, j )*
942  $ conjg( ab( i-k+ka1, k ) ) -
943  $ conjg( bb( i-k+kb1, k ) )*
944  $ ab( i-j+ka1, j ) +
945  $ REAL( AB( KA1, I ) )*
946  $ bb( i-j+kb1, j )*
947  $ conjg( bb( i-k+kb1, k ) )
948  520 continue
949  DO 530 j = i + kbt + 1, min( n, i+ka )
950  ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -
951  $ conjg( bb( i-k+kb1, k ) )*
952  $ ab( i-j+ka1, j )
953  530 continue
954  540 continue
955  DO 560 j = i1, i
956  DO 550 k = i + 1, min( j+ka, i+kbt )
957  ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -
958  $ bb( i-k+kb1, k )*ab( j-i+ka1, i )
959  550 continue
960  560 continue
961 *
962  IF( wantx ) THEN
963 *
964 * post-multiply X by inv(S(i))
965 *
966  CALL csscal( nx, one / bii, x( 1, i ), 1 )
967  IF( kbt.GT.0 )
968  $ CALL cgeru( nx, kbt, -cone, x( 1, i ), 1,
969  $ bb( kb, i+1 ), ldbb-1, x( 1, i+1 ), ldx )
970  END IF
971 *
972 * store a(i1,i) in RA1 for use in next loop over K
973 *
974  ra1 = ab( i1-i+ka1, i )
975  END IF
976 *
977 * Generate and apply vectors of rotations to chase all the
978 * existing bulges KA positions up toward the top of the band
979 *
980  DO 610 k = 1, kb - 1
981  IF( update ) THEN
982 *
983 * Determine the rotations which would annihilate the bulge
984 * which has in theory just been created
985 *
986  IF( i+k-ka1.GT.0 .AND. i+k.LT.m ) THEN
987 *
988 * generate rotation to annihilate a(i+k-ka-1,i)
989 *
990  CALL clartg( ab( k+1, i ), ra1, rwork( i+k-ka ),
991  $ work( i+k-ka ), ra )
992 *
993 * create nonzero element a(i+k-ka-1,i+k) outside the
994 * band and store it in WORK(m-kb+i+k)
995 *
996  t = -bb( kb1-k, i+k )*ra1
997  work( m-kb+i+k ) = rwork( i+k-ka )*t -
998  $ conjg( work( i+k-ka ) )*
999  $ ab( 1, i+k )
1000  ab( 1, i+k ) = work( i+k-ka )*t +
1001  $ rwork( i+k-ka )*ab( 1, i+k )
1002  ra1 = ra
1003  END IF
1004  END IF
1005  j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1006  nr = ( j2+ka-1 ) / ka1
1007  j1 = j2 - ( nr-1 )*ka1
1008  IF( update ) THEN
1009  j2t = min( j2, i-2*ka+k-1 )
1010  ELSE
1011  j2t = j2
1012  END IF
1013  nrt = ( j2t+ka-1 ) / ka1
1014  DO 570 j = j1, j2t, ka1
1015 *
1016 * create nonzero element a(j-1,j+ka) outside the band
1017 * and store it in WORK(j)
1018 *
1019  work( j ) = work( j )*ab( 1, j+ka-1 )
1020  ab( 1, j+ka-1 ) = rwork( j )*ab( 1, j+ka-1 )
1021  570 continue
1022 *
1023 * generate rotations in 1st set to annihilate elements which
1024 * have been created outside the band
1025 *
1026  IF( nrt.GT.0 )
1027  $ CALL clargv( nrt, ab( 1, j1+ka ), inca, work( j1 ), ka1,
1028  $ rwork( j1 ), ka1 )
1029  IF( nr.GT.0 ) THEN
1030 *
1031 * apply rotations in 1st set from the left
1032 *
1033  DO 580 l = 1, ka - 1
1034  CALL clartv( nr, ab( ka1-l, j1+l ), inca,
1035  $ ab( ka-l, j1+l ), inca, rwork( j1 ),
1036  $ work( j1 ), ka1 )
1037  580 continue
1038 *
1039 * apply rotations in 1st set from both sides to diagonal
1040 * blocks
1041 *
1042  CALL clar2v( nr, ab( ka1, j1 ), ab( ka1, j1-1 ),
1043  $ ab( ka, j1 ), inca, rwork( j1 ), work( j1 ),
1044  $ ka1 )
1045 *
1046  CALL clacgv( nr, work( j1 ), ka1 )
1047  END IF
1048 *
1049 * start applying rotations in 1st set from the right
1050 *
1051  DO 590 l = ka - 1, kb - k + 1, -1
1052  nrt = ( j2+l-1 ) / ka1
1053  j1t = j2 - ( nrt-1 )*ka1
1054  IF( nrt.GT.0 )
1055  $ CALL clartv( nrt, ab( l, j1t ), inca,
1056  $ ab( l+1, j1t-1 ), inca, rwork( j1t ),
1057  $ work( j1t ), ka1 )
1058  590 continue
1059 *
1060  IF( wantx ) THEN
1061 *
1062 * post-multiply X by product of rotations in 1st set
1063 *
1064  DO 600 j = j1, j2, ka1
1065  CALL crot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1066  $ rwork( j ), work( j ) )
1067  600 continue
1068  END IF
1069  610 continue
1070 *
1071  IF( update ) THEN
1072  IF( i2.GT.0 .AND. kbt.GT.0 ) THEN
1073 *
1074 * create nonzero element a(i+kbt-ka-1,i+kbt) outside the
1075 * band and store it in WORK(m-kb+i+kbt)
1076 *
1077  work( m-kb+i+kbt ) = -bb( kb1-kbt, i+kbt )*ra1
1078  END IF
1079  END IF
1080 *
1081  DO 650 k = kb, 1, -1
1082  IF( update ) THEN
1083  j2 = i + k + 1 - max( 2, k+i0-m )*ka1
1084  ELSE
1085  j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1086  END IF
1087 *
1088 * finish applying rotations in 2nd set from the right
1089 *
1090  DO 620 l = kb - k, 1, -1
1091  nrt = ( j2+ka+l-1 ) / ka1
1092  j1t = j2 - ( nrt-1 )*ka1
1093  IF( nrt.GT.0 )
1094  $ CALL clartv( nrt, ab( l, j1t+ka ), inca,
1095  $ ab( l+1, j1t+ka-1 ), inca,
1096  $ rwork( m-kb+j1t+ka ),
1097  $ work( m-kb+j1t+ka ), ka1 )
1098  620 continue
1099  nr = ( j2+ka-1 ) / ka1
1100  j1 = j2 - ( nr-1 )*ka1
1101  DO 630 j = j1, j2, ka1
1102  work( m-kb+j ) = work( m-kb+j+ka )
1103  rwork( m-kb+j ) = rwork( m-kb+j+ka )
1104  630 continue
1105  DO 640 j = j1, j2, ka1
1106 *
1107 * create nonzero element a(j-1,j+ka) outside the band
1108 * and store it in WORK(m-kb+j)
1109 *
1110  work( m-kb+j ) = work( m-kb+j )*ab( 1, j+ka-1 )
1111  ab( 1, j+ka-1 ) = rwork( m-kb+j )*ab( 1, j+ka-1 )
1112  640 continue
1113  IF( update ) THEN
1114  IF( i+k.GT.ka1 .AND. k.LE.kbt )
1115  $ work( m-kb+i+k-ka ) = work( m-kb+i+k )
1116  END IF
1117  650 continue
1118 *
1119  DO 690 k = kb, 1, -1
1120  j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1121  nr = ( j2+ka-1 ) / ka1
1122  j1 = j2 - ( nr-1 )*ka1
1123  IF( nr.GT.0 ) THEN
1124 *
1125 * generate rotations in 2nd set to annihilate elements
1126 * which have been created outside the band
1127 *
1128  CALL clargv( nr, ab( 1, j1+ka ), inca, work( m-kb+j1 ),
1129  $ ka1, rwork( m-kb+j1 ), ka1 )
1130 *
1131 * apply rotations in 2nd set from the left
1132 *
1133  DO 660 l = 1, ka - 1
1134  CALL clartv( nr, ab( ka1-l, j1+l ), inca,
1135  $ ab( ka-l, j1+l ), inca, rwork( m-kb+j1 ),
1136  $ work( m-kb+j1 ), ka1 )
1137  660 continue
1138 *
1139 * apply rotations in 2nd set from both sides to diagonal
1140 * blocks
1141 *
1142  CALL clar2v( nr, ab( ka1, j1 ), ab( ka1, j1-1 ),
1143  $ ab( ka, j1 ), inca, rwork( m-kb+j1 ),
1144  $ work( m-kb+j1 ), ka1 )
1145 *
1146  CALL clacgv( nr, work( m-kb+j1 ), ka1 )
1147  END IF
1148 *
1149 * start applying rotations in 2nd set from the right
1150 *
1151  DO 670 l = ka - 1, kb - k + 1, -1
1152  nrt = ( j2+l-1 ) / ka1
1153  j1t = j2 - ( nrt-1 )*ka1
1154  IF( nrt.GT.0 )
1155  $ CALL clartv( nrt, ab( l, j1t ), inca,
1156  $ ab( l+1, j1t-1 ), inca,
1157  $ rwork( m-kb+j1t ), work( m-kb+j1t ),
1158  $ ka1 )
1159  670 continue
1160 *
1161  IF( wantx ) THEN
1162 *
1163 * post-multiply X by product of rotations in 2nd set
1164 *
1165  DO 680 j = j1, j2, ka1
1166  CALL crot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1167  $ rwork( m-kb+j ), work( m-kb+j ) )
1168  680 continue
1169  END IF
1170  690 continue
1171 *
1172  DO 710 k = 1, kb - 1
1173  j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1174 *
1175 * finish applying rotations in 1st set from the right
1176 *
1177  DO 700 l = kb - k, 1, -1
1178  nrt = ( j2+l-1 ) / ka1
1179  j1t = j2 - ( nrt-1 )*ka1
1180  IF( nrt.GT.0 )
1181  $ CALL clartv( nrt, ab( l, j1t ), inca,
1182  $ ab( l+1, j1t-1 ), inca, rwork( j1t ),
1183  $ work( j1t ), ka1 )
1184  700 continue
1185  710 continue
1186 *
1187  IF( kb.GT.1 ) THEN
1188  DO 720 j = 2, i2 - ka
1189  rwork( j ) = rwork( j+ka )
1190  work( j ) = work( j+ka )
1191  720 continue
1192  END IF
1193 *
1194  ELSE
1195 *
1196 * Transform A, working with the lower triangle
1197 *
1198  IF( update ) THEN
1199 *
1200 * Form inv(S(i))**H * A * inv(S(i))
1201 *
1202  bii = REAL( BB( 1, I ) )
1203  ab( 1, i ) = ( REAL( AB( 1, I ) ) / bii ) / bii
1204  DO 730 j = i1, i - 1
1205  ab( i-j+1, j ) = ab( i-j+1, j ) / bii
1206  730 continue
1207  DO 740 j = i + 1, min( n, i+ka )
1208  ab( j-i+1, i ) = ab( j-i+1, i ) / bii
1209  740 continue
1210  DO 770 k = i + 1, i + kbt
1211  DO 750 j = k, i + kbt
1212  ab( j-k+1, k ) = ab( j-k+1, k ) -
1213  $ bb( j-i+1, i )*conjg( ab( k-i+1,
1214  $ i ) ) - conjg( bb( k-i+1, i ) )*
1215  $ ab( j-i+1, i ) + REAL( AB( 1, I ) )*
1216  $ bb( j-i+1, i )*conjg( bb( k-i+1,
1217  $ i ) )
1218  750 continue
1219  DO 760 j = i + kbt + 1, min( n, i+ka )
1220  ab( j-k+1, k ) = ab( j-k+1, k ) -
1221  $ conjg( bb( k-i+1, i ) )*
1222  $ ab( j-i+1, i )
1223  760 continue
1224  770 continue
1225  DO 790 j = i1, i
1226  DO 780 k = i + 1, min( j+ka, i+kbt )
1227  ab( k-j+1, j ) = ab( k-j+1, j ) -
1228  $ bb( k-i+1, i )*ab( i-j+1, j )
1229  780 continue
1230  790 continue
1231 *
1232  IF( wantx ) THEN
1233 *
1234 * post-multiply X by inv(S(i))
1235 *
1236  CALL csscal( nx, one / bii, x( 1, i ), 1 )
1237  IF( kbt.GT.0 )
1238  $ CALL cgerc( nx, kbt, -cone, x( 1, i ), 1, bb( 2, i ),
1239  $ 1, x( 1, i+1 ), ldx )
1240  END IF
1241 *
1242 * store a(i,i1) in RA1 for use in next loop over K
1243 *
1244  ra1 = ab( i-i1+1, i1 )
1245  END IF
1246 *
1247 * Generate and apply vectors of rotations to chase all the
1248 * existing bulges KA positions up toward the top of the band
1249 *
1250  DO 840 k = 1, kb - 1
1251  IF( update ) THEN
1252 *
1253 * Determine the rotations which would annihilate the bulge
1254 * which has in theory just been created
1255 *
1256  IF( i+k-ka1.GT.0 .AND. i+k.LT.m ) THEN
1257 *
1258 * generate rotation to annihilate a(i,i+k-ka-1)
1259 *
1260  CALL clartg( ab( ka1-k, i+k-ka ), ra1,
1261  $ rwork( i+k-ka ), work( i+k-ka ), ra )
1262 *
1263 * create nonzero element a(i+k,i+k-ka-1) outside the
1264 * band and store it in WORK(m-kb+i+k)
1265 *
1266  t = -bb( k+1, i )*ra1
1267  work( m-kb+i+k ) = rwork( i+k-ka )*t -
1268  $ conjg( work( i+k-ka ) )*
1269  $ ab( ka1, i+k-ka )
1270  ab( ka1, i+k-ka ) = work( i+k-ka )*t +
1271  $ rwork( i+k-ka )*ab( ka1, i+k-ka )
1272  ra1 = ra
1273  END IF
1274  END IF
1275  j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1276  nr = ( j2+ka-1 ) / ka1
1277  j1 = j2 - ( nr-1 )*ka1
1278  IF( update ) THEN
1279  j2t = min( j2, i-2*ka+k-1 )
1280  ELSE
1281  j2t = j2
1282  END IF
1283  nrt = ( j2t+ka-1 ) / ka1
1284  DO 800 j = j1, j2t, ka1
1285 *
1286 * create nonzero element a(j+ka,j-1) outside the band
1287 * and store it in WORK(j)
1288 *
1289  work( j ) = work( j )*ab( ka1, j-1 )
1290  ab( ka1, j-1 ) = rwork( j )*ab( ka1, j-1 )
1291  800 continue
1292 *
1293 * generate rotations in 1st set to annihilate elements which
1294 * have been created outside the band
1295 *
1296  IF( nrt.GT.0 )
1297  $ CALL clargv( nrt, ab( ka1, j1 ), inca, work( j1 ), ka1,
1298  $ rwork( j1 ), ka1 )
1299  IF( nr.GT.0 ) THEN
1300 *
1301 * apply rotations in 1st set from the right
1302 *
1303  DO 810 l = 1, ka - 1
1304  CALL clartv( nr, ab( l+1, j1 ), inca, ab( l+2, j1-1 ),
1305  $ inca, rwork( j1 ), work( j1 ), ka1 )
1306  810 continue
1307 *
1308 * apply rotations in 1st set from both sides to diagonal
1309 * blocks
1310 *
1311  CALL clar2v( nr, ab( 1, j1 ), ab( 1, j1-1 ),
1312  $ ab( 2, j1-1 ), inca, rwork( j1 ),
1313  $ work( j1 ), ka1 )
1314 *
1315  CALL clacgv( nr, work( j1 ), ka1 )
1316  END IF
1317 *
1318 * start applying rotations in 1st set from the left
1319 *
1320  DO 820 l = ka - 1, kb - k + 1, -1
1321  nrt = ( j2+l-1 ) / ka1
1322  j1t = j2 - ( nrt-1 )*ka1
1323  IF( nrt.GT.0 )
1324  $ CALL clartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,
1325  $ ab( ka1-l, j1t-ka1+l ), inca,
1326  $ rwork( j1t ), work( j1t ), ka1 )
1327  820 continue
1328 *
1329  IF( wantx ) THEN
1330 *
1331 * post-multiply X by product of rotations in 1st set
1332 *
1333  DO 830 j = j1, j2, ka1
1334  CALL crot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1335  $ rwork( j ), conjg( work( j ) ) )
1336  830 continue
1337  END IF
1338  840 continue
1339 *
1340  IF( update ) THEN
1341  IF( i2.GT.0 .AND. kbt.GT.0 ) THEN
1342 *
1343 * create nonzero element a(i+kbt,i+kbt-ka-1) outside the
1344 * band and store it in WORK(m-kb+i+kbt)
1345 *
1346  work( m-kb+i+kbt ) = -bb( kbt+1, i )*ra1
1347  END IF
1348  END IF
1349 *
1350  DO 880 k = kb, 1, -1
1351  IF( update ) THEN
1352  j2 = i + k + 1 - max( 2, k+i0-m )*ka1
1353  ELSE
1354  j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1355  END IF
1356 *
1357 * finish applying rotations in 2nd set from the left
1358 *
1359  DO 850 l = kb - k, 1, -1
1360  nrt = ( j2+ka+l-1 ) / ka1
1361  j1t = j2 - ( nrt-1 )*ka1
1362  IF( nrt.GT.0 )
1363  $ CALL clartv( nrt, ab( ka1-l+1, j1t+l-1 ), inca,
1364  $ ab( ka1-l, j1t+l-1 ), inca,
1365  $ rwork( m-kb+j1t+ka ),
1366  $ work( m-kb+j1t+ka ), ka1 )
1367  850 continue
1368  nr = ( j2+ka-1 ) / ka1
1369  j1 = j2 - ( nr-1 )*ka1
1370  DO 860 j = j1, j2, ka1
1371  work( m-kb+j ) = work( m-kb+j+ka )
1372  rwork( m-kb+j ) = rwork( m-kb+j+ka )
1373  860 continue
1374  DO 870 j = j1, j2, ka1
1375 *
1376 * create nonzero element a(j+ka,j-1) outside the band
1377 * and store it in WORK(m-kb+j)
1378 *
1379  work( m-kb+j ) = work( m-kb+j )*ab( ka1, j-1 )
1380  ab( ka1, j-1 ) = rwork( m-kb+j )*ab( ka1, j-1 )
1381  870 continue
1382  IF( update ) THEN
1383  IF( i+k.GT.ka1 .AND. k.LE.kbt )
1384  $ work( m-kb+i+k-ka ) = work( m-kb+i+k )
1385  END IF
1386  880 continue
1387 *
1388  DO 920 k = kb, 1, -1
1389  j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1390  nr = ( j2+ka-1 ) / ka1
1391  j1 = j2 - ( nr-1 )*ka1
1392  IF( nr.GT.0 ) THEN
1393 *
1394 * generate rotations in 2nd set to annihilate elements
1395 * which have been created outside the band
1396 *
1397  CALL clargv( nr, ab( ka1, j1 ), inca, work( m-kb+j1 ),
1398  $ ka1, rwork( m-kb+j1 ), ka1 )
1399 *
1400 * apply rotations in 2nd set from the right
1401 *
1402  DO 890 l = 1, ka - 1
1403  CALL clartv( nr, ab( l+1, j1 ), inca, ab( l+2, j1-1 ),
1404  $ inca, rwork( m-kb+j1 ), work( m-kb+j1 ),
1405  $ ka1 )
1406  890 continue
1407 *
1408 * apply rotations in 2nd set from both sides to diagonal
1409 * blocks
1410 *
1411  CALL clar2v( nr, ab( 1, j1 ), ab( 1, j1-1 ),
1412  $ ab( 2, j1-1 ), inca, rwork( m-kb+j1 ),
1413  $ work( m-kb+j1 ), ka1 )
1414 *
1415  CALL clacgv( nr, work( m-kb+j1 ), ka1 )
1416  END IF
1417 *
1418 * start applying rotations in 2nd set from the left
1419 *
1420  DO 900 l = ka - 1, kb - k + 1, -1
1421  nrt = ( j2+l-1 ) / ka1
1422  j1t = j2 - ( nrt-1 )*ka1
1423  IF( nrt.GT.0 )
1424  $ CALL clartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,
1425  $ ab( ka1-l, j1t-ka1+l ), inca,
1426  $ rwork( m-kb+j1t ), work( m-kb+j1t ),
1427  $ ka1 )
1428  900 continue
1429 *
1430  IF( wantx ) THEN
1431 *
1432 * post-multiply X by product of rotations in 2nd set
1433 *
1434  DO 910 j = j1, j2, ka1
1435  CALL crot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1436  $ rwork( m-kb+j ), conjg( work( m-kb+j ) ) )
1437  910 continue
1438  END IF
1439  920 continue
1440 *
1441  DO 940 k = 1, kb - 1
1442  j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1443 *
1444 * finish applying rotations in 1st set from the left
1445 *
1446  DO 930 l = kb - k, 1, -1
1447  nrt = ( j2+l-1 ) / ka1
1448  j1t = j2 - ( nrt-1 )*ka1
1449  IF( nrt.GT.0 )
1450  $ CALL clartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,
1451  $ ab( ka1-l, j1t-ka1+l ), inca,
1452  $ rwork( j1t ), work( j1t ), ka1 )
1453  930 continue
1454  940 continue
1455 *
1456  IF( kb.GT.1 ) THEN
1457  DO 950 j = 2, i2 - ka
1458  rwork( j ) = rwork( j+ka )
1459  work( j ) = work( j+ka )
1460  950 continue
1461  END IF
1462 *
1463  END IF
1464 *
1465  go to 490
1466 *
1467 * End of CHBGST
1468 *
1469  END