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

◆ zgghd3()

subroutine zgghd3 ( character  compq,
character  compz,
integer  n,
integer  ilo,
integer  ihi,
complex*16, dimension( lda, * )  a,
integer  lda,
complex*16, dimension( ldb, * )  b,
integer  ldb,
complex*16, dimension( ldq, * )  q,
integer  ldq,
complex*16, dimension( ldz, * )  z,
integer  ldz,
complex*16, dimension( * )  work,
integer  lwork,
integer  info 
)

ZGGHD3

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

Purpose:
 ZGGHD3 reduces a pair of complex matrices (A,B) to generalized upper
 Hessenberg form using unitary transformations, where A is a
 general matrix and B is upper triangular.  The form of the
 generalized eigenvalue problem is
    A*x = lambda*B*x,
 and B is typically made upper triangular by computing its QR
 factorization and moving the unitary matrix Q to the left side
 of the equation.

 This subroutine simultaneously reduces A to a Hessenberg matrix H:
    Q**H*A*Z = H
 and transforms B to another upper triangular matrix T:
    Q**H*B*Z = T
 in order to reduce the problem to its standard form
    H*y = lambda*T*y
 where y = Z**H*x.

 The unitary matrices Q and Z are determined as products of Givens
 rotations.  They may either be formed explicitly, or they may be
 postmultiplied into input matrices Q1 and Z1, so that
      Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H
      Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H
 If Q1 is the unitary matrix from the QR factorization of B in the
 original equation A*x = lambda*B*x, then ZGGHD3 reduces the original
 problem to generalized Hessenberg form.

 This is a blocked variant of CGGHRD, using matrix-matrix
 multiplications for parts of the computation to enhance performance.
Parameters
[in]COMPQ
          COMPQ is CHARACTER*1
          = 'N': do not compute Q;
          = 'I': Q is initialized to the unit matrix, and the
                 unitary matrix Q is returned;
          = 'V': Q must contain a unitary matrix Q1 on entry,
                 and the product Q1*Q is returned.
[in]COMPZ
          COMPZ is CHARACTER*1
          = 'N': do not compute Z;
          = 'I': Z is initialized to the unit matrix, and the
                 unitary matrix Z is returned;
          = 'V': Z must contain a unitary matrix Z1 on entry,
                 and the product Z1*Z is returned.
[in]N
          N is INTEGER
          The order of the matrices A and B.  N >= 0.
[in]ILO
          ILO is INTEGER
[in]IHI
          IHI is INTEGER

          ILO and IHI mark the rows and columns of A which are to be
          reduced.  It is assumed that A is already upper triangular
          in rows and columns 1:ILO-1 and IHI+1:N.  ILO and IHI are
          normally set by a previous call to ZGGBAL; otherwise they
          should be set to 1 and N respectively.
          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA, N)
          On entry, the N-by-N general matrix to be reduced.
          On exit, the upper triangle and the first subdiagonal of A
          are overwritten with the upper Hessenberg matrix H, and the
          rest is set to zero.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in,out]B
          B is COMPLEX*16 array, dimension (LDB, N)
          On entry, the N-by-N upper triangular matrix B.
          On exit, the upper triangular matrix T = Q**H B Z.  The
          elements below the diagonal are set to zero.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max(1,N).
[in,out]Q
          Q is COMPLEX*16 array, dimension (LDQ, N)
          On entry, if COMPQ = 'V', the unitary matrix Q1, typically
          from the QR factorization of B.
          On exit, if COMPQ='I', the unitary matrix Q, and if
          COMPQ = 'V', the product Q1*Q.
          Not referenced if COMPQ='N'.
[in]LDQ
          LDQ is INTEGER
          The leading dimension of the array Q.
          LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.
[in,out]Z
          Z is COMPLEX*16 array, dimension (LDZ, N)
          On entry, if COMPZ = 'V', the unitary matrix Z1.
          On exit, if COMPZ='I', the unitary matrix Z, and if
          COMPZ = 'V', the product Z1*Z.
          Not referenced if COMPZ='N'.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z.
          LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.
[out]WORK
          WORK is COMPLEX*16 array, dimension (LWORK)
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The length of the array WORK.  LWORK >= 1.
          For optimum performance LWORK >= 6*N*NB, where NB is the
          optimal blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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 225 of file zgghd3.f.

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