LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zgghd3 ( character  COMPQ,
character  COMPZ,
integer  N,
integer  ILO,
integer  IHI,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( ldq, * )  Q,
integer  LDQ,
complex*16, dimension( ldz, * )  Z,
integer  LDZ,
complex*16, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

ZGGHD3

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

Purpose:
 ZGGHD3 reduces a pair of complex matrices (A,B) to generalized upper
 Hessenberg form using unitary transformations, where A is a
 general matrix and B is upper triangular.  The form of the
 generalized eigenvalue problem is
    A*x = lambda*B*x,
 and B is typically made upper triangular by computing its QR
 factorization and moving the unitary matrix Q to the left side
 of the equation.

 This subroutine simultaneously reduces A to a Hessenberg matrix H:
    Q**H*A*Z = H
 and transforms B to another upper triangular matrix T:
    Q**H*B*Z = T
 in order to reduce the problem to its standard form
    H*y = lambda*T*y
 where y = Z**H*x.

 The unitary matrices Q and Z are determined as products of Givens
 rotations.  They may either be formed explicitly, or they may be
 postmultiplied into input matrices Q1 and Z1, so that
      Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H
      Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H
 If Q1 is the unitary matrix from the QR factorization of B in the
 original equation A*x = lambda*B*x, then ZGGHD3 reduces the original
 problem to generalized Hessenberg form.

 This is a blocked variant of CGGHRD, using matrix-matrix
 multiplications for parts of the computation to enhance performance.
Parameters
[in]COMPQ
          COMPQ is CHARACTER*1
          = 'N': do not compute Q;
          = 'I': Q is initialized to the unit matrix, and the
                 unitary matrix Q is returned;
          = 'V': Q must contain a unitary matrix Q1 on entry,
                 and the product Q1*Q is returned.
[in]COMPZ
          COMPZ is CHARACTER*1
          = 'N': do not compute Z;
          = 'I': Z is initialized to the unit matrix, and the
                 unitary matrix Z is returned;
          = 'V': Z must contain a unitary matrix Z1 on entry,
                 and the product Z1*Z is returned.
[in]N
          N is INTEGER
          The order of the matrices A and B.  N >= 0.
[in]ILO
          ILO is INTEGER
[in]IHI
          IHI is INTEGER

          ILO and IHI mark the rows and columns of A which are to be
          reduced.  It is assumed that A is already upper triangular
          in rows and columns 1:ILO-1 and IHI+1:N.  ILO and IHI are
          normally set by a previous call to ZGGBAL; otherwise they
          should be set to 1 and N respectively.
          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA, N)
          On entry, the N-by-N general matrix to be reduced.
          On exit, the upper triangle and the first subdiagonal of A
          are overwritten with the upper Hessenberg matrix H, and the
          rest is set to zero.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in,out]B
          B is COMPLEX*16 array, dimension (LDB, N)
          On entry, the N-by-N upper triangular matrix B.
          On exit, the upper triangular matrix T = Q**H B Z.  The
          elements below the diagonal are set to zero.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max(1,N).
[in,out]Q
          Q is COMPLEX*16 array, dimension (LDQ, N)
          On entry, if COMPQ = 'V', the unitary matrix Q1, typically
          from the QR factorization of B.
          On exit, if COMPQ='I', the unitary matrix Q, and if
          COMPQ = 'V', the product Q1*Q.
          Not referenced if COMPQ='N'.
[in]LDQ
          LDQ is INTEGER
          The leading dimension of the array Q.
          LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.
[in,out]Z
          Z is COMPLEX*16 array, dimension (LDZ, N)
          On entry, if COMPZ = 'V', the unitary matrix Z1.
          On exit, if COMPZ='I', the unitary matrix Z, and if
          COMPZ = 'V', the product Z1*Z.
          Not referenced if COMPZ='N'.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z.
          LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.
[out]WORK
          WORK is COMPLEX*16 array, dimension (LWORK)
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The length of the array WORK.  LWORK >= 1.
          For optimum performance LWORK >= 6*N*NB, where NB is the
          optimal blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[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.
Date
January 2015
Further Details:
  This routine reduces A to Hessenberg form and maintains B in
  using a blocked variant of Moler and Stewart's original algorithm,
  as described by Kagstrom, Kressner, Quintana-Orti, and Quintana-Orti
  (BIT 2008).

Definition at line 229 of file zgghd3.f.

229 *
230 * -- LAPACK computational routine (version 3.6.1) --
231 * -- LAPACK is a software package provided by Univ. of Tennessee, --
232 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
233 * January 2015
234 *
235  IMPLICIT NONE
236 *
237 * .. Scalar Arguments ..
238  CHARACTER compq, compz
239  INTEGER ihi, ilo, info, lda, ldb, ldq, ldz, n, lwork
240 * ..
241 * .. Array Arguments ..
242  COMPLEX*16 a( lda, * ), b( ldb, * ), q( ldq, * ),
243  $ z( ldz, * ), work( * )
244 * ..
245 *
246 * =====================================================================
247 *
248 * .. Parameters ..
249  COMPLEX*16 cone, czero
250  parameter ( cone = ( 1.0d+0, 0.0d+0 ),
251  $ czero = ( 0.0d+0, 0.0d+0 ) )
252 * ..
253 * .. Local Scalars ..
254  LOGICAL blk22, initq, initz, lquery, wantq, wantz
255  CHARACTER*1 compq2, compz2
256  INTEGER cola, i, ierr, j, j0, jcol, jj, jrow, k,
257  $ kacc22, len, lwkopt, n2nb, nb, nblst, nbmin,
258  $ nh, nnb, nx, ppw, ppwo, pw, top, topq
259  DOUBLE PRECISION c
260  COMPLEX*16 c1, c2, ctemp, s, s1, s2, temp, temp1, temp2,
261  $ temp3
262 * ..
263 * .. External Functions ..
264  LOGICAL lsame
265  INTEGER ilaenv
266  EXTERNAL ilaenv, lsame
267 * ..
268 * .. External Subroutines ..
269  EXTERNAL zgghrd, zlartg, zlaset, zunm22, zrot, xerbla
270 * ..
271 * .. Intrinsic Functions ..
272  INTRINSIC dble, dcmplx, dconjg, max
273 * ..
274 * .. Executable Statements ..
275 *
276 * Decode and test the input parameters.
277 *
278  info = 0
279  nb = ilaenv( 1, 'ZGGHD3', ' ', n, ilo, ihi, -1 )
280  lwkopt = max( 6*n*nb, 1 )
281  work( 1 ) = dcmplx( lwkopt )
282  initq = lsame( compq, 'I' )
283  wantq = initq .OR. lsame( compq, 'V' )
284  initz = lsame( compz, 'I' )
285  wantz = initz .OR. lsame( compz, 'V' )
286  lquery = ( lwork.EQ.-1 )
287 *
288  IF( .NOT.lsame( compq, 'N' ) .AND. .NOT.wantq ) THEN
289  info = -1
290  ELSE IF( .NOT.lsame( compz, 'N' ) .AND. .NOT.wantz ) THEN
291  info = -2
292  ELSE IF( n.LT.0 ) THEN
293  info = -3
294  ELSE IF( ilo.LT.1 ) THEN
295  info = -4
296  ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 ) THEN
297  info = -5
298  ELSE IF( lda.LT.max( 1, n ) ) THEN
299  info = -7
300  ELSE IF( ldb.LT.max( 1, n ) ) THEN
301  info = -9
302  ELSE IF( ( wantq .AND. ldq.LT.n ) .OR. ldq.LT.1 ) THEN
303  info = -11
304  ELSE IF( ( wantz .AND. ldz.LT.n ) .OR. ldz.LT.1 ) THEN
305  info = -13
306  ELSE IF( lwork.LT.1 .AND. .NOT.lquery ) THEN
307  info = -15
308  END IF
309  IF( info.NE.0 ) THEN
310  CALL xerbla( 'ZGGHD3', -info )
311  RETURN
312  ELSE IF( lquery ) THEN
313  RETURN
314  END IF
315 *
316 * Initialize Q and Z if desired.
317 *
318  IF( initq )
319  $ CALL zlaset( 'All', n, n, czero, cone, q, ldq )
320  IF( initz )
321  $ CALL zlaset( 'All', n, n, czero, cone, z, ldz )
322 *
323 * Zero out lower triangle of B.
324 *
325  IF( n.GT.1 )
326  $ CALL zlaset( 'Lower', n-1, n-1, czero, czero, b(2, 1), ldb )
327 *
328 * Quick return if possible
329 *
330  nh = ihi - ilo + 1
331  IF( nh.LE.1 ) THEN
332  work( 1 ) = cone
333  RETURN
334  END IF
335 *
336 * Determine the blocksize.
337 *
338  nbmin = ilaenv( 2, 'ZGGHD3', ' ', n, ilo, ihi, -1 )
339  IF( nb.GT.1 .AND. nb.LT.nh ) THEN
340 *
341 * Determine when to use unblocked instead of blocked code.
342 *
343  nx = max( nb, ilaenv( 3, 'ZGGHD3', ' ', n, ilo, ihi, -1 ) )
344  IF( nx.LT.nh ) THEN
345 *
346 * Determine if workspace is large enough for blocked code.
347 *
348  IF( lwork.LT.lwkopt ) THEN
349 *
350 * Not enough workspace to use optimal NB: determine the
351 * minimum value of NB, and reduce NB or force use of
352 * unblocked code.
353 *
354  nbmin = max( 2, ilaenv( 2, 'ZGGHD3', ' ', n, ilo, ihi,
355  $ -1 ) )
356  IF( lwork.GE.6*n*nbmin ) THEN
357  nb = lwork / ( 6*n )
358  ELSE
359  nb = 1
360  END IF
361  END IF
362  END IF
363  END IF
364 *
365  IF( nb.LT.nbmin .OR. nb.GE.nh ) THEN
366 *
367 * Use unblocked code below
368 *
369  jcol = ilo
370 *
371  ELSE
372 *
373 * Use blocked code
374 *
375  kacc22 = ilaenv( 16, 'ZGGHD3', ' ', n, ilo, ihi, -1 )
376  blk22 = kacc22.EQ.2
377  DO jcol = ilo, ihi-2, nb
378  nnb = min( nb, ihi-jcol-1 )
379 *
380 * Initialize small unitary factors that will hold the
381 * accumulated Givens rotations in workspace.
382 * N2NB denotes the number of 2*NNB-by-2*NNB factors
383 * NBLST denotes the (possibly smaller) order of the last
384 * factor.
385 *
386  n2nb = ( ihi-jcol-1 ) / nnb - 1
387  nblst = ihi - jcol - n2nb*nnb
388  CALL zlaset( 'All', nblst, nblst, czero, cone, work, nblst )
389  pw = nblst * nblst + 1
390  DO i = 1, n2nb
391  CALL zlaset( 'All', 2*nnb, 2*nnb, czero, cone,
392  $ work( pw ), 2*nnb )
393  pw = pw + 4*nnb*nnb
394  END DO
395 *
396 * Reduce columns JCOL:JCOL+NNB-1 of A to Hessenberg form.
397 *
398  DO j = jcol, jcol+nnb-1
399 *
400 * Reduce Jth column of A. Store cosines and sines in Jth
401 * column of A and B, respectively.
402 *
403  DO i = ihi, j+2, -1
404  temp = a( i-1, j )
405  CALL zlartg( temp, a( i, j ), c, s, a( i-1, j ) )
406  a( i, j ) = dcmplx( c )
407  b( i, j ) = s
408  END DO
409 *
410 * Accumulate Givens rotations into workspace array.
411 *
412  ppw = ( nblst + 1 )*( nblst - 2 ) - j + jcol + 1
413  len = 2 + j - jcol
414  jrow = j + n2nb*nnb + 2
415  DO i = ihi, jrow, -1
416  ctemp = a( i, j )
417  s = b( i, j )
418  DO jj = ppw, ppw+len-1
419  temp = work( jj + nblst )
420  work( jj + nblst ) = ctemp*temp - s*work( jj )
421  work( jj ) = dconjg( s )*temp + ctemp*work( jj )
422  END DO
423  len = len + 1
424  ppw = ppw - nblst - 1
425  END DO
426 *
427  ppwo = nblst*nblst + ( nnb+j-jcol-1 )*2*nnb + nnb
428  j0 = jrow - nnb
429  DO jrow = j0, j+2, -nnb
430  ppw = ppwo
431  len = 2 + j - jcol
432  DO i = jrow+nnb-1, jrow, -1
433  ctemp = a( i, j )
434  s = b( i, j )
435  DO jj = ppw, ppw+len-1
436  temp = work( jj + 2*nnb )
437  work( jj + 2*nnb ) = ctemp*temp - s*work( jj )
438  work( jj ) = dconjg( s )*temp + ctemp*work( jj )
439  END DO
440  len = len + 1
441  ppw = ppw - 2*nnb - 1
442  END DO
443  ppwo = ppwo + 4*nnb*nnb
444  END DO
445 *
446 * TOP denotes the number of top rows in A and B that will
447 * not be updated during the next steps.
448 *
449  IF( jcol.LE.2 ) THEN
450  top = 0
451  ELSE
452  top = jcol
453  END IF
454 *
455 * Propagate transformations through B and replace stored
456 * left sines/cosines by right sines/cosines.
457 *
458  DO jj = n, j+1, -1
459 *
460 * Update JJth column of B.
461 *
462  DO i = min( jj+1, ihi ), j+2, -1
463  ctemp = a( i, j )
464  s = b( i, j )
465  temp = b( i, jj )
466  b( i, jj ) = ctemp*temp - dconjg( s )*b( i-1, jj )
467  b( i-1, jj ) = s*temp + ctemp*b( i-1, jj )
468  END DO
469 *
470 * Annihilate B( JJ+1, JJ ).
471 *
472  IF( jj.LT.ihi ) THEN
473  temp = b( jj+1, jj+1 )
474  CALL zlartg( temp, b( jj+1, jj ), c, s,
475  $ b( jj+1, jj+1 ) )
476  b( jj+1, jj ) = czero
477  CALL zrot( jj-top, b( top+1, jj+1 ), 1,
478  $ b( top+1, jj ), 1, c, s )
479  a( jj+1, j ) = dcmplx( c )
480  b( jj+1, j ) = -dconjg( s )
481  END IF
482  END DO
483 *
484 * Update A by transformations from right.
485 *
486  jj = mod( ihi-j-1, 3 )
487  DO i = ihi-j-3, jj+1, -3
488  ctemp = a( j+1+i, j )
489  s = -b( j+1+i, j )
490  c1 = a( j+2+i, j )
491  s1 = -b( j+2+i, j )
492  c2 = a( j+3+i, j )
493  s2 = -b( j+3+i, j )
494 *
495  DO k = top+1, ihi
496  temp = a( k, j+i )
497  temp1 = a( k, j+i+1 )
498  temp2 = a( k, j+i+2 )
499  temp3 = a( k, j+i+3 )
500  a( k, j+i+3 ) = c2*temp3 + dconjg( s2 )*temp2
501  temp2 = -s2*temp3 + c2*temp2
502  a( k, j+i+2 ) = c1*temp2 + dconjg( s1 )*temp1
503  temp1 = -s1*temp2 + c1*temp1
504  a( k, j+i+1 ) = ctemp*temp1 + dconjg( s )*temp
505  a( k, j+i ) = -s*temp1 + ctemp*temp
506  END DO
507  END DO
508 *
509  IF( jj.GT.0 ) THEN
510  DO i = jj, 1, -1
511  c = dble( a( j+1+i, j ) )
512  CALL zrot( ihi-top, a( top+1, j+i+1 ), 1,
513  $ a( top+1, j+i ), 1, c,
514  $ -dconjg( b( j+1+i, j ) ) )
515  END DO
516  END IF
517 *
518 * Update (J+1)th column of A by transformations from left.
519 *
520  IF ( j .LT. jcol + nnb - 1 ) THEN
521  len = 1 + j - jcol
522 *
523 * Multiply with the trailing accumulated unitary
524 * matrix, which takes the form
525 *
526 * [ U11 U12 ]
527 * U = [ ],
528 * [ U21 U22 ]
529 *
530 * where U21 is a LEN-by-LEN matrix and U12 is lower
531 * triangular.
532 *
533  jrow = ihi - nblst + 1
534  CALL zgemv( 'Conjugate', nblst, len, cone, work,
535  $ nblst, a( jrow, j+1 ), 1, czero,
536  $ work( pw ), 1 )
537  ppw = pw + len
538  DO i = jrow, jrow+nblst-len-1
539  work( ppw ) = a( i, j+1 )
540  ppw = ppw + 1
541  END DO
542  CALL ztrmv( 'Lower', 'Conjugate', 'Non-unit',
543  $ nblst-len, work( len*nblst + 1 ), nblst,
544  $ work( pw+len ), 1 )
545  CALL zgemv( 'Conjugate', len, nblst-len, cone,
546  $ work( (len+1)*nblst - len + 1 ), nblst,
547  $ a( jrow+nblst-len, j+1 ), 1, cone,
548  $ work( pw+len ), 1 )
549  ppw = pw
550  DO i = jrow, jrow+nblst-1
551  a( i, j+1 ) = work( ppw )
552  ppw = ppw + 1
553  END DO
554 *
555 * Multiply with the other accumulated unitary
556 * matrices, which take the form
557 *
558 * [ U11 U12 0 ]
559 * [ ]
560 * U = [ U21 U22 0 ],
561 * [ ]
562 * [ 0 0 I ]
563 *
564 * where I denotes the (NNB-LEN)-by-(NNB-LEN) identity
565 * matrix, U21 is a LEN-by-LEN upper triangular matrix
566 * and U12 is an NNB-by-NNB lower triangular matrix.
567 *
568  ppwo = 1 + nblst*nblst
569  j0 = jrow - nnb
570  DO jrow = j0, jcol+1, -nnb
571  ppw = pw + len
572  DO i = jrow, jrow+nnb-1
573  work( ppw ) = a( i, j+1 )
574  ppw = ppw + 1
575  END DO
576  ppw = pw
577  DO i = jrow+nnb, jrow+nnb+len-1
578  work( ppw ) = a( i, j+1 )
579  ppw = ppw + 1
580  END DO
581  CALL ztrmv( 'Upper', 'Conjugate', 'Non-unit', len,
582  $ work( ppwo + nnb ), 2*nnb, work( pw ),
583  $ 1 )
584  CALL ztrmv( 'Lower', 'Conjugate', 'Non-unit', nnb,
585  $ work( ppwo + 2*len*nnb ),
586  $ 2*nnb, work( pw + len ), 1 )
587  CALL zgemv( 'Conjugate', nnb, len, cone,
588  $ work( ppwo ), 2*nnb, a( jrow, j+1 ), 1,
589  $ cone, work( pw ), 1 )
590  CALL zgemv( 'Conjugate', len, nnb, cone,
591  $ work( ppwo + 2*len*nnb + nnb ), 2*nnb,
592  $ a( jrow+nnb, j+1 ), 1, cone,
593  $ work( pw+len ), 1 )
594  ppw = pw
595  DO i = jrow, jrow+len+nnb-1
596  a( i, j+1 ) = work( ppw )
597  ppw = ppw + 1
598  END DO
599  ppwo = ppwo + 4*nnb*nnb
600  END DO
601  END IF
602  END DO
603 *
604 * Apply accumulated unitary matrices to A.
605 *
606  cola = n - jcol - nnb + 1
607  j = ihi - nblst + 1
608  CALL zgemm( 'Conjugate', 'No Transpose', nblst,
609  $ cola, nblst, cone, work, nblst,
610  $ a( j, jcol+nnb ), lda, czero, work( pw ),
611  $ nblst )
612  CALL zlacpy( 'All', nblst, cola, work( pw ), nblst,
613  $ a( j, jcol+nnb ), lda )
614  ppwo = nblst*nblst + 1
615  j0 = j - nnb
616  DO j = j0, jcol+1, -nnb
617  IF ( blk22 ) THEN
618 *
619 * Exploit the structure of
620 *
621 * [ U11 U12 ]
622 * U = [ ]
623 * [ U21 U22 ],
624 *
625 * where all blocks are NNB-by-NNB, U21 is upper
626 * triangular and U12 is lower triangular.
627 *
628  CALL zunm22( 'Left', 'Conjugate', 2*nnb, cola, nnb,
629  $ nnb, work( ppwo ), 2*nnb,
630  $ a( j, jcol+nnb ), lda, work( pw ),
631  $ lwork-pw+1, ierr )
632  ELSE
633 *
634 * Ignore the structure of U.
635 *
636  CALL zgemm( 'Conjugate', 'No Transpose', 2*nnb,
637  $ cola, 2*nnb, cone, work( ppwo ), 2*nnb,
638  $ a( j, jcol+nnb ), lda, czero, work( pw ),
639  $ 2*nnb )
640  CALL zlacpy( 'All', 2*nnb, cola, work( pw ), 2*nnb,
641  $ a( j, jcol+nnb ), lda )
642  END IF
643  ppwo = ppwo + 4*nnb*nnb
644  END DO
645 *
646 * Apply accumulated unitary matrices to Q.
647 *
648  IF( wantq ) THEN
649  j = ihi - nblst + 1
650  IF ( initq ) THEN
651  topq = max( 2, j - jcol + 1 )
652  nh = ihi - topq + 1
653  ELSE
654  topq = 1
655  nh = n
656  END IF
657  CALL zgemm( 'No Transpose', 'No Transpose', nh,
658  $ nblst, nblst, cone, q( topq, j ), ldq,
659  $ work, nblst, czero, work( pw ), nh )
660  CALL zlacpy( 'All', nh, nblst, work( pw ), nh,
661  $ q( topq, j ), ldq )
662  ppwo = nblst*nblst + 1
663  j0 = j - nnb
664  DO j = j0, jcol+1, -nnb
665  IF ( initq ) THEN
666  topq = max( 2, j - jcol + 1 )
667  nh = ihi - topq + 1
668  END IF
669  IF ( blk22 ) THEN
670 *
671 * Exploit the structure of U.
672 *
673  CALL zunm22( 'Right', 'No Transpose', nh, 2*nnb,
674  $ nnb, nnb, work( ppwo ), 2*nnb,
675  $ q( topq, j ), ldq, work( pw ),
676  $ lwork-pw+1, ierr )
677  ELSE
678 *
679 * Ignore the structure of U.
680 *
681  CALL zgemm( 'No Transpose', 'No Transpose', nh,
682  $ 2*nnb, 2*nnb, cone, q( topq, j ), ldq,
683  $ work( ppwo ), 2*nnb, czero, work( pw ),
684  $ nh )
685  CALL zlacpy( 'All', nh, 2*nnb, work( pw ), nh,
686  $ q( topq, j ), ldq )
687  END IF
688  ppwo = ppwo + 4*nnb*nnb
689  END DO
690  END IF
691 *
692 * Accumulate right Givens rotations if required.
693 *
694  IF ( wantz .OR. top.GT.0 ) THEN
695 *
696 * Initialize small unitary factors that will hold the
697 * accumulated Givens rotations in workspace.
698 *
699  CALL zlaset( 'All', nblst, nblst, czero, cone, work,
700  $ nblst )
701  pw = nblst * nblst + 1
702  DO i = 1, n2nb
703  CALL zlaset( 'All', 2*nnb, 2*nnb, czero, cone,
704  $ work( pw ), 2*nnb )
705  pw = pw + 4*nnb*nnb
706  END DO
707 *
708 * Accumulate Givens rotations into workspace array.
709 *
710  DO j = jcol, jcol+nnb-1
711  ppw = ( nblst + 1 )*( nblst - 2 ) - j + jcol + 1
712  len = 2 + j - jcol
713  jrow = j + n2nb*nnb + 2
714  DO i = ihi, jrow, -1
715  ctemp = a( i, j )
716  a( i, j ) = czero
717  s = b( i, j )
718  b( i, j ) = czero
719  DO jj = ppw, ppw+len-1
720  temp = work( jj + nblst )
721  work( jj + nblst ) = ctemp*temp -
722  $ dconjg( s )*work( jj )
723  work( jj ) = s*temp + ctemp*work( jj )
724  END DO
725  len = len + 1
726  ppw = ppw - nblst - 1
727  END DO
728 *
729  ppwo = nblst*nblst + ( nnb+j-jcol-1 )*2*nnb + nnb
730  j0 = jrow - nnb
731  DO jrow = j0, j+2, -nnb
732  ppw = ppwo
733  len = 2 + j - jcol
734  DO i = jrow+nnb-1, jrow, -1
735  ctemp = a( i, j )
736  a( i, j ) = czero
737  s = b( i, j )
738  b( i, j ) = czero
739  DO jj = ppw, ppw+len-1
740  temp = work( jj + 2*nnb )
741  work( jj + 2*nnb ) = ctemp*temp -
742  $ dconjg( s )*work( jj )
743  work( jj ) = s*temp + ctemp*work( jj )
744  END DO
745  len = len + 1
746  ppw = ppw - 2*nnb - 1
747  END DO
748  ppwo = ppwo + 4*nnb*nnb
749  END DO
750  END DO
751  ELSE
752 *
753  CALL zlaset( 'Lower', ihi - jcol - 1, nnb, czero, czero,
754  $ a( jcol + 2, jcol ), lda )
755  CALL zlaset( 'Lower', ihi - jcol - 1, nnb, czero, czero,
756  $ b( jcol + 2, jcol ), ldb )
757  END IF
758 *
759 * Apply accumulated unitary matrices to A and B.
760 *
761  IF ( top.GT.0 ) THEN
762  j = ihi - nblst + 1
763  CALL zgemm( 'No Transpose', 'No Transpose', top,
764  $ nblst, nblst, cone, a( 1, j ), lda,
765  $ work, nblst, czero, work( pw ), top )
766  CALL zlacpy( 'All', top, nblst, work( pw ), top,
767  $ a( 1, j ), lda )
768  ppwo = nblst*nblst + 1
769  j0 = j - nnb
770  DO j = j0, jcol+1, -nnb
771  IF ( blk22 ) THEN
772 *
773 * Exploit the structure of U.
774 *
775  CALL zunm22( 'Right', 'No Transpose', top, 2*nnb,
776  $ nnb, nnb, work( ppwo ), 2*nnb,
777  $ a( 1, j ), lda, work( pw ),
778  $ lwork-pw+1, ierr )
779  ELSE
780 *
781 * Ignore the structure of U.
782 *
783  CALL zgemm( 'No Transpose', 'No Transpose', top,
784  $ 2*nnb, 2*nnb, cone, a( 1, j ), lda,
785  $ work( ppwo ), 2*nnb, czero,
786  $ work( pw ), top )
787  CALL zlacpy( 'All', top, 2*nnb, work( pw ), top,
788  $ a( 1, j ), lda )
789  END IF
790  ppwo = ppwo + 4*nnb*nnb
791  END DO
792 *
793  j = ihi - nblst + 1
794  CALL zgemm( 'No Transpose', 'No Transpose', top,
795  $ nblst, nblst, cone, b( 1, j ), ldb,
796  $ work, nblst, czero, work( pw ), top )
797  CALL zlacpy( 'All', top, nblst, work( pw ), top,
798  $ b( 1, j ), ldb )
799  ppwo = nblst*nblst + 1
800  j0 = j - nnb
801  DO j = j0, jcol+1, -nnb
802  IF ( blk22 ) THEN
803 *
804 * Exploit the structure of U.
805 *
806  CALL zunm22( 'Right', 'No Transpose', top, 2*nnb,
807  $ nnb, nnb, work( ppwo ), 2*nnb,
808  $ b( 1, j ), ldb, work( pw ),
809  $ lwork-pw+1, ierr )
810  ELSE
811 *
812 * Ignore the structure of U.
813 *
814  CALL zgemm( 'No Transpose', 'No Transpose', top,
815  $ 2*nnb, 2*nnb, cone, b( 1, j ), ldb,
816  $ work( ppwo ), 2*nnb, czero,
817  $ work( pw ), top )
818  CALL zlacpy( 'All', top, 2*nnb, work( pw ), top,
819  $ b( 1, j ), ldb )
820  END IF
821  ppwo = ppwo + 4*nnb*nnb
822  END DO
823  END IF
824 *
825 * Apply accumulated unitary matrices to Z.
826 *
827  IF( wantz ) THEN
828  j = ihi - nblst + 1
829  IF ( initq ) THEN
830  topq = max( 2, j - jcol + 1 )
831  nh = ihi - topq + 1
832  ELSE
833  topq = 1
834  nh = n
835  END IF
836  CALL zgemm( 'No Transpose', 'No Transpose', nh,
837  $ nblst, nblst, cone, z( topq, j ), ldz,
838  $ work, nblst, czero, work( pw ), nh )
839  CALL zlacpy( 'All', nh, nblst, work( pw ), nh,
840  $ z( topq, j ), ldz )
841  ppwo = nblst*nblst + 1
842  j0 = j - nnb
843  DO j = j0, jcol+1, -nnb
844  IF ( initq ) THEN
845  topq = max( 2, j - jcol + 1 )
846  nh = ihi - topq + 1
847  END IF
848  IF ( blk22 ) THEN
849 *
850 * Exploit the structure of U.
851 *
852  CALL zunm22( 'Right', 'No Transpose', nh, 2*nnb,
853  $ nnb, nnb, work( ppwo ), 2*nnb,
854  $ z( topq, j ), ldz, work( pw ),
855  $ lwork-pw+1, ierr )
856  ELSE
857 *
858 * Ignore the structure of U.
859 *
860  CALL zgemm( 'No Transpose', 'No Transpose', nh,
861  $ 2*nnb, 2*nnb, cone, z( topq, j ), ldz,
862  $ work( ppwo ), 2*nnb, czero, work( pw ),
863  $ nh )
864  CALL zlacpy( 'All', nh, 2*nnb, work( pw ), nh,
865  $ z( topq, j ), ldz )
866  END IF
867  ppwo = ppwo + 4*nnb*nnb
868  END DO
869  END IF
870  END DO
871  END IF
872 *
873 * Use unblocked code to reduce the rest of the matrix
874 * Avoid re-initialization of modified Q and Z.
875 *
876  compq2 = compq
877  compz2 = compz
878  IF ( jcol.NE.ilo ) THEN
879  IF ( wantq )
880  $ compq2 = 'V'
881  IF ( wantz )
882  $ compz2 = 'V'
883  END IF
884 *
885  IF ( jcol.LT.ihi )
886  $ CALL zgghrd( compq2, compz2, n, jcol, ihi, a, lda, b, ldb, q,
887  $ ldq, z, ldz, ierr )
888  work( 1 ) = dcmplx( lwkopt )
889 *
890  RETURN
891 *
892 * End of ZGGHD3
893 *
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine zgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
ZGGHRD
Definition: zgghrd.f:206
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zunm22(SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC, WORK, LWORK, INFO)
ZUNM22 multiplies a general matrix by a banded unitary matrix.
Definition: zunm22.f:164
subroutine ztrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
ZTRMV
Definition: ztrmv.f:149
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine zrot(N, CX, INCX, CY, INCY, C, S)
ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors...
Definition: zrot.f:105
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zlartg(F, G, CS, SN, R)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition: zlartg.f:105

Here is the call graph for this function:

Here is the caller graph for this function: