LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ zhbgst()

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

ZHBGST

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

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