LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
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 233 of file cgghd3.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: