LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ zgghd3()

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.
Further Details:
  This routine reduces A to Hessenberg form and maintains B in triangular form
  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 225 of file zgghd3.f.

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