LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ chbgst()

subroutine chbgst ( character  VECT,
character  UPLO,
integer  N,
integer  KA,
integer  KB,
complex, dimension( ldab, * )  AB,
integer  LDAB,
complex, dimension( ldbb, * )  BB,
integer  LDBB,
complex, dimension( ldx, * )  X,
integer  LDX,
complex, dimension( * )  WORK,
real, dimension( * )  RWORK,
integer  INFO 
)

CHBGST

Download CHBGST + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 CHBGST reduces a complex Hermitian-definite banded generalized
 eigenproblem  A*x = lambda*B*x  to standard form  C*y = lambda*y,
 such that C has the same bandwidth as A.

 B must have been previously factorized as S**H*S by CPBSTF, using a
 split Cholesky factorization. A is overwritten by C = X**H*A*X, where
 X = S**(-1)*Q and Q is a unitary matrix chosen to preserve the
 bandwidth of A.
Parameters
[in]VECT
          VECT is CHARACTER*1
          = 'N':  do not form the transformation matrix X;
          = 'V':  form X.
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
          The order of the matrices A and B.  N >= 0.
[in]KA
          KA is INTEGER
          The number of superdiagonals of the matrix A if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KA >= 0.
[in]KB
          KB is INTEGER
          The number of superdiagonals of the matrix B if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KA >= KB >= 0.
[in,out]AB
          AB is COMPLEX array, dimension (LDAB,N)
          On entry, the upper or lower triangle of the Hermitian band
          matrix A, stored in the first ka+1 rows of the array.  The
          j-th column of A is stored in the j-th column of the array AB
          as follows:
          if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).

          On exit, the transformed matrix X**H*A*X, stored in the same
          format as A.
[in]LDAB
          LDAB is INTEGER
          The leading dimension of the array AB.  LDAB >= KA+1.
[in]BB
          BB is COMPLEX array, dimension (LDBB,N)
          The banded factor S from the split Cholesky factorization of
          B, as returned by CPBSTF, stored in the first kb+1 rows of
          the array.
[in]LDBB
          LDBB is INTEGER
          The leading dimension of the array BB.  LDBB >= KB+1.
[out]X
          X is COMPLEX array, dimension (LDX,N)
          If VECT = 'V', the n-by-n matrix X.
          If VECT = 'N', the array X is not referenced.
[in]LDX
          LDX is INTEGER
          The leading dimension of the array X.
          LDX >= max(1,N) if VECT = 'V'; LDX >= 1 otherwise.
[out]WORK
          WORK is COMPLEX array, dimension (N)
[out]RWORK
          RWORK is REAL array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 163 of file chbgst.f.

165 *
166 * -- LAPACK computational routine --
167 * -- LAPACK is a software package provided by Univ. of Tennessee, --
168 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
169 *
170 * .. Scalar Arguments ..
171  CHARACTER UPLO, VECT
172  INTEGER INFO, KA, KB, LDAB, LDBB, LDX, N
173 * ..
174 * .. Array Arguments ..
175  REAL RWORK( * )
176  COMPLEX AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
177  $ X( LDX, * )
178 * ..
179 *
180 * =====================================================================
181 *
182 * .. Parameters ..
183  COMPLEX CZERO, CONE
184  REAL ONE
185  parameter( czero = ( 0.0e+0, 0.0e+0 ),
186  $ cone = ( 1.0e+0, 0.0e+0 ), one = 1.0e+0 )
187 * ..
188 * .. Local Scalars ..
189  LOGICAL UPDATE, UPPER, WANTX
190  INTEGER I, I0, I1, I2, INCA, J, J1, J1T, J2, J2T, K,
191  $ KA1, KB1, KBT, L, M, NR, NRT, NX
192  REAL BII
193  COMPLEX RA, RA1, T
194 * ..
195 * .. External Functions ..
196  LOGICAL LSAME
197  EXTERNAL lsame
198 * ..
199 * .. External Subroutines ..
200  EXTERNAL cgerc, cgeru, clacgv, clar2v, clargv, clartg,
202 * ..
203 * .. Intrinsic Functions ..
204  INTRINSIC conjg, max, min, real
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( 'CHBGST', -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 claset( '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 CPBSTF. 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 = real( bb( kb1, i ) )
346  ab( ka1, i ) = ( real( 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  $ conjg( ab( k-i+ka1, i ) ) -
358  $ conjg( bb( k-i+kb1, i ) )*
359  $ ab( j-i+ka1, i ) +
360  $ real( ab( ka1, i ) )*
361  $ bb( j-i+kb1, i )*
362  $ conjg( 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  $ conjg( 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 csscal( n-m, one / bii, x( m+1, i ), 1 )
382  IF( kbt.GT.0 )
383  $ CALL cgerc( 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 clartg( 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  $ conjg( 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 clargv( nrt, ab( 1, j2t ), inca, work( j2t-m ), ka1,
445  $ rwork( j2t-m ), ka1 )
446  IF( nr.GT.0 ) THEN
447 *
448 * apply rotations in 1st set from the right
449 *
450  DO 100 l = 1, ka - 1
451  CALL clartv( nr, ab( ka1-l, j2 ), inca,
452  $ ab( ka-l, j2+1 ), inca, rwork( j2-m ),
453  $ work( j2-m ), ka1 )
454  100 CONTINUE
455 *
456 * apply rotations in 1st set from both sides to diagonal
457 * blocks
458 *
459  CALL clar2v( nr, ab( ka1, j2 ), ab( ka1, j2+1 ),
460  $ ab( ka, j2+1 ), inca, rwork( j2-m ),
461  $ work( j2-m ), ka1 )
462 *
463  CALL clacgv( nr, work( j2-m ), ka1 )
464  END IF
465 *
466 * start applying rotations in 1st set from the left
467 *
468  DO 110 l = ka - 1, kb - k + 1, -1
469  nrt = ( n-j2+l ) / ka1
470  IF( nrt.GT.0 )
471  $ CALL clartv( nrt, ab( l, j2+ka1-l ), inca,
472  $ ab( l+1, j2+ka1-l ), inca, rwork( j2-m ),
473  $ work( j2-m ), ka1 )
474  110 CONTINUE
475 *
476  IF( wantx ) THEN
477 *
478 * post-multiply X by product of rotations in 1st set
479 *
480  DO 120 j = j2, j1, ka1
481  CALL crot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
482  $ rwork( j-m ), conjg( work( j-m ) ) )
483  120 CONTINUE
484  END IF
485  130 CONTINUE
486 *
487  IF( update ) THEN
488  IF( i2.LE.n .AND. kbt.GT.0 ) THEN
489 *
490 * create nonzero element a(i-kbt,i-kbt+ka+1) outside the
491 * band and store it in WORK(i-kbt)
492 *
493  work( i-kbt ) = -bb( kb1-kbt, i )*ra1
494  END IF
495  END IF
496 *
497  DO 170 k = kb, 1, -1
498  IF( update ) THEN
499  j2 = i - k - 1 + max( 2, k-i0+1 )*ka1
500  ELSE
501  j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
502  END IF
503 *
504 * finish applying rotations in 2nd set from the left
505 *
506  DO 140 l = kb - k, 1, -1
507  nrt = ( n-j2+ka+l ) / ka1
508  IF( nrt.GT.0 )
509  $ CALL clartv( nrt, ab( l, j2-l+1 ), inca,
510  $ ab( l+1, j2-l+1 ), inca, rwork( j2-ka ),
511  $ work( j2-ka ), ka1 )
512  140 CONTINUE
513  nr = ( n-j2+ka ) / ka1
514  j1 = j2 + ( nr-1 )*ka1
515  DO 150 j = j1, j2, -ka1
516  work( j ) = work( j-ka )
517  rwork( j ) = rwork( j-ka )
518  150 CONTINUE
519  DO 160 j = j2, j1, ka1
520 *
521 * create nonzero element a(j-ka,j+1) outside the band
522 * and store it in WORK(j)
523 *
524  work( j ) = work( j )*ab( 1, j+1 )
525  ab( 1, j+1 ) = rwork( j )*ab( 1, j+1 )
526  160 CONTINUE
527  IF( update ) THEN
528  IF( i-k.LT.n-ka .AND. k.LE.kbt )
529  $ work( i-k+ka ) = work( i-k )
530  END IF
531  170 CONTINUE
532 *
533  DO 210 k = kb, 1, -1
534  j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
535  nr = ( n-j2+ka ) / ka1
536  j1 = j2 + ( nr-1 )*ka1
537  IF( nr.GT.0 ) THEN
538 *
539 * generate rotations in 2nd set to annihilate elements
540 * which have been created outside the band
541 *
542  CALL clargv( nr, ab( 1, j2 ), inca, work( j2 ), ka1,
543  $ rwork( j2 ), ka1 )
544 *
545 * apply rotations in 2nd set from the right
546 *
547  DO 180 l = 1, ka - 1
548  CALL clartv( nr, ab( ka1-l, j2 ), inca,
549  $ ab( ka-l, j2+1 ), inca, rwork( j2 ),
550  $ work( j2 ), ka1 )
551  180 CONTINUE
552 *
553 * apply rotations in 2nd set from both sides to diagonal
554 * blocks
555 *
556  CALL clar2v( nr, ab( ka1, j2 ), ab( ka1, j2+1 ),
557  $ ab( ka, j2+1 ), inca, rwork( j2 ),
558  $ work( j2 ), ka1 )
559 *
560  CALL clacgv( nr, work( j2 ), ka1 )
561  END IF
562 *
563 * start applying rotations in 2nd set from the left
564 *
565  DO 190 l = ka - 1, kb - k + 1, -1
566  nrt = ( n-j2+l ) / ka1
567  IF( nrt.GT.0 )
568  $ CALL clartv( nrt, ab( l, j2+ka1-l ), inca,
569  $ ab( l+1, j2+ka1-l ), inca, rwork( j2 ),
570  $ work( j2 ), ka1 )
571  190 CONTINUE
572 *
573  IF( wantx ) THEN
574 *
575 * post-multiply X by product of rotations in 2nd set
576 *
577  DO 200 j = j2, j1, ka1
578  CALL crot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
579  $ rwork( j ), conjg( work( j ) ) )
580  200 CONTINUE
581  END IF
582  210 CONTINUE
583 *
584  DO 230 k = 1, kb - 1
585  j2 = i - k - 1 + max( 1, k-i0+2 )*ka1
586 *
587 * finish applying rotations in 1st set from the left
588 *
589  DO 220 l = kb - k, 1, -1
590  nrt = ( n-j2+l ) / ka1
591  IF( nrt.GT.0 )
592  $ CALL clartv( nrt, ab( l, j2+ka1-l ), inca,
593  $ ab( l+1, j2+ka1-l ), inca, rwork( j2-m ),
594  $ work( j2-m ), ka1 )
595  220 CONTINUE
596  230 CONTINUE
597 *
598  IF( kb.GT.1 ) THEN
599  DO 240 j = n - 1, j2 + ka, -1
600  rwork( j-m ) = rwork( j-ka-m )
601  work( j-m ) = work( j-ka-m )
602  240 CONTINUE
603  END IF
604 *
605  ELSE
606 *
607 * Transform A, working with the lower triangle
608 *
609  IF( update ) THEN
610 *
611 * Form inv(S(i))**H * A * inv(S(i))
612 *
613  bii = real( bb( 1, i ) )
614  ab( 1, i ) = ( real( ab( 1, i ) ) / bii ) / bii
615  DO 250 j = i + 1, i1
616  ab( j-i+1, i ) = ab( j-i+1, i ) / bii
617  250 CONTINUE
618  DO 260 j = max( 1, i-ka ), i - 1
619  ab( i-j+1, j ) = ab( i-j+1, j ) / bii
620  260 CONTINUE
621  DO 290 k = i - kbt, i - 1
622  DO 270 j = i - kbt, k
623  ab( k-j+1, j ) = ab( k-j+1, j ) -
624  $ bb( i-j+1, j )*conjg( ab( i-k+1,
625  $ k ) ) - conjg( bb( i-k+1, k ) )*
626  $ ab( i-j+1, j ) + real( ab( 1, i ) )*
627  $ bb( i-j+1, j )*conjg( bb( i-k+1,
628  $ k ) )
629  270 CONTINUE
630  DO 280 j = max( 1, i-ka ), i - kbt - 1
631  ab( k-j+1, j ) = ab( k-j+1, j ) -
632  $ conjg( bb( i-k+1, k ) )*
633  $ ab( i-j+1, j )
634  280 CONTINUE
635  290 CONTINUE
636  DO 310 j = i, i1
637  DO 300 k = max( j-ka, i-kbt ), i - 1
638  ab( j-k+1, k ) = ab( j-k+1, k ) -
639  $ bb( i-k+1, k )*ab( j-i+1, i )
640  300 CONTINUE
641  310 CONTINUE
642 *
643  IF( wantx ) THEN
644 *
645 * post-multiply X by inv(S(i))
646 *
647  CALL csscal( n-m, one / bii, x( m+1, i ), 1 )
648  IF( kbt.GT.0 )
649  $ CALL cgeru( n-m, kbt, -cone, x( m+1, i ), 1,
650  $ bb( kbt+1, i-kbt ), ldbb-1,
651  $ x( m+1, i-kbt ), ldx )
652  END IF
653 *
654 * store a(i1,i) in RA1 for use in next loop over K
655 *
656  ra1 = ab( i1-i+1, i )
657  END IF
658 *
659 * Generate and apply vectors of rotations to chase all the
660 * existing bulges KA positions down toward the bottom of the
661 * band
662 *
663  DO 360 k = 1, kb - 1
664  IF( update ) THEN
665 *
666 * Determine the rotations which would annihilate the bulge
667 * which has in theory just been created
668 *
669  IF( i-k+ka.LT.n .AND. i-k.GT.1 ) THEN
670 *
671 * generate rotation to annihilate a(i-k+ka+1,i)
672 *
673  CALL clartg( ab( ka1-k, i ), ra1, rwork( i-k+ka-m ),
674  $ work( i-k+ka-m ), ra )
675 *
676 * create nonzero element a(i-k+ka+1,i-k) outside the
677 * band and store it in WORK(i-k)
678 *
679  t = -bb( k+1, i-k )*ra1
680  work( i-k ) = rwork( i-k+ka-m )*t -
681  $ conjg( work( i-k+ka-m ) )*ab( ka1, i-k )
682  ab( ka1, i-k ) = work( i-k+ka-m )*t +
683  $ rwork( i-k+ka-m )*ab( ka1, i-k )
684  ra1 = ra
685  END IF
686  END IF
687  j2 = i - k - 1 + max( 1, k-i0+2 )*ka1
688  nr = ( n-j2+ka ) / ka1
689  j1 = j2 + ( nr-1 )*ka1
690  IF( update ) THEN
691  j2t = max( j2, i+2*ka-k+1 )
692  ELSE
693  j2t = j2
694  END IF
695  nrt = ( n-j2t+ka ) / ka1
696  DO 320 j = j2t, j1, ka1
697 *
698 * create nonzero element a(j+1,j-ka) outside the band
699 * and store it in WORK(j-m)
700 *
701  work( j-m ) = work( j-m )*ab( ka1, j-ka+1 )
702  ab( ka1, j-ka+1 ) = rwork( j-m )*ab( ka1, j-ka+1 )
703  320 CONTINUE
704 *
705 * generate rotations in 1st set to annihilate elements which
706 * have been created outside the band
707 *
708  IF( nrt.GT.0 )
709  $ CALL clargv( nrt, ab( ka1, j2t-ka ), inca, work( j2t-m ),
710  $ ka1, rwork( j2t-m ), ka1 )
711  IF( nr.GT.0 ) THEN
712 *
713 * apply rotations in 1st set from the left
714 *
715  DO 330 l = 1, ka - 1
716  CALL clartv( nr, ab( l+1, j2-l ), inca,
717  $ ab( l+2, j2-l ), inca, rwork( j2-m ),
718  $ work( j2-m ), ka1 )
719  330 CONTINUE
720 *
721 * apply rotations in 1st set from both sides to diagonal
722 * blocks
723 *
724  CALL clar2v( nr, ab( 1, j2 ), ab( 1, j2+1 ), ab( 2, j2 ),
725  $ inca, rwork( j2-m ), work( j2-m ), ka1 )
726 *
727  CALL clacgv( nr, work( j2-m ), ka1 )
728  END IF
729 *
730 * start applying rotations in 1st set from the right
731 *
732  DO 340 l = ka - 1, kb - k + 1, -1
733  nrt = ( n-j2+l ) / ka1
734  IF( nrt.GT.0 )
735  $ CALL clartv( nrt, ab( ka1-l+1, j2 ), inca,
736  $ ab( ka1-l, j2+1 ), inca, rwork( j2-m ),
737  $ work( j2-m ), ka1 )
738  340 CONTINUE
739 *
740  IF( wantx ) THEN
741 *
742 * post-multiply X by product of rotations in 1st set
743 *
744  DO 350 j = j2, j1, ka1
745  CALL crot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
746  $ rwork( j-m ), work( j-m ) )
747  350 CONTINUE
748  END IF
749  360 CONTINUE
750 *
751  IF( update ) THEN
752  IF( i2.LE.n .AND. kbt.GT.0 ) THEN
753 *
754 * create nonzero element a(i-kbt+ka+1,i-kbt) outside the
755 * band and store it in WORK(i-kbt)
756 *
757  work( i-kbt ) = -bb( kbt+1, i-kbt )*ra1
758  END IF
759  END IF
760 *
761  DO 400 k = kb, 1, -1
762  IF( update ) THEN
763  j2 = i - k - 1 + max( 2, k-i0+1 )*ka1
764  ELSE
765  j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
766  END IF
767 *
768 * finish applying rotations in 2nd set from the right
769 *
770  DO 370 l = kb - k, 1, -1
771  nrt = ( n-j2+ka+l ) / ka1
772  IF( nrt.GT.0 )
773  $ CALL clartv( nrt, ab( ka1-l+1, j2-ka ), inca,
774  $ ab( ka1-l, j2-ka+1 ), inca,
775  $ rwork( j2-ka ), work( j2-ka ), ka1 )
776  370 CONTINUE
777  nr = ( n-j2+ka ) / ka1
778  j1 = j2 + ( nr-1 )*ka1
779  DO 380 j = j1, j2, -ka1
780  work( j ) = work( j-ka )
781  rwork( j ) = rwork( j-ka )
782  380 CONTINUE
783  DO 390 j = j2, j1, ka1
784 *
785 * create nonzero element a(j+1,j-ka) outside the band
786 * and store it in WORK(j)
787 *
788  work( j ) = work( j )*ab( ka1, j-ka+1 )
789  ab( ka1, j-ka+1 ) = rwork( j )*ab( ka1, j-ka+1 )
790  390 CONTINUE
791  IF( update ) THEN
792  IF( i-k.LT.n-ka .AND. k.LE.kbt )
793  $ work( i-k+ka ) = work( i-k )
794  END IF
795  400 CONTINUE
796 *
797  DO 440 k = kb, 1, -1
798  j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
799  nr = ( n-j2+ka ) / ka1
800  j1 = j2 + ( nr-1 )*ka1
801  IF( nr.GT.0 ) THEN
802 *
803 * generate rotations in 2nd set to annihilate elements
804 * which have been created outside the band
805 *
806  CALL clargv( nr, ab( ka1, j2-ka ), inca, work( j2 ), ka1,
807  $ rwork( j2 ), ka1 )
808 *
809 * apply rotations in 2nd set from the left
810 *
811  DO 410 l = 1, ka - 1
812  CALL clartv( nr, ab( l+1, j2-l ), inca,
813  $ ab( l+2, j2-l ), inca, rwork( j2 ),
814  $ work( j2 ), ka1 )
815  410 CONTINUE
816 *
817 * apply rotations in 2nd set from both sides to diagonal
818 * blocks
819 *
820  CALL clar2v( nr, ab( 1, j2 ), ab( 1, j2+1 ), ab( 2, j2 ),
821  $ inca, rwork( j2 ), work( j2 ), ka1 )
822 *
823  CALL clacgv( nr, work( j2 ), ka1 )
824  END IF
825 *
826 * start applying rotations in 2nd set from the right
827 *
828  DO 420 l = ka - 1, kb - k + 1, -1
829  nrt = ( n-j2+l ) / ka1
830  IF( nrt.GT.0 )
831  $ CALL clartv( nrt, ab( ka1-l+1, j2 ), inca,
832  $ ab( ka1-l, j2+1 ), inca, rwork( j2 ),
833  $ work( j2 ), ka1 )
834  420 CONTINUE
835 *
836  IF( wantx ) THEN
837 *
838 * post-multiply X by product of rotations in 2nd set
839 *
840  DO 430 j = j2, j1, ka1
841  CALL crot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
842  $ rwork( j ), work( j ) )
843  430 CONTINUE
844  END IF
845  440 CONTINUE
846 *
847  DO 460 k = 1, kb - 1
848  j2 = i - k - 1 + max( 1, k-i0+2 )*ka1
849 *
850 * finish applying rotations in 1st set from the right
851 *
852  DO 450 l = kb - k, 1, -1
853  nrt = ( n-j2+l ) / ka1
854  IF( nrt.GT.0 )
855  $ CALL clartv( nrt, ab( ka1-l+1, j2 ), inca,
856  $ ab( ka1-l, j2+1 ), inca, rwork( j2-m ),
857  $ work( j2-m ), ka1 )
858  450 CONTINUE
859  460 CONTINUE
860 *
861  IF( kb.GT.1 ) THEN
862  DO 470 j = n - 1, j2 + ka, -1
863  rwork( j-m ) = rwork( j-ka-m )
864  work( j-m ) = work( j-ka-m )
865  470 CONTINUE
866  END IF
867 *
868  END IF
869 *
870  GO TO 10
871 *
872  480 CONTINUE
873 *
874 * **************************** Phase 2 *****************************
875 *
876 * The logical structure of this phase is:
877 *
878 * UPDATE = .TRUE.
879 * DO I = 1, M
880 * use S(i) to update A and create a new bulge
881 * apply rotations to push all bulges KA positions upward
882 * END DO
883 * UPDATE = .FALSE.
884 * DO I = M - KA - 1, 2, -1
885 * apply rotations to push all bulges KA positions upward
886 * END DO
887 *
888 * To avoid duplicating code, the two loops are merged.
889 *
890  update = .true.
891  i = 0
892  490 CONTINUE
893  IF( update ) THEN
894  i = i + 1
895  kbt = min( kb, m-i )
896  i0 = i + 1
897  i1 = max( 1, i-ka )
898  i2 = i + kbt - ka1
899  IF( i.GT.m ) THEN
900  update = .false.
901  i = i - 1
902  i0 = m + 1
903  IF( ka.EQ.0 )
904  $ RETURN
905  GO TO 490
906  END IF
907  ELSE
908  i = i - ka
909  IF( i.LT.2 )
910  $ RETURN
911  END IF
912 *
913  IF( i.LT.m-kbt ) THEN
914  nx = m
915  ELSE
916  nx = n
917  END IF
918 *
919  IF( upper ) THEN
920 *
921 * Transform A, working with the upper triangle
922 *
923  IF( update ) THEN
924 *
925 * Form inv(S(i))**H * A * inv(S(i))
926 *
927  bii = real( bb( kb1, i ) )
928  ab( ka1, i ) = ( real( ab( ka1, i ) ) / bii ) / bii
929  DO 500 j = i1, i - 1
930  ab( j-i+ka1, i ) = ab( j-i+ka1, i ) / bii
931  500 CONTINUE
932  DO 510 j = i + 1, min( n, i+ka )
933  ab( i-j+ka1, j ) = ab( i-j+ka1, j ) / bii
934  510 CONTINUE
935  DO 540 k = i + 1, i + kbt
936  DO 520 j = k, i + kbt
937  ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -
938  $ bb( i-j+kb1, j )*
939  $ conjg( ab( i-k+ka1, k ) ) -
940  $ conjg( bb( i-k+kb1, k ) )*
941  $ ab( i-j+ka1, j ) +
942  $ real( ab( ka1, i ) )*
943  $ bb( i-j+kb1, j )*
944  $ conjg( bb( i-k+kb1, k ) )
945  520 CONTINUE
946  DO 530 j = i + kbt + 1, min( n, i+ka )
947  ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -
948  $ conjg( bb( i-k+kb1, k ) )*
949  $ ab( i-j+ka1, j )
950  530 CONTINUE
951  540 CONTINUE
952  DO 560 j = i1, i
953  DO 550 k = i + 1, min( j+ka, i+kbt )
954  ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -
955  $ bb( i-k+kb1, k )*ab( j-i+ka1, i )
956  550 CONTINUE
957  560 CONTINUE
958 *
959  IF( wantx ) THEN
960 *
961 * post-multiply X by inv(S(i))
962 *
963  CALL csscal( nx, one / bii, x( 1, i ), 1 )
964  IF( kbt.GT.0 )
965  $ CALL cgeru( nx, kbt, -cone, x( 1, i ), 1,
966  $ bb( kb, i+1 ), ldbb-1, x( 1, i+1 ), ldx )
967  END IF
968 *
969 * store a(i1,i) in RA1 for use in next loop over K
970 *
971  ra1 = ab( i1-i+ka1, i )
972  END IF
973 *
974 * Generate and apply vectors of rotations to chase all the
975 * existing bulges KA positions up toward the top of the band
976 *
977  DO 610 k = 1, kb - 1
978  IF( update ) THEN
979 *
980 * Determine the rotations which would annihilate the bulge
981 * which has in theory just been created
982 *
983  IF( i+k-ka1.GT.0 .AND. i+k.LT.m ) THEN
984 *
985 * generate rotation to annihilate a(i+k-ka-1,i)
986 *
987  CALL clartg( ab( k+1, i ), ra1, rwork( i+k-ka ),
988  $ work( i+k-ka ), ra )
989 *
990 * create nonzero element a(i+k-ka-1,i+k) outside the
991 * band and store it in WORK(m-kb+i+k)
992 *
993  t = -bb( kb1-k, i+k )*ra1
994  work( m-kb+i+k ) = rwork( i+k-ka )*t -
995  $ conjg( work( i+k-ka ) )*
996  $ ab( 1, i+k )
997  ab( 1, i+k ) = work( i+k-ka )*t +
998  $ rwork( i+k-ka )*ab( 1, i+k )
999  ra1 = ra
1000  END IF
1001  END IF
1002  j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1003  nr = ( j2+ka-1 ) / ka1
1004  j1 = j2 - ( nr-1 )*ka1
1005  IF( update ) THEN
1006  j2t = min( j2, i-2*ka+k-1 )
1007  ELSE
1008  j2t = j2
1009  END IF
1010  nrt = ( j2t+ka-1 ) / ka1
1011  DO 570 j = j1, j2t, ka1
1012 *
1013 * create nonzero element a(j-1,j+ka) outside the band
1014 * and store it in WORK(j)
1015 *
1016  work( j ) = work( j )*ab( 1, j+ka-1 )
1017  ab( 1, j+ka-1 ) = rwork( j )*ab( 1, j+ka-1 )
1018  570 CONTINUE
1019 *
1020 * generate rotations in 1st set to annihilate elements which
1021 * have been created outside the band
1022 *
1023  IF( nrt.GT.0 )
1024  $ CALL clargv( nrt, ab( 1, j1+ka ), inca, work( j1 ), ka1,
1025  $ rwork( j1 ), ka1 )
1026  IF( nr.GT.0 ) THEN
1027 *
1028 * apply rotations in 1st set from the left
1029 *
1030  DO 580 l = 1, ka - 1
1031  CALL clartv( nr, ab( ka1-l, j1+l ), inca,
1032  $ ab( ka-l, j1+l ), inca, rwork( j1 ),
1033  $ work( j1 ), ka1 )
1034  580 CONTINUE
1035 *
1036 * apply rotations in 1st set from both sides to diagonal
1037 * blocks
1038 *
1039  CALL clar2v( nr, ab( ka1, j1 ), ab( ka1, j1-1 ),
1040  $ ab( ka, j1 ), inca, rwork( j1 ), work( j1 ),
1041  $ ka1 )
1042 *
1043  CALL clacgv( nr, work( j1 ), ka1 )
1044  END IF
1045 *
1046 * start applying rotations in 1st set from the right
1047 *
1048  DO 590 l = ka - 1, kb - k + 1, -1
1049  nrt = ( j2+l-1 ) / ka1
1050  j1t = j2 - ( nrt-1 )*ka1
1051  IF( nrt.GT.0 )
1052  $ CALL clartv( nrt, ab( l, j1t ), inca,
1053  $ ab( l+1, j1t-1 ), inca, rwork( j1t ),
1054  $ work( j1t ), ka1 )
1055  590 CONTINUE
1056 *
1057  IF( wantx ) THEN
1058 *
1059 * post-multiply X by product of rotations in 1st set
1060 *
1061  DO 600 j = j1, j2, ka1
1062  CALL crot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1063  $ rwork( j ), work( j ) )
1064  600 CONTINUE
1065  END IF
1066  610 CONTINUE
1067 *
1068  IF( update ) THEN
1069  IF( i2.GT.0 .AND. kbt.GT.0 ) THEN
1070 *
1071 * create nonzero element a(i+kbt-ka-1,i+kbt) outside the
1072 * band and store it in WORK(m-kb+i+kbt)
1073 *
1074  work( m-kb+i+kbt ) = -bb( kb1-kbt, i+kbt )*ra1
1075  END IF
1076  END IF
1077 *
1078  DO 650 k = kb, 1, -1
1079  IF( update ) THEN
1080  j2 = i + k + 1 - max( 2, k+i0-m )*ka1
1081  ELSE
1082  j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1083  END IF
1084 *
1085 * finish applying rotations in 2nd set from the right
1086 *
1087  DO 620 l = kb - k, 1, -1
1088  nrt = ( j2+ka+l-1 ) / ka1
1089  j1t = j2 - ( nrt-1 )*ka1
1090  IF( nrt.GT.0 )
1091  $ CALL clartv( nrt, ab( l, j1t+ka ), inca,
1092  $ ab( l+1, j1t+ka-1 ), inca,
1093  $ rwork( m-kb+j1t+ka ),
1094  $ work( m-kb+j1t+ka ), ka1 )
1095  620 CONTINUE
1096  nr = ( j2+ka-1 ) / ka1
1097  j1 = j2 - ( nr-1 )*ka1
1098  DO 630 j = j1, j2, ka1
1099  work( m-kb+j ) = work( m-kb+j+ka )
1100  rwork( m-kb+j ) = rwork( m-kb+j+ka )
1101  630 CONTINUE
1102  DO 640 j = j1, j2, ka1
1103 *
1104 * create nonzero element a(j-1,j+ka) outside the band
1105 * and store it in WORK(m-kb+j)
1106 *
1107  work( m-kb+j ) = work( m-kb+j )*ab( 1, j+ka-1 )
1108  ab( 1, j+ka-1 ) = rwork( m-kb+j )*ab( 1, j+ka-1 )
1109  640 CONTINUE
1110  IF( update ) THEN
1111  IF( i+k.GT.ka1 .AND. k.LE.kbt )
1112  $ work( m-kb+i+k-ka ) = work( m-kb+i+k )
1113  END IF
1114  650 CONTINUE
1115 *
1116  DO 690 k = kb, 1, -1
1117  j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1118  nr = ( j2+ka-1 ) / ka1
1119  j1 = j2 - ( nr-1 )*ka1
1120  IF( nr.GT.0 ) THEN
1121 *
1122 * generate rotations in 2nd set to annihilate elements
1123 * which have been created outside the band
1124 *
1125  CALL clargv( nr, ab( 1, j1+ka ), inca, work( m-kb+j1 ),
1126  $ ka1, rwork( m-kb+j1 ), ka1 )
1127 *
1128 * apply rotations in 2nd set from the left
1129 *
1130  DO 660 l = 1, ka - 1
1131  CALL clartv( nr, ab( ka1-l, j1+l ), inca,
1132  $ ab( ka-l, j1+l ), inca, rwork( m-kb+j1 ),
1133  $ work( m-kb+j1 ), ka1 )
1134  660 CONTINUE
1135 *
1136 * apply rotations in 2nd set from both sides to diagonal
1137 * blocks
1138 *
1139  CALL clar2v( nr, ab( ka1, j1 ), ab( ka1, j1-1 ),
1140  $ ab( ka, j1 ), inca, rwork( m-kb+j1 ),
1141  $ work( m-kb+j1 ), ka1 )
1142 *
1143  CALL clacgv( nr, work( m-kb+j1 ), ka1 )
1144  END IF
1145 *
1146 * start applying rotations in 2nd set from the right
1147 *
1148  DO 670 l = ka - 1, kb - k + 1, -1
1149  nrt = ( j2+l-1 ) / ka1
1150  j1t = j2 - ( nrt-1 )*ka1
1151  IF( nrt.GT.0 )
1152  $ CALL clartv( nrt, ab( l, j1t ), inca,
1153  $ ab( l+1, j1t-1 ), inca,
1154  $ rwork( m-kb+j1t ), work( m-kb+j1t ),
1155  $ ka1 )
1156  670 CONTINUE
1157 *
1158  IF( wantx ) THEN
1159 *
1160 * post-multiply X by product of rotations in 2nd set
1161 *
1162  DO 680 j = j1, j2, ka1
1163  CALL crot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1164  $ rwork( m-kb+j ), work( m-kb+j ) )
1165  680 CONTINUE
1166  END IF
1167  690 CONTINUE
1168 *
1169  DO 710 k = 1, kb - 1
1170  j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1171 *
1172 * finish applying rotations in 1st set from the right
1173 *
1174  DO 700 l = kb - k, 1, -1
1175  nrt = ( j2+l-1 ) / ka1
1176  j1t = j2 - ( nrt-1 )*ka1
1177  IF( nrt.GT.0 )
1178  $ CALL clartv( nrt, ab( l, j1t ), inca,
1179  $ ab( l+1, j1t-1 ), inca, rwork( j1t ),
1180  $ work( j1t ), ka1 )
1181  700 CONTINUE
1182  710 CONTINUE
1183 *
1184  IF( kb.GT.1 ) THEN
1185  DO 720 j = 2, i2 - ka
1186  rwork( j ) = rwork( j+ka )
1187  work( j ) = work( j+ka )
1188  720 CONTINUE
1189  END IF
1190 *
1191  ELSE
1192 *
1193 * Transform A, working with the lower triangle
1194 *
1195  IF( update ) THEN
1196 *
1197 * Form inv(S(i))**H * A * inv(S(i))
1198 *
1199  bii = real( bb( 1, i ) )
1200  ab( 1, i ) = ( real( ab( 1, i ) ) / bii ) / bii
1201  DO 730 j = i1, i - 1
1202  ab( i-j+1, j ) = ab( i-j+1, j ) / bii
1203  730 CONTINUE
1204  DO 740 j = i + 1, min( n, i+ka )
1205  ab( j-i+1, i ) = ab( j-i+1, i ) / bii
1206  740 CONTINUE
1207  DO 770 k = i + 1, i + kbt
1208  DO 750 j = k, i + kbt
1209  ab( j-k+1, k ) = ab( j-k+1, k ) -
1210  $ bb( j-i+1, i )*conjg( ab( k-i+1,
1211  $ i ) ) - conjg( bb( k-i+1, i ) )*
1212  $ ab( j-i+1, i ) + real( ab( 1, i ) )*
1213  $ bb( j-i+1, i )*conjg( bb( k-i+1,
1214  $ i ) )
1215  750 CONTINUE
1216  DO 760 j = i + kbt + 1, min( n, i+ka )
1217  ab( j-k+1, k ) = ab( j-k+1, k ) -
1218  $ conjg( bb( k-i+1, i ) )*
1219  $ ab( j-i+1, i )
1220  760 CONTINUE
1221  770 CONTINUE
1222  DO 790 j = i1, i
1223  DO 780 k = i + 1, min( j+ka, i+kbt )
1224  ab( k-j+1, j ) = ab( k-j+1, j ) -
1225  $ bb( k-i+1, i )*ab( i-j+1, j )
1226  780 CONTINUE
1227  790 CONTINUE
1228 *
1229  IF( wantx ) THEN
1230 *
1231 * post-multiply X by inv(S(i))
1232 *
1233  CALL csscal( nx, one / bii, x( 1, i ), 1 )
1234  IF( kbt.GT.0 )
1235  $ CALL cgerc( nx, kbt, -cone, x( 1, i ), 1, bb( 2, i ),
1236  $ 1, x( 1, i+1 ), ldx )
1237  END IF
1238 *
1239 * store a(i,i1) in RA1 for use in next loop over K
1240 *
1241  ra1 = ab( i-i1+1, i1 )
1242  END IF
1243 *
1244 * Generate and apply vectors of rotations to chase all the
1245 * existing bulges KA positions up toward the top of the band
1246 *
1247  DO 840 k = 1, kb - 1
1248  IF( update ) THEN
1249 *
1250 * Determine the rotations which would annihilate the bulge
1251 * which has in theory just been created
1252 *
1253  IF( i+k-ka1.GT.0 .AND. i+k.LT.m ) THEN
1254 *
1255 * generate rotation to annihilate a(i,i+k-ka-1)
1256 *
1257  CALL clartg( ab( ka1-k, i+k-ka ), ra1,
1258  $ rwork( i+k-ka ), work( i+k-ka ), ra )
1259 *
1260 * create nonzero element a(i+k,i+k-ka-1) outside the
1261 * band and store it in WORK(m-kb+i+k)
1262 *
1263  t = -bb( k+1, i )*ra1
1264  work( m-kb+i+k ) = rwork( i+k-ka )*t -
1265  $ conjg( work( i+k-ka ) )*
1266  $ ab( ka1, i+k-ka )
1267  ab( ka1, i+k-ka ) = work( i+k-ka )*t +
1268  $ rwork( i+k-ka )*ab( ka1, i+k-ka )
1269  ra1 = ra
1270  END IF
1271  END IF
1272  j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1273  nr = ( j2+ka-1 ) / ka1
1274  j1 = j2 - ( nr-1 )*ka1
1275  IF( update ) THEN
1276  j2t = min( j2, i-2*ka+k-1 )
1277  ELSE
1278  j2t = j2
1279  END IF
1280  nrt = ( j2t+ka-1 ) / ka1
1281  DO 800 j = j1, j2t, ka1
1282 *
1283 * create nonzero element a(j+ka,j-1) outside the band
1284 * and store it in WORK(j)
1285 *
1286  work( j ) = work( j )*ab( ka1, j-1 )
1287  ab( ka1, j-1 ) = rwork( j )*ab( ka1, j-1 )
1288  800 CONTINUE
1289 *
1290 * generate rotations in 1st set to annihilate elements which
1291 * have been created outside the band
1292 *
1293  IF( nrt.GT.0 )
1294  $ CALL clargv( nrt, ab( ka1, j1 ), inca, work( j1 ), ka1,
1295  $ rwork( j1 ), ka1 )
1296  IF( nr.GT.0 ) THEN
1297 *
1298 * apply rotations in 1st set from the right
1299 *
1300  DO 810 l = 1, ka - 1
1301  CALL clartv( nr, ab( l+1, j1 ), inca, ab( l+2, j1-1 ),
1302  $ inca, rwork( j1 ), work( j1 ), ka1 )
1303  810 CONTINUE
1304 *
1305 * apply rotations in 1st set from both sides to diagonal
1306 * blocks
1307 *
1308  CALL clar2v( nr, ab( 1, j1 ), ab( 1, j1-1 ),
1309  $ ab( 2, j1-1 ), inca, rwork( j1 ),
1310  $ work( j1 ), ka1 )
1311 *
1312  CALL clacgv( nr, work( j1 ), ka1 )
1313  END IF
1314 *
1315 * start applying rotations in 1st set from the left
1316 *
1317  DO 820 l = ka - 1, kb - k + 1, -1
1318  nrt = ( j2+l-1 ) / ka1
1319  j1t = j2 - ( nrt-1 )*ka1
1320  IF( nrt.GT.0 )
1321  $ CALL clartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,
1322  $ ab( ka1-l, j1t-ka1+l ), inca,
1323  $ rwork( j1t ), work( j1t ), ka1 )
1324  820 CONTINUE
1325 *
1326  IF( wantx ) THEN
1327 *
1328 * post-multiply X by product of rotations in 1st set
1329 *
1330  DO 830 j = j1, j2, ka1
1331  CALL crot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1332  $ rwork( j ), conjg( work( j ) ) )
1333  830 CONTINUE
1334  END IF
1335  840 CONTINUE
1336 *
1337  IF( update ) THEN
1338  IF( i2.GT.0 .AND. kbt.GT.0 ) THEN
1339 *
1340 * create nonzero element a(i+kbt,i+kbt-ka-1) outside the
1341 * band and store it in WORK(m-kb+i+kbt)
1342 *
1343  work( m-kb+i+kbt ) = -bb( kbt+1, i )*ra1
1344  END IF
1345  END IF
1346 *
1347  DO 880 k = kb, 1, -1
1348  IF( update ) THEN
1349  j2 = i + k + 1 - max( 2, k+i0-m )*ka1
1350  ELSE
1351  j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1352  END IF
1353 *
1354 * finish applying rotations in 2nd set from the left
1355 *
1356  DO 850 l = kb - k, 1, -1
1357  nrt = ( j2+ka+l-1 ) / ka1
1358  j1t = j2 - ( nrt-1 )*ka1
1359  IF( nrt.GT.0 )
1360  $ CALL clartv( nrt, ab( ka1-l+1, j1t+l-1 ), inca,
1361  $ ab( ka1-l, j1t+l-1 ), inca,
1362  $ rwork( m-kb+j1t+ka ),
1363  $ work( m-kb+j1t+ka ), ka1 )
1364  850 CONTINUE
1365  nr = ( j2+ka-1 ) / ka1
1366  j1 = j2 - ( nr-1 )*ka1
1367  DO 860 j = j1, j2, ka1
1368  work( m-kb+j ) = work( m-kb+j+ka )
1369  rwork( m-kb+j ) = rwork( m-kb+j+ka )
1370  860 CONTINUE
1371  DO 870 j = j1, j2, ka1
1372 *
1373 * create nonzero element a(j+ka,j-1) outside the band
1374 * and store it in WORK(m-kb+j)
1375 *
1376  work( m-kb+j ) = work( m-kb+j )*ab( ka1, j-1 )
1377  ab( ka1, j-1 ) = rwork( m-kb+j )*ab( ka1, j-1 )
1378  870 CONTINUE
1379  IF( update ) THEN
1380  IF( i+k.GT.ka1 .AND. k.LE.kbt )
1381  $ work( m-kb+i+k-ka ) = work( m-kb+i+k )
1382  END IF
1383  880 CONTINUE
1384 *
1385  DO 920 k = kb, 1, -1
1386  j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1387  nr = ( j2+ka-1 ) / ka1
1388  j1 = j2 - ( nr-1 )*ka1
1389  IF( nr.GT.0 ) THEN
1390 *
1391 * generate rotations in 2nd set to annihilate elements
1392 * which have been created outside the band
1393 *
1394  CALL clargv( nr, ab( ka1, j1 ), inca, work( m-kb+j1 ),
1395  $ ka1, rwork( m-kb+j1 ), ka1 )
1396 *
1397 * apply rotations in 2nd set from the right
1398 *
1399  DO 890 l = 1, ka - 1
1400  CALL clartv( nr, ab( l+1, j1 ), inca, ab( l+2, j1-1 ),
1401  $ inca, rwork( m-kb+j1 ), work( m-kb+j1 ),
1402  $ ka1 )
1403  890 CONTINUE
1404 *
1405 * apply rotations in 2nd set from both sides to diagonal
1406 * blocks
1407 *
1408  CALL clar2v( nr, ab( 1, j1 ), ab( 1, j1-1 ),
1409  $ ab( 2, j1-1 ), inca, rwork( m-kb+j1 ),
1410  $ work( m-kb+j1 ), ka1 )
1411 *
1412  CALL clacgv( nr, work( m-kb+j1 ), ka1 )
1413  END IF
1414 *
1415 * start applying rotations in 2nd set from the left
1416 *
1417  DO 900 l = ka - 1, kb - k + 1, -1
1418  nrt = ( j2+l-1 ) / ka1
1419  j1t = j2 - ( nrt-1 )*ka1
1420  IF( nrt.GT.0 )
1421  $ CALL clartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,
1422  $ ab( ka1-l, j1t-ka1+l ), inca,
1423  $ rwork( m-kb+j1t ), work( m-kb+j1t ),
1424  $ ka1 )
1425  900 CONTINUE
1426 *
1427  IF( wantx ) THEN
1428 *
1429 * post-multiply X by product of rotations in 2nd set
1430 *
1431  DO 910 j = j1, j2, ka1
1432  CALL crot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1433  $ rwork( m-kb+j ), conjg( work( m-kb+j ) ) )
1434  910 CONTINUE
1435  END IF
1436  920 CONTINUE
1437 *
1438  DO 940 k = 1, kb - 1
1439  j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1440 *
1441 * finish applying rotations in 1st set from the left
1442 *
1443  DO 930 l = kb - k, 1, -1
1444  nrt = ( j2+l-1 ) / ka1
1445  j1t = j2 - ( nrt-1 )*ka1
1446  IF( nrt.GT.0 )
1447  $ CALL clartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,
1448  $ ab( ka1-l, j1t-ka1+l ), inca,
1449  $ rwork( j1t ), work( j1t ), ka1 )
1450  930 CONTINUE
1451  940 CONTINUE
1452 *
1453  IF( kb.GT.1 ) THEN
1454  DO 950 j = 2, i2 - ka
1455  rwork( j ) = rwork( j+ka )
1456  work( j ) = work( j+ka )
1457  950 CONTINUE
1458  END IF
1459 *
1460  END IF
1461 *
1462  GO TO 490
1463 *
1464 * End of CHBGST
1465 *
subroutine clartg(f, g, c, s, r)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition: clartg.f90:118
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:78
subroutine cgerc(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
CGERC
Definition: cgerc.f:130
subroutine cgeru(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
CGERU
Definition: cgeru.f:130
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:74
subroutine clartv(N, X, INCX, Y, INCY, C, S, INCC)
CLARTV applies a vector of plane rotations with real cosines and complex sines to the elements of a p...
Definition: clartv.f:107
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: claset.f:106
subroutine crot(N, CX, INCX, CY, INCY, C, S)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition: crot.f:103
subroutine clargv(N, X, INCX, Y, INCY, C, INCC)
CLARGV generates a vector of plane rotations with real cosines and complex sines.
Definition: clargv.f:122
subroutine clar2v(N, X, Y, Z, INCX, C, S, INCC)
CLAR2V applies a vector of plane rotations with real cosines and complex sines from both sides to a s...
Definition: clar2v.f:111
Here is the call graph for this function:
Here is the caller graph for this function: