LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ cgghd3()

subroutine cgghd3 ( character  COMPQ,
character  COMPZ,
integer  N,
integer  ILO,
integer  IHI,
complex, dimension( lda, * )  A,
integer  LDA,
complex, dimension( ldb, * )  B,
integer  LDB,
complex, dimension( ldq, * )  Q,
integer  LDQ,
complex, dimension( ldz, * )  Z,
integer  LDZ,
complex, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

CGGHD3

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

Purpose:
 CGGHD3 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 CGGHD3 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 CGGBAL; 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 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 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 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 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 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 229 of file cgghd3.f.

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