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