LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ chgeqz()

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

CHGEQZ

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

Purpose:
 CHGEQZ 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 CGGHRD.

 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 CGGHRD 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 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 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 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 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 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 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 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 REAL 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 chgeqz.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  REAL RWORK( * )
295  COMPLEX ALPHA( * ), BETA( * ), H( LDH, * ),
296  $ Q( LDQ, * ), T( LDT, * ), WORK( * ),
297  $ Z( LDZ, * )
298 * ..
299 *
300 * =====================================================================
301 *
302 * .. Parameters ..
303  COMPLEX CZERO, CONE
304  parameter( czero = ( 0.0e+0, 0.0e+0 ),
305  $ cone = ( 1.0e+0, 0.0e+0 ) )
306  REAL ZERO, ONE
307  parameter( zero = 0.0e+0, one = 1.0e+0 )
308  REAL HALF
309  parameter( half = 0.5e+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  REAL ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,
317  $ C, SAFMIN, TEMP, TEMP2, TEMPR, ULP
318  COMPLEX ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2,
319  $ CTEMP3, ESHIFT, S, SHIFT, SIGNBC,
320  $ U12, X, ABI12, Y
321 * ..
322 * .. External Functions ..
323  COMPLEX CLADIV
324  LOGICAL LSAME
325  REAL CLANHS, SLAMCH
326  EXTERNAL cladiv, lsame, clanhs, slamch
327 * ..
328 * .. External Subroutines ..
329  EXTERNAL clartg, claset, crot, cscal, xerbla
330 * ..
331 * .. Intrinsic Functions ..
332  INTRINSIC abs, aimag, cmplx, conjg, max, min, real, sqrt
333 * ..
334 * .. Statement Functions ..
335  REAL ABS1
336 * ..
337 * .. Statement Function definitions ..
338  abs1( x ) = abs( real( x ) ) + abs( aimag( x ) )
339 * ..
340 * .. Executable Statements ..
341 *
342 * Decode JOB, COMPQ, COMPZ
343 *
344  IF( lsame( job, 'E' ) ) THEN
345  ilschr = .false.
346  ischur = 1
347  ELSE IF( lsame( job, 'S' ) ) THEN
348  ilschr = .true.
349  ischur = 2
350  ELSE
351  ilschr = .true.
352  ischur = 0
353  END IF
354 *
355  IF( lsame( compq, 'N' ) ) THEN
356  ilq = .false.
357  icompq = 1
358  ELSE IF( lsame( compq, 'V' ) ) THEN
359  ilq = .true.
360  icompq = 2
361  ELSE IF( lsame( compq, 'I' ) ) THEN
362  ilq = .true.
363  icompq = 3
364  ELSE
365  ilq = .true.
366  icompq = 0
367  END IF
368 *
369  IF( lsame( compz, 'N' ) ) THEN
370  ilz = .false.
371  icompz = 1
372  ELSE IF( lsame( compz, 'V' ) ) THEN
373  ilz = .true.
374  icompz = 2
375  ELSE IF( lsame( compz, 'I' ) ) THEN
376  ilz = .true.
377  icompz = 3
378  ELSE
379  ilz = .true.
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( 'CHGEQZ', -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 ) = cmplx( 1 )
423  RETURN
424  END IF
425 *
426 * Initialize Q and Z
427 *
428  IF( icompq.EQ.3 )
429  $ CALL claset( 'Full', n, n, czero, cone, q, ldq )
430  IF( icompz.EQ.3 )
431  $ CALL claset( 'Full', n, n, czero, cone, z, ldz )
432 *
433 * Machine Constants
434 *
435  in = ihi + 1 - ilo
436  safmin = slamch( 'S' )
437  ulp = slamch( 'E' )*slamch( 'B' )
438  anorm = clanhs( 'F', in, h( ilo, ilo ), ldh, rwork )
439  bnorm = clanhs( '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 = conjg( t( j, j ) / absb )
452  t( j, j ) = absb
453  IF( ilschr ) THEN
454  CALL cscal( j-1, signbc, t( 1, j ), 1 )
455  CALL cscal( j, signbc, h( 1, j ), 1 )
456  ELSE
457  CALL cscal( 1, signbc, h( j, j ), 1 )
458  END IF
459  IF( ilz )
460  $ CALL cscal( 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.max( safmin, ulp*(
519  $ abs1( h( ilast, ilast ) ) + abs1( h( ilast-1, ilast-1 )
520  $ ) ) ) ) THEN
521  h( ilast, ilast-1 ) = czero
522  GO TO 60
523  END IF
524  END IF
525 *
526  IF( abs( t( ilast, ilast ) ).LE.max( safmin, ulp*(
527  $ abs( t( ilast - 1, ilast ) ) + abs( t( ilast-1, ilast-1 )
528  $ ) ) ) ) THEN
529  t( ilast, ilast ) = czero
530  GO TO 50
531  END IF
532 *
533 * General case: j<ILAST
534 *
535  DO 40 j = ilast - 1, ilo, -1
536 *
537 * Test 1: for H(j,j-1)=0 or j=ILO
538 *
539  IF( j.EQ.ilo ) THEN
540  ilazro = .true.
541  ELSE
542  IF( abs1( h( j, j-1 ) ).LE.max( safmin, ulp*(
543  $ abs1( h( j, j ) ) + abs1( h( j-1, j-1 ) )
544  $ ) ) ) THEN
545  h( j, j-1 ) = czero
546  ilazro = .true.
547  ELSE
548  ilazro = .false.
549  END IF
550  END IF
551 *
552 * Test 2: for T(j,j)=0
553 *
554  temp = abs( t( j, j + 1 ) )
555  IF ( j .GT. ilo )
556  $ temp = temp + abs( t( j - 1, j ) )
557  IF( abs( t( j, j ) ).LT.max( safmin,ulp*temp ) ) THEN
558  t( j, j ) = czero
559 *
560 * Test 1a: Check for 2 consecutive small subdiagonals in A
561 *
562  ilazr2 = .false.
563  IF( .NOT.ilazro ) THEN
564  IF( abs1( h( j, j-1 ) )*( ascale*abs1( h( j+1,
565  $ j ) ) ).LE.abs1( h( j, j ) )*( ascale*atol ) )
566  $ ilazr2 = .true.
567  END IF
568 *
569 * If both tests pass (1 & 2), i.e., the leading diagonal
570 * element of B in the block is zero, split a 1x1 block off
571 * at the top. (I.e., at the J-th row/column) The leading
572 * diagonal element of the remainder can also be zero, so
573 * this may have to be done repeatedly.
574 *
575  IF( ilazro .OR. ilazr2 ) THEN
576  DO 20 jch = j, ilast - 1
577  ctemp = h( jch, jch )
578  CALL clartg( ctemp, h( jch+1, jch ), c, s,
579  $ h( jch, jch ) )
580  h( jch+1, jch ) = czero
581  CALL crot( ilastm-jch, h( jch, jch+1 ), ldh,
582  $ h( jch+1, jch+1 ), ldh, c, s )
583  CALL crot( ilastm-jch, t( jch, jch+1 ), ldt,
584  $ t( jch+1, jch+1 ), ldt, c, s )
585  IF( ilq )
586  $ CALL crot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
587  $ c, conjg( s ) )
588  IF( ilazr2 )
589  $ h( jch, jch-1 ) = h( jch, jch-1 )*c
590  ilazr2 = .false.
591  IF( abs1( t( jch+1, jch+1 ) ).GE.btol ) THEN
592  IF( jch+1.GE.ilast ) THEN
593  GO TO 60
594  ELSE
595  ifirst = jch + 1
596  GO TO 70
597  END IF
598  END IF
599  t( jch+1, jch+1 ) = czero
600  20 CONTINUE
601  GO TO 50
602  ELSE
603 *
604 * Only test 2 passed -- chase the zero to T(ILAST,ILAST)
605 * Then process as in the case T(ILAST,ILAST)=0
606 *
607  DO 30 jch = j, ilast - 1
608  ctemp = t( jch, jch+1 )
609  CALL clartg( ctemp, t( jch+1, jch+1 ), c, s,
610  $ t( jch, jch+1 ) )
611  t( jch+1, jch+1 ) = czero
612  IF( jch.LT.ilastm-1 )
613  $ CALL crot( ilastm-jch-1, t( jch, jch+2 ), ldt,
614  $ t( jch+1, jch+2 ), ldt, c, s )
615  CALL crot( ilastm-jch+2, h( jch, jch-1 ), ldh,
616  $ h( jch+1, jch-1 ), ldh, c, s )
617  IF( ilq )
618  $ CALL crot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
619  $ c, conjg( s ) )
620  ctemp = h( jch+1, jch )
621  CALL clartg( ctemp, h( jch+1, jch-1 ), c, s,
622  $ h( jch+1, jch ) )
623  h( jch+1, jch-1 ) = czero
624  CALL crot( jch+1-ifrstm, h( ifrstm, jch ), 1,
625  $ h( ifrstm, jch-1 ), 1, c, s )
626  CALL crot( jch-ifrstm, t( ifrstm, jch ), 1,
627  $ t( ifrstm, jch-1 ), 1, c, s )
628  IF( ilz )
629  $ CALL crot( n, z( 1, jch ), 1, z( 1, jch-1 ), 1,
630  $ c, s )
631  30 CONTINUE
632  GO TO 50
633  END IF
634  ELSE IF( ilazro ) THEN
635 *
636 * Only test 1 passed -- work on J:ILAST
637 *
638  ifirst = j
639  GO TO 70
640  END IF
641 *
642 * Neither test passed -- try next J
643 *
644  40 CONTINUE
645 *
646 * (Drop-through is "impossible")
647 *
648  info = 2*n + 1
649  GO TO 210
650 *
651 * T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
652 * 1x1 block.
653 *
654  50 CONTINUE
655  ctemp = h( ilast, ilast )
656  CALL clartg( ctemp, h( ilast, ilast-1 ), c, s,
657  $ h( ilast, ilast ) )
658  h( ilast, ilast-1 ) = czero
659  CALL crot( ilast-ifrstm, h( ifrstm, ilast ), 1,
660  $ h( ifrstm, ilast-1 ), 1, c, s )
661  CALL crot( ilast-ifrstm, t( ifrstm, ilast ), 1,
662  $ t( ifrstm, ilast-1 ), 1, c, s )
663  IF( ilz )
664  $ CALL crot( n, z( 1, ilast ), 1, z( 1, ilast-1 ), 1, c, s )
665 *
666 * H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA
667 *
668  60 CONTINUE
669  absb = abs( t( ilast, ilast ) )
670  IF( absb.GT.safmin ) THEN
671  signbc = conjg( t( ilast, ilast ) / absb )
672  t( ilast, ilast ) = absb
673  IF( ilschr ) THEN
674  CALL cscal( ilast-ifrstm, signbc, t( ifrstm, ilast ), 1 )
675  CALL cscal( ilast+1-ifrstm, signbc, h( ifrstm, ilast ),
676  $ 1 )
677  ELSE
678  CALL cscal( 1, signbc, h( ilast, ilast ), 1 )
679  END IF
680  IF( ilz )
681  $ CALL cscal( n, signbc, z( 1, ilast ), 1 )
682  ELSE
683  t( ilast, ilast ) = czero
684  END IF
685  alpha( ilast ) = h( ilast, ilast )
686  beta( ilast ) = t( ilast, ilast )
687 *
688 * Go to next block -- exit if finished.
689 *
690  ilast = ilast - 1
691  IF( ilast.LT.ilo )
692  $ GO TO 190
693 *
694 * Reset counters
695 *
696  iiter = 0
697  eshift = czero
698  IF( .NOT.ilschr ) THEN
699  ilastm = ilast
700  IF( ifrstm.GT.ilast )
701  $ ifrstm = ilo
702  END IF
703  GO TO 160
704 *
705 * QZ step
706 *
707 * This iteration only involves rows/columns IFIRST:ILAST. We
708 * assume IFIRST < ILAST, and that the diagonal of B is non-zero.
709 *
710  70 CONTINUE
711  iiter = iiter + 1
712  IF( .NOT.ilschr ) THEN
713  ifrstm = ifirst
714  END IF
715 *
716 * Compute the Shift.
717 *
718 * At this point, IFIRST < ILAST, and the diagonal elements of
719 * T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
720 * magnitude)
721 *
722  IF( ( iiter / 10 )*10.NE.iiter ) THEN
723 *
724 * The Wilkinson shift (AEP p.512), i.e., the eigenvalue of
725 * the bottom-right 2x2 block of A inv(B) which is nearest to
726 * the bottom-right element.
727 *
728 * We factor B as U*D, where U has unit diagonals, and
729 * compute (A*inv(D))*inv(U).
730 *
731  u12 = ( bscale*t( ilast-1, ilast ) ) /
732  $ ( bscale*t( ilast, ilast ) )
733  ad11 = ( ascale*h( ilast-1, ilast-1 ) ) /
734  $ ( bscale*t( ilast-1, ilast-1 ) )
735  ad21 = ( ascale*h( ilast, ilast-1 ) ) /
736  $ ( bscale*t( ilast-1, ilast-1 ) )
737  ad12 = ( ascale*h( ilast-1, ilast ) ) /
738  $ ( bscale*t( ilast, ilast ) )
739  ad22 = ( ascale*h( ilast, ilast ) ) /
740  $ ( bscale*t( ilast, ilast ) )
741  abi22 = ad22 - u12*ad21
742  abi12 = ad12 - u12*ad11
743 *
744  shift = abi22
745  ctemp = sqrt( abi12 )*sqrt( ad21 )
746  temp = abs1( ctemp )
747  IF( ctemp.NE.zero ) THEN
748  x = half*( ad11-shift )
749  temp2 = abs1( x )
750  temp = max( temp, abs1( x ) )
751  y = temp*sqrt( ( x / temp )**2+( ctemp / temp )**2 )
752  IF( temp2.GT.zero ) THEN
753  IF( real( x / temp2 )*real( y )+
754  $ aimag( x / temp2 )*aimag( y ).LT.zero )y = -y
755  END IF
756  shift = shift - ctemp*cladiv( ctemp, ( x+y ) )
757  END IF
758  ELSE
759 *
760 * Exceptional shift. Chosen for no particularly good reason.
761 *
762  IF( ( iiter / 20 )*20.EQ.iiter .AND.
763  $ bscale*abs1(t( ilast, ilast )).GT.safmin ) THEN
764  eshift = eshift + ( ascale*h( ilast,
765  $ ilast ) )/( bscale*t( ilast, ilast ) )
766  ELSE
767  eshift = eshift + ( ascale*h( ilast,
768  $ ilast-1 ) )/( bscale*t( ilast-1, ilast-1 ) )
769  END IF
770  shift = eshift
771  END IF
772 *
773 * Now check for two consecutive small subdiagonals.
774 *
775  DO 80 j = ilast - 1, ifirst + 1, -1
776  istart = j
777  ctemp = ascale*h( j, j ) - shift*( bscale*t( j, j ) )
778  temp = abs1( ctemp )
779  temp2 = ascale*abs1( h( j+1, j ) )
780  tempr = max( temp, temp2 )
781  IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
782  temp = temp / tempr
783  temp2 = temp2 / tempr
784  END IF
785  IF( abs1( h( j, j-1 ) )*temp2.LE.temp*atol )
786  $ GO TO 90
787  80 CONTINUE
788 *
789  istart = ifirst
790  ctemp = ascale*h( ifirst, ifirst ) -
791  $ shift*( bscale*t( ifirst, ifirst ) )
792  90 CONTINUE
793 *
794 * Do an implicit-shift QZ sweep.
795 *
796 * Initial Q
797 *
798  ctemp2 = ascale*h( istart+1, istart )
799  CALL clartg( ctemp, ctemp2, c, s, ctemp3 )
800 *
801 * Sweep
802 *
803  DO 150 j = istart, ilast - 1
804  IF( j.GT.istart ) THEN
805  ctemp = h( j, j-1 )
806  CALL clartg( ctemp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
807  h( j+1, j-1 ) = czero
808  END IF
809 *
810  DO 100 jc = j, ilastm
811  ctemp = c*h( j, jc ) + s*h( j+1, jc )
812  h( j+1, jc ) = -conjg( s )*h( j, jc ) + c*h( j+1, jc )
813  h( j, jc ) = ctemp
814  ctemp2 = c*t( j, jc ) + s*t( j+1, jc )
815  t( j+1, jc ) = -conjg( s )*t( j, jc ) + c*t( j+1, jc )
816  t( j, jc ) = ctemp2
817  100 CONTINUE
818  IF( ilq ) THEN
819  DO 110 jr = 1, n
820  ctemp = c*q( jr, j ) + conjg( s )*q( jr, j+1 )
821  q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
822  q( jr, j ) = ctemp
823  110 CONTINUE
824  END IF
825 *
826  ctemp = t( j+1, j+1 )
827  CALL clartg( ctemp, t( j+1, j ), c, s, t( j+1, j+1 ) )
828  t( j+1, j ) = czero
829 *
830  DO 120 jr = ifrstm, min( j+2, ilast )
831  ctemp = c*h( jr, j+1 ) + s*h( jr, j )
832  h( jr, j ) = -conjg( s )*h( jr, j+1 ) + c*h( jr, j )
833  h( jr, j+1 ) = ctemp
834  120 CONTINUE
835  DO 130 jr = ifrstm, j
836  ctemp = c*t( jr, j+1 ) + s*t( jr, j )
837  t( jr, j ) = -conjg( s )*t( jr, j+1 ) + c*t( jr, j )
838  t( jr, j+1 ) = ctemp
839  130 CONTINUE
840  IF( ilz ) THEN
841  DO 140 jr = 1, n
842  ctemp = c*z( jr, j+1 ) + s*z( jr, j )
843  z( jr, j ) = -conjg( s )*z( jr, j+1 ) + c*z( jr, j )
844  z( jr, j+1 ) = ctemp
845  140 CONTINUE
846  END IF
847  150 CONTINUE
848 *
849  160 CONTINUE
850 *
851  170 CONTINUE
852 *
853 * Drop-through = non-convergence
854 *
855  180 CONTINUE
856  info = ilast
857  GO TO 210
858 *
859 * Successful completion of all QZ steps
860 *
861  190 CONTINUE
862 *
863 * Set Eigenvalues 1:ILO-1
864 *
865  DO 200 j = 1, ilo - 1
866  absb = abs( t( j, j ) )
867  IF( absb.GT.safmin ) THEN
868  signbc = conjg( t( j, j ) / absb )
869  t( j, j ) = absb
870  IF( ilschr ) THEN
871  CALL cscal( j-1, signbc, t( 1, j ), 1 )
872  CALL cscal( j, signbc, h( 1, j ), 1 )
873  ELSE
874  CALL cscal( 1, signbc, h( j, j ), 1 )
875  END IF
876  IF( ilz )
877  $ CALL cscal( n, signbc, z( 1, j ), 1 )
878  ELSE
879  t( j, j ) = czero
880  END IF
881  alpha( j ) = h( j, j )
882  beta( j ) = t( j, j )
883  200 CONTINUE
884 *
885 * Normal Termination
886 *
887  info = 0
888 *
889 * Exit (other than argument error) -- return optimal workspace size
890 *
891  210 CONTINUE
892  work( 1 ) = cmplx( n )
893  RETURN
894 *
895 * End of CHGEQZ
896 *
subroutine clartg(f, g, c, s, r)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition: clartg.f90:118
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
real function clanhs(NORM, N, A, LDA, WORK)
CLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: clanhs.f:109
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
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
complex function cladiv(X, Y)
CLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition: cladiv.f:64
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: