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

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