LAPACK  3.8.0
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.
Date
April 2012
Further Details:
  We assume that complex ABS works as long as its value is less than
  overflow.

Definition at line 286 of file zhgeqz.f.

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