LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ zhgeqz()

subroutine zhgeqz ( character  JOB,
character  COMPQ,
character  COMPZ,
integer  N,
integer  ILO,
integer  IHI,
complex*16, dimension( ldh, * )  H,
integer  LDH,
complex*16, dimension( ldt, * )  T,
integer  LDT,
complex*16, dimension( * )  ALPHA,
complex*16, dimension( * )  BETA,
complex*16, dimension( ldq, * )  Q,
integer  LDQ,
complex*16, dimension( ldz, * )  Z,
integer  LDZ,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZHGEQZ

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

Purpose:
 ZHGEQZ computes the eigenvalues of a complex matrix pair (H,T),
 where H is an upper Hessenberg matrix and T is upper triangular,
 using the single-shift QZ method.
 Matrix pairs of this type are produced by the reduction to
 generalized upper Hessenberg form of a complex matrix pair (A,B):

    A = Q1*H*Z1**H,  B = Q1*T*Z1**H,

 as computed by ZGGHRD.

 If JOB='S', then the Hessenberg-triangular pair (H,T) is
 also reduced to generalized Schur form,

    H = Q*S*Z**H,  T = Q*P*Z**H,

 where Q and Z are unitary matrices and S and P are upper triangular.

 Optionally, the unitary matrix Q from the generalized Schur
 factorization may be postmultiplied into an input matrix Q1, and the
 unitary matrix Z may be postmultiplied into an input matrix Z1.
 If Q1 and Z1 are the unitary matrices from ZGGHRD that reduced
 the matrix pair (A,B) to generalized Hessenberg form, then the output
 matrices Q1*Q and Z1*Z are the unitary factors from the generalized
 Schur factorization of (A,B):

    A = (Q1*Q)*S*(Z1*Z)**H,  B = (Q1*Q)*P*(Z1*Z)**H.

 To avoid overflow, eigenvalues of the matrix pair (H,T)
 (equivalently, of (A,B)) are computed as a pair of complex values
 (alpha,beta).  If beta is nonzero, lambda = alpha / beta is an
 eigenvalue of the generalized nonsymmetric eigenvalue problem (GNEP)
    A*x = lambda*B*x
 and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the
 alternate form of the GNEP
    mu*A*y = B*y.
 The values of alpha and beta for the i-th eigenvalue can be read
 directly from the generalized Schur form:  alpha = S(i,i),
 beta = P(i,i).

 Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
      Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
      pp. 241--256.
Parameters
[in]JOB
          JOB is CHARACTER*1
          = 'E': Compute eigenvalues only;
          = 'S': Computer eigenvalues and the Schur form.
[in]COMPQ
          COMPQ is CHARACTER*1
          = 'N': Left Schur vectors (Q) are not computed;
          = 'I': Q is initialized to the unit matrix and the matrix Q
                 of left Schur vectors of (H,T) 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': Right Schur vectors (Z) are not computed;
          = 'I': Q is initialized to the unit matrix and the matrix Z
                 of right Schur vectors of (H,T) 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 H, T, Q, and Z.  N >= 0.
[in]ILO
          ILO is INTEGER
[in]IHI
          IHI is INTEGER
          ILO and IHI mark the rows and columns of H which are in
          Hessenberg form.  It is assumed that A is already upper
          triangular in rows and columns 1:ILO-1 and IHI+1:N.
          If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0.
[in,out]H
          H is COMPLEX*16 array, dimension (LDH, N)
          On entry, the N-by-N upper Hessenberg matrix H.
          On exit, if JOB = 'S', H contains the upper triangular
          matrix S from the generalized Schur factorization.
          If JOB = 'E', the diagonal of H matches that of S, but
          the rest of H is unspecified.
[in]LDH
          LDH is INTEGER
          The leading dimension of the array H.  LDH >= max( 1, N ).
[in,out]T
          T is COMPLEX*16 array, dimension (LDT, N)
          On entry, the N-by-N upper triangular matrix T.
          On exit, if JOB = 'S', T contains the upper triangular
          matrix P from the generalized Schur factorization.
          If JOB = 'E', the diagonal of T matches that of P, but
          the rest of T is unspecified.
[in]LDT
          LDT is INTEGER
          The leading dimension of the array T.  LDT >= max( 1, N ).
[out]ALPHA
          ALPHA is COMPLEX*16 array, dimension (N)
          The complex scalars alpha that define the eigenvalues of
          GNEP.  ALPHA(i) = S(i,i) in the generalized Schur
          factorization.
[out]BETA
          BETA is COMPLEX*16 array, dimension (N)
          The real non-negative scalars beta that define the
          eigenvalues of GNEP.  BETA(i) = P(i,i) in the generalized
          Schur factorization.

          Together, the quantities alpha = ALPHA(j) and beta = BETA(j)
          represent the j-th eigenvalue of the matrix pair (A,B), in
          one of the forms lambda = alpha/beta or mu = beta/alpha.
          Since either lambda or mu may overflow, they should not,
          in general, be computed.
[in,out]Q
          Q is COMPLEX*16 array, dimension (LDQ, N)
          On entry, if COMPQ = 'V', the unitary matrix Q1 used in the
          reduction of (A,B) to generalized Hessenberg form.
          On exit, if COMPQ = 'I', the unitary matrix of left Schur
          vectors of (H,T), and if COMPQ = 'V', the unitary matrix of
          left Schur vectors of (A,B).
          Not referenced if COMPQ = 'N'.
[in]LDQ
          LDQ is INTEGER
          The leading dimension of the array Q.  LDQ >= 1.
          If COMPQ='V' or 'I', then LDQ >= N.
[in,out]Z
          Z is COMPLEX*16 array, dimension (LDZ, N)
          On entry, if COMPZ = 'V', the unitary matrix Z1 used in the
          reduction of (A,B) to generalized Hessenberg form.
          On exit, if COMPZ = 'I', the unitary matrix of right Schur
          vectors of (H,T), and if COMPZ = 'V', the unitary matrix of
          right Schur vectors of (A,B).
          Not referenced if COMPZ = 'N'.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z.  LDZ >= 1.
          If COMPZ='V' or 'I', then LDZ >= N.
[out]WORK
          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
          On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.  LWORK >= max(1,N).

          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]RWORK
          RWORK is DOUBLE PRECISION array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
          = 1,...,N: the QZ iteration did not converge.  (H,T) is not
                     in Schur form, but ALPHA(i) and BETA(i),
                     i=INFO+1,...,N should be correct.
          = N+1,...,2*N: the shift calculation failed.  (H,T) is not
                     in Schur form, but ALPHA(i) and BETA(i),
                     i=INFO-N+1,...,N should be correct.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
  We assume that complex ABS works as long as its value is less than
  overflow.

Definition at line 281 of file zhgeqz.f.

284 *
285 * -- LAPACK computational routine --
286 * -- LAPACK is a software package provided by Univ. of Tennessee, --
287 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
288 *
289 * .. Scalar Arguments ..
290  CHARACTER COMPQ, COMPZ, JOB
291  INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
292 * ..
293 * .. Array Arguments ..
294  DOUBLE PRECISION RWORK( * )
295  COMPLEX*16 ALPHA( * ), BETA( * ), H( LDH, * ),
296  $ Q( LDQ, * ), T( LDT, * ), WORK( * ),
297  $ Z( LDZ, * )
298 * ..
299 *
300 * =====================================================================
301 *
302 * .. Parameters ..
303  COMPLEX*16 CZERO, CONE
304  parameter( czero = ( 0.0d+0, 0.0d+0 ),
305  $ cone = ( 1.0d+0, 0.0d+0 ) )
306  DOUBLE PRECISION ZERO, ONE
307  parameter( zero = 0.0d+0, one = 1.0d+0 )
308  DOUBLE PRECISION HALF
309  parameter( half = 0.5d+0 )
310 * ..
311 * .. Local Scalars ..
312  LOGICAL ILAZR2, ILAZRO, ILQ, ILSCHR, ILZ, LQUERY
313  INTEGER ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
314  $ ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER,
315  $ JR, MAXIT
316  DOUBLE PRECISION ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,
317  $ C, SAFMIN, TEMP, TEMP2, TEMPR, ULP
318  COMPLEX*16 ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2,
319  $ CTEMP3, ESHIFT, S, SHIFT, SIGNBC,
320  $ U12, X, ABI12, Y
321 * ..
322 * .. External Functions ..
323  COMPLEX*16 ZLADIV
324  LOGICAL LSAME
325  DOUBLE PRECISION DLAMCH, ZLANHS
326  EXTERNAL zladiv, lsame, dlamch, zlanhs
327 * ..
328 * .. External Subroutines ..
329  EXTERNAL xerbla, zlartg, zlaset, zrot, zscal
330 * ..
331 * .. Intrinsic Functions ..
332  INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, min,
333  $ sqrt
334 * ..
335 * .. Statement Functions ..
336  DOUBLE PRECISION ABS1
337 * ..
338 * .. Statement Function definitions ..
339  abs1( x ) = abs( dble( x ) ) + abs( dimag( x ) )
340 * ..
341 * .. Executable Statements ..
342 *
343 * Decode JOB, COMPQ, COMPZ
344 *
345  IF( lsame( job, 'E' ) ) THEN
346  ilschr = .false.
347  ischur = 1
348  ELSE IF( lsame( job, 'S' ) ) THEN
349  ilschr = .true.
350  ischur = 2
351  ELSE
352  ilschr = .true.
353  ischur = 0
354  END IF
355 *
356  IF( lsame( compq, 'N' ) ) THEN
357  ilq = .false.
358  icompq = 1
359  ELSE IF( lsame( compq, 'V' ) ) THEN
360  ilq = .true.
361  icompq = 2
362  ELSE IF( lsame( compq, 'I' ) ) THEN
363  ilq = .true.
364  icompq = 3
365  ELSE
366  ilq = .true.
367  icompq = 0
368  END IF
369 *
370  IF( lsame( compz, 'N' ) ) THEN
371  ilz = .false.
372  icompz = 1
373  ELSE IF( lsame( compz, 'V' ) ) THEN
374  ilz = .true.
375  icompz = 2
376  ELSE IF( lsame( compz, 'I' ) ) THEN
377  ilz = .true.
378  icompz = 3
379  ELSE
380  ilz = .true.
381  icompz = 0
382  END IF
383 *
384 * Check Argument Values
385 *
386  info = 0
387  work( 1 ) = max( 1, n )
388  lquery = ( lwork.EQ.-1 )
389  IF( ischur.EQ.0 ) THEN
390  info = -1
391  ELSE IF( icompq.EQ.0 ) THEN
392  info = -2
393  ELSE IF( icompz.EQ.0 ) THEN
394  info = -3
395  ELSE IF( n.LT.0 ) THEN
396  info = -4
397  ELSE IF( ilo.LT.1 ) THEN
398  info = -5
399  ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 ) THEN
400  info = -6
401  ELSE IF( ldh.LT.n ) THEN
402  info = -8
403  ELSE IF( ldt.LT.n ) THEN
404  info = -10
405  ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) ) THEN
406  info = -14
407  ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) ) THEN
408  info = -16
409  ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
410  info = -18
411  END IF
412  IF( info.NE.0 ) THEN
413  CALL xerbla( 'ZHGEQZ', -info )
414  RETURN
415  ELSE IF( lquery ) THEN
416  RETURN
417  END IF
418 *
419 * Quick return if possible
420 *
421 * WORK( 1 ) = CMPLX( 1 )
422  IF( n.LE.0 ) THEN
423  work( 1 ) = dcmplx( 1 )
424  RETURN
425  END IF
426 *
427 * Initialize Q and Z
428 *
429  IF( icompq.EQ.3 )
430  $ CALL zlaset( 'Full', n, n, czero, cone, q, ldq )
431  IF( icompz.EQ.3 )
432  $ CALL zlaset( 'Full', n, n, czero, cone, z, ldz )
433 *
434 * Machine Constants
435 *
436  in = ihi + 1 - ilo
437  safmin = dlamch( 'S' )
438  ulp = dlamch( 'E' )*dlamch( 'B' )
439  anorm = zlanhs( 'F', in, h( ilo, ilo ), ldh, rwork )
440  bnorm = zlanhs( 'F', in, t( ilo, ilo ), ldt, rwork )
441  atol = max( safmin, ulp*anorm )
442  btol = max( safmin, ulp*bnorm )
443  ascale = one / max( safmin, anorm )
444  bscale = one / max( safmin, bnorm )
445 *
446 *
447 * Set Eigenvalues IHI+1:N
448 *
449  DO 10 j = ihi + 1, n
450  absb = abs( t( j, j ) )
451  IF( absb.GT.safmin ) THEN
452  signbc = dconjg( t( j, j ) / absb )
453  t( j, j ) = absb
454  IF( ilschr ) THEN
455  CALL zscal( j-1, signbc, t( 1, j ), 1 )
456  CALL zscal( j, signbc, h( 1, j ), 1 )
457  ELSE
458  CALL zscal( 1, signbc, h( j, j ), 1 )
459  END IF
460  IF( ilz )
461  $ CALL zscal( n, signbc, z( 1, j ), 1 )
462  ELSE
463  t( j, j ) = czero
464  END IF
465  alpha( j ) = h( j, j )
466  beta( j ) = t( j, j )
467  10 CONTINUE
468 *
469 * If IHI < ILO, skip QZ steps
470 *
471  IF( ihi.LT.ilo )
472  $ GO TO 190
473 *
474 * MAIN QZ ITERATION LOOP
475 *
476 * Initialize dynamic indices
477 *
478 * Eigenvalues ILAST+1:N have been found.
479 * Column operations modify rows IFRSTM:whatever
480 * Row operations modify columns whatever:ILASTM
481 *
482 * If only eigenvalues are being computed, then
483 * IFRSTM is the row of the last splitting row above row ILAST;
484 * this is always at least ILO.
485 * IITER counts iterations since the last eigenvalue was found,
486 * to tell when to use an extraordinary shift.
487 * MAXIT is the maximum number of QZ sweeps allowed.
488 *
489  ilast = ihi
490  IF( ilschr ) THEN
491  ifrstm = 1
492  ilastm = n
493  ELSE
494  ifrstm = ilo
495  ilastm = ihi
496  END IF
497  iiter = 0
498  eshift = czero
499  maxit = 30*( ihi-ilo+1 )
500 *
501  DO 170 jiter = 1, maxit
502 *
503 * Check for too many iterations.
504 *
505  IF( jiter.GT.maxit )
506  $ GO TO 180
507 *
508 * Split the matrix if possible.
509 *
510 * Two tests:
511 * 1: H(j,j-1)=0 or j=ILO
512 * 2: T(j,j)=0
513 *
514 * Special case: j=ILAST
515 *
516  IF( ilast.EQ.ilo ) THEN
517  GO TO 60
518  ELSE
519  IF( abs1( h( ilast, ilast-1 ) ).LE.max( safmin, ulp*(
520  $ abs1( h( ilast, ilast ) ) + abs1( h( ilast-1, ilast-1 )
521  $ ) ) ) ) THEN
522  h( ilast, ilast-1 ) = czero
523  GO TO 60
524  END IF
525  END IF
526 *
527  IF( abs( t( ilast, ilast ) ).LE.max( safmin, ulp*(
528  $ abs( t( ilast - 1, ilast ) ) + abs( t( ilast-1, ilast-1 )
529  $ ) ) ) ) THEN
530  t( ilast, ilast ) = czero
531  GO TO 50
532  END IF
533 *
534 * General case: j<ILAST
535 *
536  DO 40 j = ilast - 1, ilo, -1
537 *
538 * Test 1: for H(j,j-1)=0 or j=ILO
539 *
540  IF( j.EQ.ilo ) THEN
541  ilazro = .true.
542  ELSE
543  IF( abs1( h( j, j-1 ) ).LE.max( safmin, ulp*(
544  $ abs1( h( j, j ) ) + abs1( h( j-1, j-1 ) )
545  $ ) ) ) THEN
546  h( j, j-1 ) = czero
547  ilazro = .true.
548  ELSE
549  ilazro = .false.
550  END IF
551  END IF
552 *
553 * Test 2: for T(j,j)=0
554 *
555  temp = abs( t( j, j + 1 ) )
556  IF ( j .GT. ilo )
557  $ temp = temp + abs( t( j - 1, j ) )
558  IF( abs( t( j, j ) ).LT.max( safmin,ulp*temp ) ) THEN
559  t( j, j ) = czero
560 *
561 * Test 1a: Check for 2 consecutive small subdiagonals in A
562 *
563  ilazr2 = .false.
564  IF( .NOT.ilazro ) THEN
565  IF( abs1( h( j, j-1 ) )*( ascale*abs1( h( j+1,
566  $ j ) ) ).LE.abs1( h( j, j ) )*( ascale*atol ) )
567  $ ilazr2 = .true.
568  END IF
569 *
570 * If both tests pass (1 & 2), i.e., the leading diagonal
571 * element of B in the block is zero, split a 1x1 block off
572 * at the top. (I.e., at the J-th row/column) The leading
573 * diagonal element of the remainder can also be zero, so
574 * this may have to be done repeatedly.
575 *
576  IF( ilazro .OR. ilazr2 ) THEN
577  DO 20 jch = j, ilast - 1
578  ctemp = h( jch, jch )
579  CALL zlartg( ctemp, h( jch+1, jch ), c, s,
580  $ h( jch, jch ) )
581  h( jch+1, jch ) = czero
582  CALL zrot( ilastm-jch, h( jch, jch+1 ), ldh,
583  $ h( jch+1, jch+1 ), ldh, c, s )
584  CALL zrot( ilastm-jch, t( jch, jch+1 ), ldt,
585  $ t( jch+1, jch+1 ), ldt, c, s )
586  IF( ilq )
587  $ CALL zrot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
588  $ c, dconjg( s ) )
589  IF( ilazr2 )
590  $ h( jch, jch-1 ) = h( jch, jch-1 )*c
591  ilazr2 = .false.
592  IF( abs1( t( jch+1, jch+1 ) ).GE.btol ) THEN
593  IF( jch+1.GE.ilast ) THEN
594  GO TO 60
595  ELSE
596  ifirst = jch + 1
597  GO TO 70
598  END IF
599  END IF
600  t( jch+1, jch+1 ) = czero
601  20 CONTINUE
602  GO TO 50
603  ELSE
604 *
605 * Only test 2 passed -- chase the zero to T(ILAST,ILAST)
606 * Then process as in the case T(ILAST,ILAST)=0
607 *
608  DO 30 jch = j, ilast - 1
609  ctemp = t( jch, jch+1 )
610  CALL zlartg( ctemp, t( jch+1, jch+1 ), c, s,
611  $ t( jch, jch+1 ) )
612  t( jch+1, jch+1 ) = czero
613  IF( jch.LT.ilastm-1 )
614  $ CALL zrot( ilastm-jch-1, t( jch, jch+2 ), ldt,
615  $ t( jch+1, jch+2 ), ldt, c, s )
616  CALL zrot( ilastm-jch+2, h( jch, jch-1 ), ldh,
617  $ h( jch+1, jch-1 ), ldh, c, s )
618  IF( ilq )
619  $ CALL zrot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
620  $ c, dconjg( s ) )
621  ctemp = h( jch+1, jch )
622  CALL zlartg( ctemp, h( jch+1, jch-1 ), c, s,
623  $ h( jch+1, jch ) )
624  h( jch+1, jch-1 ) = czero
625  CALL zrot( jch+1-ifrstm, h( ifrstm, jch ), 1,
626  $ h( ifrstm, jch-1 ), 1, c, s )
627  CALL zrot( jch-ifrstm, t( ifrstm, jch ), 1,
628  $ t( ifrstm, jch-1 ), 1, c, s )
629  IF( ilz )
630  $ CALL zrot( n, z( 1, jch ), 1, z( 1, jch-1 ), 1,
631  $ c, s )
632  30 CONTINUE
633  GO TO 50
634  END IF
635  ELSE IF( ilazro ) THEN
636 *
637 * Only test 1 passed -- work on J:ILAST
638 *
639  ifirst = j
640  GO TO 70
641  END IF
642 *
643 * Neither test passed -- try next J
644 *
645  40 CONTINUE
646 *
647 * (Drop-through is "impossible")
648 *
649  info = 2*n + 1
650  GO TO 210
651 *
652 * T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
653 * 1x1 block.
654 *
655  50 CONTINUE
656  ctemp = h( ilast, ilast )
657  CALL zlartg( ctemp, h( ilast, ilast-1 ), c, s,
658  $ h( ilast, ilast ) )
659  h( ilast, ilast-1 ) = czero
660  CALL zrot( ilast-ifrstm, h( ifrstm, ilast ), 1,
661  $ h( ifrstm, ilast-1 ), 1, c, s )
662  CALL zrot( ilast-ifrstm, t( ifrstm, ilast ), 1,
663  $ t( ifrstm, ilast-1 ), 1, c, s )
664  IF( ilz )
665  $ CALL zrot( n, z( 1, ilast ), 1, z( 1, ilast-1 ), 1, c, s )
666 *
667 * H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA
668 *
669  60 CONTINUE
670  absb = abs( t( ilast, ilast ) )
671  IF( absb.GT.safmin ) THEN
672  signbc = dconjg( t( ilast, ilast ) / absb )
673  t( ilast, ilast ) = absb
674  IF( ilschr ) THEN
675  CALL zscal( ilast-ifrstm, signbc, t( ifrstm, ilast ), 1 )
676  CALL zscal( ilast+1-ifrstm, signbc, h( ifrstm, ilast ),
677  $ 1 )
678  ELSE
679  CALL zscal( 1, signbc, h( ilast, ilast ), 1 )
680  END IF
681  IF( ilz )
682  $ CALL zscal( n, signbc, z( 1, ilast ), 1 )
683  ELSE
684  t( ilast, ilast ) = czero
685  END IF
686  alpha( ilast ) = h( ilast, ilast )
687  beta( ilast ) = t( ilast, ilast )
688 *
689 * Go to next block -- exit if finished.
690 *
691  ilast = ilast - 1
692  IF( ilast.LT.ilo )
693  $ GO TO 190
694 *
695 * Reset counters
696 *
697  iiter = 0
698  eshift = czero
699  IF( .NOT.ilschr ) THEN
700  ilastm = ilast
701  IF( ifrstm.GT.ilast )
702  $ ifrstm = ilo
703  END IF
704  GO TO 160
705 *
706 * QZ step
707 *
708 * This iteration only involves rows/columns IFIRST:ILAST. We
709 * assume IFIRST < ILAST, and that the diagonal of B is non-zero.
710 *
711  70 CONTINUE
712  iiter = iiter + 1
713  IF( .NOT.ilschr ) THEN
714  ifrstm = ifirst
715  END IF
716 *
717 * Compute the Shift.
718 *
719 * At this point, IFIRST < ILAST, and the diagonal elements of
720 * T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
721 * magnitude)
722 *
723  IF( ( iiter / 10 )*10.NE.iiter ) THEN
724 *
725 * The Wilkinson shift (AEP p.512), i.e., the eigenvalue of
726 * the bottom-right 2x2 block of A inv(B) which is nearest to
727 * the bottom-right element.
728 *
729 * We factor B as U*D, where U has unit diagonals, and
730 * compute (A*inv(D))*inv(U).
731 *
732  u12 = ( bscale*t( ilast-1, ilast ) ) /
733  $ ( bscale*t( ilast, ilast ) )
734  ad11 = ( ascale*h( ilast-1, ilast-1 ) ) /
735  $ ( bscale*t( ilast-1, ilast-1 ) )
736  ad21 = ( ascale*h( ilast, ilast-1 ) ) /
737  $ ( bscale*t( ilast-1, ilast-1 ) )
738  ad12 = ( ascale*h( ilast-1, ilast ) ) /
739  $ ( bscale*t( ilast, ilast ) )
740  ad22 = ( ascale*h( ilast, ilast ) ) /
741  $ ( bscale*t( ilast, ilast ) )
742  abi22 = ad22 - u12*ad21
743  abi12 = ad12 - u12*ad11
744 *
745  shift = abi22
746  ctemp = sqrt( abi12 )*sqrt( ad21 )
747  temp = abs1( ctemp )
748  IF( ctemp.NE.zero ) THEN
749  x = half*( ad11-shift )
750  temp2 = abs1( x )
751  temp = max( temp, abs1( x ) )
752  y = temp*sqrt( ( x / temp )**2+( ctemp / temp )**2 )
753  IF( temp2.GT.zero ) THEN
754  IF( dble( x / temp2 )*dble( y )+
755  $ dimag( x / temp2 )*dimag( y ).LT.zero )y = -y
756  END IF
757  shift = shift - ctemp*zladiv( ctemp, ( x+y ) )
758  END IF
759  ELSE
760 *
761 * Exceptional shift. Chosen for no particularly good reason.
762 *
763  IF( ( iiter / 20 )*20.EQ.iiter .AND.
764  $ bscale*abs1(t( ilast, ilast )).GT.safmin ) THEN
765  eshift = eshift + ( ascale*h( ilast,
766  $ ilast ) )/( bscale*t( ilast, ilast ) )
767  ELSE
768  eshift = eshift + ( ascale*h( ilast,
769  $ ilast-1 ) )/( bscale*t( ilast-1, ilast-1 ) )
770  END IF
771  shift = eshift
772  END IF
773 *
774 * Now check for two consecutive small subdiagonals.
775 *
776  DO 80 j = ilast - 1, ifirst + 1, -1
777  istart = j
778  ctemp = ascale*h( j, j ) - shift*( bscale*t( j, j ) )
779  temp = abs1( ctemp )
780  temp2 = ascale*abs1( h( j+1, j ) )
781  tempr = max( temp, temp2 )
782  IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
783  temp = temp / tempr
784  temp2 = temp2 / tempr
785  END IF
786  IF( abs1( h( j, j-1 ) )*temp2.LE.temp*atol )
787  $ GO TO 90
788  80 CONTINUE
789 *
790  istart = ifirst
791  ctemp = ascale*h( ifirst, ifirst ) -
792  $ shift*( bscale*t( ifirst, ifirst ) )
793  90 CONTINUE
794 *
795 * Do an implicit-shift QZ sweep.
796 *
797 * Initial Q
798 *
799  ctemp2 = ascale*h( istart+1, istart )
800  CALL zlartg( ctemp, ctemp2, c, s, ctemp3 )
801 *
802 * Sweep
803 *
804  DO 150 j = istart, ilast - 1
805  IF( j.GT.istart ) THEN
806  ctemp = h( j, j-1 )
807  CALL zlartg( ctemp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
808  h( j+1, j-1 ) = czero
809  END IF
810 *
811  DO 100 jc = j, ilastm
812  ctemp = c*h( j, jc ) + s*h( j+1, jc )
813  h( j+1, jc ) = -dconjg( s )*h( j, jc ) + c*h( j+1, jc )
814  h( j, jc ) = ctemp
815  ctemp2 = c*t( j, jc ) + s*t( j+1, jc )
816  t( j+1, jc ) = -dconjg( s )*t( j, jc ) + c*t( j+1, jc )
817  t( j, jc ) = ctemp2
818  100 CONTINUE
819  IF( ilq ) THEN
820  DO 110 jr = 1, n
821  ctemp = c*q( jr, j ) + dconjg( s )*q( jr, j+1 )
822  q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
823  q( jr, j ) = ctemp
824  110 CONTINUE
825  END IF
826 *
827  ctemp = t( j+1, j+1 )
828  CALL zlartg( ctemp, t( j+1, j ), c, s, t( j+1, j+1 ) )
829  t( j+1, j ) = czero
830 *
831  DO 120 jr = ifrstm, min( j+2, ilast )
832  ctemp = c*h( jr, j+1 ) + s*h( jr, j )
833  h( jr, j ) = -dconjg( s )*h( jr, j+1 ) + c*h( jr, j )
834  h( jr, j+1 ) = ctemp
835  120 CONTINUE
836  DO 130 jr = ifrstm, j
837  ctemp = c*t( jr, j+1 ) + s*t( jr, j )
838  t( jr, j ) = -dconjg( s )*t( jr, j+1 ) + c*t( jr, j )
839  t( jr, j+1 ) = ctemp
840  130 CONTINUE
841  IF( ilz ) THEN
842  DO 140 jr = 1, n
843  ctemp = c*z( jr, j+1 ) + s*z( jr, j )
844  z( jr, j ) = -dconjg( s )*z( jr, j+1 ) + c*z( jr, j )
845  z( jr, j+1 ) = ctemp
846  140 CONTINUE
847  END IF
848  150 CONTINUE
849 *
850  160 CONTINUE
851 *
852  170 CONTINUE
853 *
854 * Drop-through = non-convergence
855 *
856  180 CONTINUE
857  info = ilast
858  GO TO 210
859 *
860 * Successful completion of all QZ steps
861 *
862  190 CONTINUE
863 *
864 * Set Eigenvalues 1:ILO-1
865 *
866  DO 200 j = 1, ilo - 1
867  absb = abs( t( j, j ) )
868  IF( absb.GT.safmin ) THEN
869  signbc = dconjg( t( j, j ) / absb )
870  t( j, j ) = absb
871  IF( ilschr ) THEN
872  CALL zscal( j-1, signbc, t( 1, j ), 1 )
873  CALL zscal( j, signbc, h( 1, j ), 1 )
874  ELSE
875  CALL zscal( 1, signbc, h( j, j ), 1 )
876  END IF
877  IF( ilz )
878  $ CALL zscal( n, signbc, z( 1, j ), 1 )
879  ELSE
880  t( j, j ) = czero
881  END IF
882  alpha( j ) = h( j, j )
883  beta( j ) = t( j, j )
884  200 CONTINUE
885 *
886 * Normal Termination
887 *
888  info = 0
889 *
890 * Exit (other than argument error) -- return optimal workspace size
891 *
892  210 CONTINUE
893  work( 1 ) = dcmplx( n )
894  RETURN
895 *
896 * End of ZHGEQZ
897 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition: zlartg.f90:118
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:78
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
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
complex *16 function zladiv(X, Y)
ZLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition: zladiv.f:64
double precision function zlanhs(NORM, N, A, LDA, WORK)
ZLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: zlanhs.f:109
Here is the call graph for this function:
Here is the caller graph for this function: