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