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

◆ 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.btol ) THEN
527 t( ilast, ilast ) = czero
528 GO TO 50
529 END IF
530*
531* General case: j<ILAST
532*
533 DO 40 j = ilast - 1, ilo, -1
534*
535* Test 1: for H(j,j-1)=0 or j=ILO
536*
537 IF( j.EQ.ilo ) THEN
538 ilazro = .true.
539 ELSE
540 IF( abs1( h( j, j-1 ) ).LE.max( safmin, ulp*(
541 $ abs1( h( j, j ) ) + abs1( h( j-1, j-1 ) )
542 $ ) ) ) THEN
543 h( j, j-1 ) = czero
544 ilazro = .true.
545 ELSE
546 ilazro = .false.
547 END IF
548 END IF
549*
550* Test 2: for T(j,j)=0
551*
552 IF( abs( t( j, j ) ).LT.btol ) THEN
553 t( j, j ) = czero
554*
555* Test 1a: Check for 2 consecutive small subdiagonals in A
556*
557 ilazr2 = .false.
558 IF( .NOT.ilazro ) THEN
559 IF( abs1( h( j, j-1 ) )*( ascale*abs1( h( j+1,
560 $ j ) ) ).LE.abs1( h( j, j ) )*( ascale*atol ) )
561 $ ilazr2 = .true.
562 END IF
563*
564* If both tests pass (1 & 2), i.e., the leading diagonal
565* element of B in the block is zero, split a 1x1 block off
566* at the top. (I.e., at the J-th row/column) The leading
567* diagonal element of the remainder can also be zero, so
568* this may have to be done repeatedly.
569*
570 IF( ilazro .OR. ilazr2 ) THEN
571 DO 20 jch = j, ilast - 1
572 ctemp = h( jch, jch )
573 CALL clartg( ctemp, h( jch+1, jch ), c, s,
574 $ h( jch, jch ) )
575 h( jch+1, jch ) = czero
576 CALL crot( ilastm-jch, h( jch, jch+1 ), ldh,
577 $ h( jch+1, jch+1 ), ldh, c, s )
578 CALL crot( ilastm-jch, t( jch, jch+1 ), ldt,
579 $ t( jch+1, jch+1 ), ldt, c, s )
580 IF( ilq )
581 $ CALL crot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
582 $ c, conjg( s ) )
583 IF( ilazr2 )
584 $ h( jch, jch-1 ) = h( jch, jch-1 )*c
585 ilazr2 = .false.
586 IF( abs1( t( jch+1, jch+1 ) ).GE.btol ) THEN
587 IF( jch+1.GE.ilast ) THEN
588 GO TO 60
589 ELSE
590 ifirst = jch + 1
591 GO TO 70
592 END IF
593 END IF
594 t( jch+1, jch+1 ) = czero
595 20 CONTINUE
596 GO TO 50
597 ELSE
598*
599* Only test 2 passed -- chase the zero to T(ILAST,ILAST)
600* Then process as in the case T(ILAST,ILAST)=0
601*
602 DO 30 jch = j, ilast - 1
603 ctemp = t( jch, jch+1 )
604 CALL clartg( ctemp, t( jch+1, jch+1 ), c, s,
605 $ t( jch, jch+1 ) )
606 t( jch+1, jch+1 ) = czero
607 IF( jch.LT.ilastm-1 )
608 $ CALL crot( ilastm-jch-1, t( jch, jch+2 ), ldt,
609 $ t( jch+1, jch+2 ), ldt, c, s )
610 CALL crot( ilastm-jch+2, h( jch, jch-1 ), ldh,
611 $ h( jch+1, jch-1 ), ldh, c, s )
612 IF( ilq )
613 $ CALL crot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
614 $ c, conjg( s ) )
615 ctemp = h( jch+1, jch )
616 CALL clartg( ctemp, h( jch+1, jch-1 ), c, s,
617 $ h( jch+1, jch ) )
618 h( jch+1, jch-1 ) = czero
619 CALL crot( jch+1-ifrstm, h( ifrstm, jch ), 1,
620 $ h( ifrstm, jch-1 ), 1, c, s )
621 CALL crot( jch-ifrstm, t( ifrstm, jch ), 1,
622 $ t( ifrstm, jch-1 ), 1, c, s )
623 IF( ilz )
624 $ CALL crot( n, z( 1, jch ), 1, z( 1, jch-1 ), 1,
625 $ c, s )
626 30 CONTINUE
627 GO TO 50
628 END IF
629 ELSE IF( ilazro ) THEN
630*
631* Only test 1 passed -- work on J:ILAST
632*
633 ifirst = j
634 GO TO 70
635 END IF
636*
637* Neither test passed -- try next J
638*
639 40 CONTINUE
640*
641* (Drop-through is "impossible")
642*
643 info = 2*n + 1
644 GO TO 210
645*
646* T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
647* 1x1 block.
648*
649 50 CONTINUE
650 ctemp = h( ilast, ilast )
651 CALL clartg( ctemp, h( ilast, ilast-1 ), c, s,
652 $ h( ilast, ilast ) )
653 h( ilast, ilast-1 ) = czero
654 CALL crot( ilast-ifrstm, h( ifrstm, ilast ), 1,
655 $ h( ifrstm, ilast-1 ), 1, c, s )
656 CALL crot( ilast-ifrstm, t( ifrstm, ilast ), 1,
657 $ t( ifrstm, ilast-1 ), 1, c, s )
658 IF( ilz )
659 $ CALL crot( n, z( 1, ilast ), 1, z( 1, ilast-1 ), 1, c, s )
660*
661* H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA
662*
663 60 CONTINUE
664 absb = abs( t( ilast, ilast ) )
665 IF( absb.GT.safmin ) THEN
666 signbc = conjg( t( ilast, ilast ) / absb )
667 t( ilast, ilast ) = absb
668 IF( ilschr ) THEN
669 CALL cscal( ilast-ifrstm, signbc, t( ifrstm, ilast ), 1 )
670 CALL cscal( ilast+1-ifrstm, signbc, h( ifrstm, ilast ),
671 $ 1 )
672 ELSE
673 CALL cscal( 1, signbc, h( ilast, ilast ), 1 )
674 END IF
675 IF( ilz )
676 $ CALL cscal( n, signbc, z( 1, ilast ), 1 )
677 ELSE
678 t( ilast, ilast ) = czero
679 END IF
680 alpha( ilast ) = h( ilast, ilast )
681 beta( ilast ) = t( ilast, ilast )
682*
683* Go to next block -- exit if finished.
684*
685 ilast = ilast - 1
686 IF( ilast.LT.ilo )
687 $ GO TO 190
688*
689* Reset counters
690*
691 iiter = 0
692 eshift = czero
693 IF( .NOT.ilschr ) THEN
694 ilastm = ilast
695 IF( ifrstm.GT.ilast )
696 $ ifrstm = ilo
697 END IF
698 GO TO 160
699*
700* QZ step
701*
702* This iteration only involves rows/columns IFIRST:ILAST. We
703* assume IFIRST < ILAST, and that the diagonal of B is non-zero.
704*
705 70 CONTINUE
706 iiter = iiter + 1
707 IF( .NOT.ilschr ) THEN
708 ifrstm = ifirst
709 END IF
710*
711* Compute the Shift.
712*
713* At this point, IFIRST < ILAST, and the diagonal elements of
714* T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
715* magnitude)
716*
717 IF( ( iiter / 10 )*10.NE.iiter ) THEN
718*
719* The Wilkinson shift (AEP p.512), i.e., the eigenvalue of
720* the bottom-right 2x2 block of A inv(B) which is nearest to
721* the bottom-right element.
722*
723* We factor B as U*D, where U has unit diagonals, and
724* compute (A*inv(D))*inv(U).
725*
726 u12 = ( bscale*t( ilast-1, ilast ) ) /
727 $ ( bscale*t( ilast, ilast ) )
728 ad11 = ( ascale*h( ilast-1, ilast-1 ) ) /
729 $ ( bscale*t( ilast-1, ilast-1 ) )
730 ad21 = ( ascale*h( ilast, ilast-1 ) ) /
731 $ ( bscale*t( ilast-1, ilast-1 ) )
732 ad12 = ( ascale*h( ilast-1, ilast ) ) /
733 $ ( bscale*t( ilast, ilast ) )
734 ad22 = ( ascale*h( ilast, ilast ) ) /
735 $ ( bscale*t( ilast, ilast ) )
736 abi22 = ad22 - u12*ad21
737 abi12 = ad12 - u12*ad11
738*
739 shift = abi22
740 ctemp = sqrt( abi12 )*sqrt( ad21 )
741 temp = abs1( ctemp )
742 IF( ctemp.NE.zero ) THEN
743 x = half*( ad11-shift )
744 temp2 = abs1( x )
745 temp = max( temp, abs1( x ) )
746 y = temp*sqrt( ( x / temp )**2+( ctemp / temp )**2 )
747 IF( temp2.GT.zero ) THEN
748 IF( real( x / temp2 )*real( y )+
749 $ aimag( x / temp2 )*aimag( y ).LT.zero )y = -y
750 END IF
751 shift = shift - ctemp*cladiv( ctemp, ( x+y ) )
752 END IF
753 ELSE
754*
755* Exceptional shift. Chosen for no particularly good reason.
756*
757 IF( ( iiter / 20 )*20.EQ.iiter .AND.
758 $ bscale*abs1(t( ilast, ilast )).GT.safmin ) THEN
759 eshift = eshift + ( ascale*h( ilast,
760 $ ilast ) )/( bscale*t( ilast, ilast ) )
761 ELSE
762 eshift = eshift + ( ascale*h( ilast,
763 $ ilast-1 ) )/( bscale*t( ilast-1, ilast-1 ) )
764 END IF
765 shift = eshift
766 END IF
767*
768* Now check for two consecutive small subdiagonals.
769*
770 DO 80 j = ilast - 1, ifirst + 1, -1
771 istart = j
772 ctemp = ascale*h( j, j ) - shift*( bscale*t( j, j ) )
773 temp = abs1( ctemp )
774 temp2 = ascale*abs1( h( j+1, j ) )
775 tempr = max( temp, temp2 )
776 IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
777 temp = temp / tempr
778 temp2 = temp2 / tempr
779 END IF
780 IF( abs1( h( j, j-1 ) )*temp2.LE.temp*atol )
781 $ GO TO 90
782 80 CONTINUE
783*
784 istart = ifirst
785 ctemp = ascale*h( ifirst, ifirst ) -
786 $ shift*( bscale*t( ifirst, ifirst ) )
787 90 CONTINUE
788*
789* Do an implicit-shift QZ sweep.
790*
791* Initial Q
792*
793 ctemp2 = ascale*h( istart+1, istart )
794 CALL clartg( ctemp, ctemp2, c, s, ctemp3 )
795*
796* Sweep
797*
798 DO 150 j = istart, ilast - 1
799 IF( j.GT.istart ) THEN
800 ctemp = h( j, j-1 )
801 CALL clartg( ctemp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
802 h( j+1, j-1 ) = czero
803 END IF
804*
805 DO 100 jc = j, ilastm
806 ctemp = c*h( j, jc ) + s*h( j+1, jc )
807 h( j+1, jc ) = -conjg( s )*h( j, jc ) + c*h( j+1, jc )
808 h( j, jc ) = ctemp
809 ctemp2 = c*t( j, jc ) + s*t( j+1, jc )
810 t( j+1, jc ) = -conjg( s )*t( j, jc ) + c*t( j+1, jc )
811 t( j, jc ) = ctemp2
812 100 CONTINUE
813 IF( ilq ) THEN
814 DO 110 jr = 1, n
815 ctemp = c*q( jr, j ) + conjg( s )*q( jr, j+1 )
816 q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
817 q( jr, j ) = ctemp
818 110 CONTINUE
819 END IF
820*
821 ctemp = t( j+1, j+1 )
822 CALL clartg( ctemp, t( j+1, j ), c, s, t( j+1, j+1 ) )
823 t( j+1, j ) = czero
824*
825 DO 120 jr = ifrstm, min( j+2, ilast )
826 ctemp = c*h( jr, j+1 ) + s*h( jr, j )
827 h( jr, j ) = -conjg( s )*h( jr, j+1 ) + c*h( jr, j )
828 h( jr, j+1 ) = ctemp
829 120 CONTINUE
830 DO 130 jr = ifrstm, j
831 ctemp = c*t( jr, j+1 ) + s*t( jr, j )
832 t( jr, j ) = -conjg( s )*t( jr, j+1 ) + c*t( jr, j )
833 t( jr, j+1 ) = ctemp
834 130 CONTINUE
835 IF( ilz ) THEN
836 DO 140 jr = 1, n
837 ctemp = c*z( jr, j+1 ) + s*z( jr, j )
838 z( jr, j ) = -conjg( s )*z( jr, j+1 ) + c*z( jr, j )
839 z( jr, j+1 ) = ctemp
840 140 CONTINUE
841 END IF
842 150 CONTINUE
843*
844 160 CONTINUE
845*
846 170 CONTINUE
847*
848* Drop-through = non-convergence
849*
850 180 CONTINUE
851 info = ilast
852 GO TO 210
853*
854* Successful completion of all QZ steps
855*
856 190 CONTINUE
857*
858* Set Eigenvalues 1:ILO-1
859*
860 DO 200 j = 1, ilo - 1
861 absb = abs( t( j, j ) )
862 IF( absb.GT.safmin ) THEN
863 signbc = conjg( t( j, j ) / absb )
864 t( j, j ) = absb
865 IF( ilschr ) THEN
866 CALL cscal( j-1, signbc, t( 1, j ), 1 )
867 CALL cscal( j, signbc, h( 1, j ), 1 )
868 ELSE
869 CALL cscal( 1, signbc, h( j, j ), 1 )
870 END IF
871 IF( ilz )
872 $ CALL cscal( n, signbc, z( 1, j ), 1 )
873 ELSE
874 t( j, j ) = czero
875 END IF
876 alpha( j ) = h( j, j )
877 beta( j ) = t( j, j )
878 200 CONTINUE
879*
880* Normal Termination
881*
882 info = 0
883*
884* Exit (other than argument error) -- return optimal workspace size
885*
886 210 CONTINUE
887 work( 1 ) = cmplx( n )
888 RETURN
889*
890* End of CHGEQZ
891*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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
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 clartg(f, g, c, s, r)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition clartg.f90:116
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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
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
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
Here is the call graph for this function:
Here is the caller graph for this function: