LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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 = real( 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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
cgghrd
subroutine cgghrd(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, info)
CGGHRD
Definition cgghrd.f:204
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
clacpy
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
clartg
subroutine clartg(f, g, c, s, r)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition clartg.f90:116
claset
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
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
crot
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
ctrmv
subroutine ctrmv(uplo, trans, diag, n, a, lda, x, incx)
CTRMV
Definition ctrmv.f:147
cunm22
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
Here is the call graph for this function:
Here is the caller graph for this function: