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

◆ dhgeqz()

subroutine dhgeqz ( character  job,
character  compq,
character  compz,
integer  n,
integer  ilo,
integer  ihi,
double precision, dimension( ldh, * )  h,
integer  ldh,
double precision, dimension( ldt, * )  t,
integer  ldt,
double precision, dimension( * )  alphar,
double precision, dimension( * )  alphai,
double precision, dimension( * )  beta,
double precision, dimension( ldq, * )  q,
integer  ldq,
double precision, dimension( ldz, * )  z,
integer  ldz,
double precision, dimension( * )  work,
integer  lwork,
integer  info 
)

DHGEQZ

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

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

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

 as computed by DGGHRD.

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

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

 where Q and Z are orthogonal matrices, P is an upper triangular
 matrix, and S is a quasi-triangular matrix with 1-by-1 and 2-by-2
 diagonal blocks.

 The 1-by-1 blocks correspond to real eigenvalues of the matrix pair
 (H,T) and the 2-by-2 blocks correspond to complex conjugate pairs of
 eigenvalues.

 Additionally, the 2-by-2 upper triangular diagonal blocks of P
 corresponding to 2-by-2 blocks of S are reduced to positive diagonal
 form, i.e., if S(j+1,j) is non-zero, then P(j+1,j) = P(j,j+1) = 0,
 P(j,j) > 0, and P(j+1,j+1) > 0.

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

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

 To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently,
 of (A,B)) are computed as a pair of values (alpha,beta), where alpha is
 complex and beta real.
 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.
 Real eigenvalues 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': Compute 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 an orthogonal 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': Z is initialized to the unit matrix and the matrix Z
                 of right Schur vectors of (H,T) is returned;
          = 'V': Z must contain an orthogonal 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 DOUBLE PRECISION array, dimension (LDH, N)
          On entry, the N-by-N upper Hessenberg matrix H.
          On exit, if JOB = 'S', H contains the upper quasi-triangular
          matrix S from the generalized Schur factorization.
          If JOB = 'E', the diagonal blocks of H match those 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 DOUBLE PRECISION 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;
          2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks of S
          are reduced to positive diagonal form, i.e., if H(j+1,j) is
          non-zero, then T(j+1,j) = T(j,j+1) = 0, T(j,j) > 0, and
          T(j+1,j+1) > 0.
          If JOB = 'E', the diagonal blocks of T match those 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]ALPHAR
          ALPHAR is DOUBLE PRECISION array, dimension (N)
          The real parts of each scalar alpha defining an eigenvalue
          of GNEP.
[out]ALPHAI
          ALPHAI is DOUBLE PRECISION array, dimension (N)
          The imaginary parts of each scalar alpha defining an
          eigenvalue of GNEP.
          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
          positive, then the j-th and (j+1)-st eigenvalues are a
          complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j).
[out]BETA
          BETA is DOUBLE PRECISION array, dimension (N)
          The scalars beta that define the eigenvalues of GNEP.
          Together, the quantities alpha = (ALPHAR(j),ALPHAI(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 DOUBLE PRECISION array, dimension (LDQ, N)
          On entry, if COMPQ = 'V', the orthogonal matrix Q1 used in
          the reduction of (A,B) to generalized Hessenberg form.
          On exit, if COMPQ = 'I', the orthogonal matrix of left Schur
          vectors of (H,T), and if COMPQ = 'V', the orthogonal 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 DOUBLE PRECISION array, dimension (LDZ, N)
          On entry, if COMPZ = 'V', the orthogonal matrix Z1 used in
          the reduction of (A,B) to generalized Hessenberg form.
          On exit, if COMPZ = 'I', the orthogonal matrix of
          right Schur vectors of (H,T), and if COMPZ = 'V', the
          orthogonal 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 DOUBLE PRECISION 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]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 ALPHAR(i), ALPHAI(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 ALPHAR(i), ALPHAI(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:
  Iteration counters:

  JITER  -- counts iterations.
  IITER  -- counts iterations run since ILAST was last
            changed.  This is therefore reset only when a 1-by-1 or
            2-by-2 block deflates off the bottom.

Definition at line 301 of file dhgeqz.f.

304*
305* -- LAPACK computational routine --
306* -- LAPACK is a software package provided by Univ. of Tennessee, --
307* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
308*
309* .. Scalar Arguments ..
310 CHARACTER COMPQ, COMPZ, JOB
311 INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
312* ..
313* .. Array Arguments ..
314 DOUBLE PRECISION ALPHAI( * ), ALPHAR( * ), BETA( * ),
315 $ H( LDH, * ), Q( LDQ, * ), T( LDT, * ),
316 $ WORK( * ), Z( LDZ, * )
317* ..
318*
319* =====================================================================
320*
321* .. Parameters ..
322* $ SAFETY = 1.0E+0 )
323 DOUBLE PRECISION HALF, ZERO, ONE, SAFETY
324 parameter( half = 0.5d+0, zero = 0.0d+0, one = 1.0d+0,
325 $ safety = 1.0d+2 )
326* ..
327* .. Local Scalars ..
328 LOGICAL ILAZR2, ILAZRO, ILPIVT, ILQ, ILSCHR, ILZ,
329 $ LQUERY
330 INTEGER ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
331 $ ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER,
332 $ JR, MAXIT
333 DOUBLE PRECISION A11, A12, A1I, A1R, A21, A22, A2I, A2R, AD11,
334 $ AD11L, AD12, AD12L, AD21, AD21L, AD22, AD22L,
335 $ AD32L, AN, ANORM, ASCALE, ATOL, B11, B1A, B1I,
336 $ B1R, B22, B2A, B2I, B2R, BN, BNORM, BSCALE,
337 $ BTOL, C, C11I, C11R, C12, C21, C22I, C22R, CL,
338 $ CQ, CR, CZ, ESHIFT, S, S1, S1INV, S2, SAFMAX,
339 $ SAFMIN, SCALE, SL, SQI, SQR, SR, SZI, SZR, T1,
340 $ T2, T3, TAU, TEMP, TEMP2, TEMPI, TEMPR, U1,
341 $ U12, U12L, U2, ULP, VS, W11, W12, W21, W22,
342 $ WABS, WI, WR, WR2
343* ..
344* .. Local Arrays ..
345 DOUBLE PRECISION V( 3 )
346* ..
347* .. External Functions ..
348 LOGICAL LSAME
349 DOUBLE PRECISION DLAMCH, DLANHS, DLAPY2, DLAPY3
350 EXTERNAL lsame, dlamch, dlanhs, dlapy2, dlapy3
351* ..
352* .. External Subroutines ..
353 EXTERNAL dlag2, dlarfg, dlartg, dlaset, dlasv2, drot,
354 $ xerbla
355* ..
356* .. Intrinsic Functions ..
357 INTRINSIC abs, dble, max, min, sqrt
358* ..
359* .. Executable Statements ..
360*
361* Decode JOB, COMPQ, COMPZ
362*
363 IF( lsame( job, 'E' ) ) THEN
364 ilschr = .false.
365 ischur = 1
366 ELSE IF( lsame( job, 'S' ) ) THEN
367 ilschr = .true.
368 ischur = 2
369 ELSE
370 ischur = 0
371 END IF
372*
373 IF( lsame( compq, 'N' ) ) THEN
374 ilq = .false.
375 icompq = 1
376 ELSE IF( lsame( compq, 'V' ) ) THEN
377 ilq = .true.
378 icompq = 2
379 ELSE IF( lsame( compq, 'I' ) ) THEN
380 ilq = .true.
381 icompq = 3
382 ELSE
383 icompq = 0
384 END IF
385*
386 IF( lsame( compz, 'N' ) ) THEN
387 ilz = .false.
388 icompz = 1
389 ELSE IF( lsame( compz, 'V' ) ) THEN
390 ilz = .true.
391 icompz = 2
392 ELSE IF( lsame( compz, 'I' ) ) THEN
393 ilz = .true.
394 icompz = 3
395 ELSE
396 icompz = 0
397 END IF
398*
399* Check Argument Values
400*
401 info = 0
402 work( 1 ) = max( 1, n )
403 lquery = ( lwork.EQ.-1 )
404 IF( ischur.EQ.0 ) THEN
405 info = -1
406 ELSE IF( icompq.EQ.0 ) THEN
407 info = -2
408 ELSE IF( icompz.EQ.0 ) THEN
409 info = -3
410 ELSE IF( n.LT.0 ) THEN
411 info = -4
412 ELSE IF( ilo.LT.1 ) THEN
413 info = -5
414 ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 ) THEN
415 info = -6
416 ELSE IF( ldh.LT.n ) THEN
417 info = -8
418 ELSE IF( ldt.LT.n ) THEN
419 info = -10
420 ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) ) THEN
421 info = -15
422 ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) ) THEN
423 info = -17
424 ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
425 info = -19
426 END IF
427 IF( info.NE.0 ) THEN
428 CALL xerbla( 'DHGEQZ', -info )
429 RETURN
430 ELSE IF( lquery ) THEN
431 RETURN
432 END IF
433*
434* Quick return if possible
435*
436 IF( n.LE.0 ) THEN
437 work( 1 ) = dble( 1 )
438 RETURN
439 END IF
440*
441* Initialize Q and Z
442*
443 IF( icompq.EQ.3 )
444 $ CALL dlaset( 'Full', n, n, zero, one, q, ldq )
445 IF( icompz.EQ.3 )
446 $ CALL dlaset( 'Full', n, n, zero, one, z, ldz )
447*
448* Machine Constants
449*
450 in = ihi + 1 - ilo
451 safmin = dlamch( 'S' )
452 safmax = one / safmin
453 ulp = dlamch( 'E' )*dlamch( 'B' )
454 anorm = dlanhs( 'F', in, h( ilo, ilo ), ldh, work )
455 bnorm = dlanhs( 'F', in, t( ilo, ilo ), ldt, work )
456 atol = max( safmin, ulp*anorm )
457 btol = max( safmin, ulp*bnorm )
458 ascale = one / max( safmin, anorm )
459 bscale = one / max( safmin, bnorm )
460*
461* Set Eigenvalues IHI+1:N
462*
463 DO 30 j = ihi + 1, n
464 IF( t( j, j ).LT.zero ) THEN
465 IF( ilschr ) THEN
466 DO 10 jr = 1, j
467 h( jr, j ) = -h( jr, j )
468 t( jr, j ) = -t( jr, j )
469 10 CONTINUE
470 ELSE
471 h( j, j ) = -h( j, j )
472 t( j, j ) = -t( j, j )
473 END IF
474 IF( ilz ) THEN
475 DO 20 jr = 1, n
476 z( jr, j ) = -z( jr, j )
477 20 CONTINUE
478 END IF
479 END IF
480 alphar( j ) = h( j, j )
481 alphai( j ) = zero
482 beta( j ) = t( j, j )
483 30 CONTINUE
484*
485* If IHI < ILO, skip QZ steps
486*
487 IF( ihi.LT.ilo )
488 $ GO TO 380
489*
490* MAIN QZ ITERATION LOOP
491*
492* Initialize dynamic indices
493*
494* Eigenvalues ILAST+1:N have been found.
495* Column operations modify rows IFRSTM:whatever.
496* Row operations modify columns whatever:ILASTM.
497*
498* If only eigenvalues are being computed, then
499* IFRSTM is the row of the last splitting row above row ILAST;
500* this is always at least ILO.
501* IITER counts iterations since the last eigenvalue was found,
502* to tell when to use an extraordinary shift.
503* MAXIT is the maximum number of QZ sweeps allowed.
504*
505 ilast = ihi
506 IF( ilschr ) THEN
507 ifrstm = 1
508 ilastm = n
509 ELSE
510 ifrstm = ilo
511 ilastm = ihi
512 END IF
513 iiter = 0
514 eshift = zero
515 maxit = 30*( ihi-ilo+1 )
516*
517 DO 360 jiter = 1, maxit
518*
519* Split the matrix if possible.
520*
521* Two tests:
522* 1: H(j,j-1)=0 or j=ILO
523* 2: T(j,j)=0
524*
525 IF( ilast.EQ.ilo ) THEN
526*
527* Special case: j=ILAST
528*
529 GO TO 80
530 ELSE
531 IF( abs( h( ilast, ilast-1 ) ).LE.max( safmin, ulp*(
532 $ abs( h( ilast, ilast ) ) + abs( h( ilast-1, ilast-1 ) )
533 $ ) ) ) THEN
534 h( ilast, ilast-1 ) = zero
535 GO TO 80
536 END IF
537 END IF
538*
539 IF( abs( t( ilast, ilast ) ).LE.btol ) THEN
540 t( ilast, ilast ) = zero
541 GO TO 70
542 END IF
543*
544* General case: j<ILAST
545*
546 DO 60 j = ilast - 1, ilo, -1
547*
548* Test 1: for H(j,j-1)=0 or j=ILO
549*
550 IF( j.EQ.ilo ) THEN
551 ilazro = .true.
552 ELSE
553 IF( abs( h( j, j-1 ) ).LE.max( safmin, ulp*(
554 $ abs( h( j, j ) ) + abs( h( j-1, j-1 ) )
555 $ ) ) ) THEN
556 h( j, j-1 ) = zero
557 ilazro = .true.
558 ELSE
559 ilazro = .false.
560 END IF
561 END IF
562*
563* Test 2: for T(j,j)=0
564*
565 IF( abs( t( j, j ) ).LT.btol ) THEN
566 t( j, j ) = zero
567*
568* Test 1a: Check for 2 consecutive small subdiagonals in A
569*
570 ilazr2 = .false.
571 IF( .NOT.ilazro ) THEN
572 temp = abs( h( j, j-1 ) )
573 temp2 = abs( h( j, j ) )
574 tempr = max( temp, temp2 )
575 IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
576 temp = temp / tempr
577 temp2 = temp2 / tempr
578 END IF
579 IF( temp*( ascale*abs( h( j+1, j ) ) ).LE.temp2*
580 $ ( ascale*atol ) )ilazr2 = .true.
581 END IF
582*
583* If both tests pass (1 & 2), i.e., the leading diagonal
584* element of B in the block is zero, split a 1x1 block off
585* at the top. (I.e., at the J-th row/column) The leading
586* diagonal element of the remainder can also be zero, so
587* this may have to be done repeatedly.
588*
589 IF( ilazro .OR. ilazr2 ) THEN
590 DO 40 jch = j, ilast - 1
591 temp = h( jch, jch )
592 CALL dlartg( temp, h( jch+1, jch ), c, s,
593 $ h( jch, jch ) )
594 h( jch+1, jch ) = zero
595 CALL drot( ilastm-jch, h( jch, jch+1 ), ldh,
596 $ h( jch+1, jch+1 ), ldh, c, s )
597 CALL drot( ilastm-jch, t( jch, jch+1 ), ldt,
598 $ t( jch+1, jch+1 ), ldt, c, s )
599 IF( ilq )
600 $ CALL drot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
601 $ c, s )
602 IF( ilazr2 )
603 $ h( jch, jch-1 ) = h( jch, jch-1 )*c
604 ilazr2 = .false.
605 IF( abs( t( jch+1, jch+1 ) ).GE.btol ) THEN
606 IF( jch+1.GE.ilast ) THEN
607 GO TO 80
608 ELSE
609 ifirst = jch + 1
610 GO TO 110
611 END IF
612 END IF
613 t( jch+1, jch+1 ) = zero
614 40 CONTINUE
615 GO TO 70
616 ELSE
617*
618* Only test 2 passed -- chase the zero to T(ILAST,ILAST)
619* Then process as in the case T(ILAST,ILAST)=0
620*
621 DO 50 jch = j, ilast - 1
622 temp = t( jch, jch+1 )
623 CALL dlartg( temp, t( jch+1, jch+1 ), c, s,
624 $ t( jch, jch+1 ) )
625 t( jch+1, jch+1 ) = zero
626 IF( jch.LT.ilastm-1 )
627 $ CALL drot( ilastm-jch-1, t( jch, jch+2 ), ldt,
628 $ t( jch+1, jch+2 ), ldt, c, s )
629 CALL drot( ilastm-jch+2, h( jch, jch-1 ), ldh,
630 $ h( jch+1, jch-1 ), ldh, c, s )
631 IF( ilq )
632 $ CALL drot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
633 $ c, s )
634 temp = h( jch+1, jch )
635 CALL dlartg( temp, h( jch+1, jch-1 ), c, s,
636 $ h( jch+1, jch ) )
637 h( jch+1, jch-1 ) = zero
638 CALL drot( jch+1-ifrstm, h( ifrstm, jch ), 1,
639 $ h( ifrstm, jch-1 ), 1, c, s )
640 CALL drot( jch-ifrstm, t( ifrstm, jch ), 1,
641 $ t( ifrstm, jch-1 ), 1, c, s )
642 IF( ilz )
643 $ CALL drot( n, z( 1, jch ), 1, z( 1, jch-1 ), 1,
644 $ c, s )
645 50 CONTINUE
646 GO TO 70
647 END IF
648 ELSE IF( ilazro ) THEN
649*
650* Only test 1 passed -- work on J:ILAST
651*
652 ifirst = j
653 GO TO 110
654 END IF
655*
656* Neither test passed -- try next J
657*
658 60 CONTINUE
659*
660* (Drop-through is "impossible")
661*
662 info = n + 1
663 GO TO 420
664*
665* T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
666* 1x1 block.
667*
668 70 CONTINUE
669 temp = h( ilast, ilast )
670 CALL dlartg( temp, h( ilast, ilast-1 ), c, s,
671 $ h( ilast, ilast ) )
672 h( ilast, ilast-1 ) = zero
673 CALL drot( ilast-ifrstm, h( ifrstm, ilast ), 1,
674 $ h( ifrstm, ilast-1 ), 1, c, s )
675 CALL drot( ilast-ifrstm, t( ifrstm, ilast ), 1,
676 $ t( ifrstm, ilast-1 ), 1, c, s )
677 IF( ilz )
678 $ CALL drot( n, z( 1, ilast ), 1, z( 1, ilast-1 ), 1, c, s )
679*
680* H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI,
681* and BETA
682*
683 80 CONTINUE
684 IF( t( ilast, ilast ).LT.zero ) THEN
685 IF( ilschr ) THEN
686 DO 90 j = ifrstm, ilast
687 h( j, ilast ) = -h( j, ilast )
688 t( j, ilast ) = -t( j, ilast )
689 90 CONTINUE
690 ELSE
691 h( ilast, ilast ) = -h( ilast, ilast )
692 t( ilast, ilast ) = -t( ilast, ilast )
693 END IF
694 IF( ilz ) THEN
695 DO 100 j = 1, n
696 z( j, ilast ) = -z( j, ilast )
697 100 CONTINUE
698 END IF
699 END IF
700 alphar( ilast ) = h( ilast, ilast )
701 alphai( ilast ) = zero
702 beta( ilast ) = t( ilast, ilast )
703*
704* Go to next block -- exit if finished.
705*
706 ilast = ilast - 1
707 IF( ilast.LT.ilo )
708 $ GO TO 380
709*
710* Reset counters
711*
712 iiter = 0
713 eshift = zero
714 IF( .NOT.ilschr ) THEN
715 ilastm = ilast
716 IF( ifrstm.GT.ilast )
717 $ ifrstm = ilo
718 END IF
719 GO TO 350
720*
721* QZ step
722*
723* This iteration only involves rows/columns IFIRST:ILAST. We
724* assume IFIRST < ILAST, and that the diagonal of B is non-zero.
725*
726 110 CONTINUE
727 iiter = iiter + 1
728 IF( .NOT.ilschr ) THEN
729 ifrstm = ifirst
730 END IF
731*
732* Compute single shifts.
733*
734* At this point, IFIRST < ILAST, and the diagonal elements of
735* T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
736* magnitude)
737*
738 IF( ( iiter / 10 )*10.EQ.iiter ) THEN
739*
740* Exceptional shift. Chosen for no particularly good reason.
741* (Single shift only.)
742*
743 IF( ( dble( maxit )*safmin )*abs( h( ilast, ilast-1 ) ).LT.
744 $ abs( t( ilast-1, ilast-1 ) ) ) THEN
745 eshift = h( ilast, ilast-1 ) /
746 $ t( ilast-1, ilast-1 )
747 ELSE
748 eshift = eshift + one / ( safmin*dble( maxit ) )
749 END IF
750 s1 = one
751 wr = eshift
752*
753 ELSE
754*
755* Shifts based on the generalized eigenvalues of the
756* bottom-right 2x2 block of A and B. The first eigenvalue
757* returned by DLAG2 is the Wilkinson shift (AEP p.512),
758*
759 CALL dlag2( h( ilast-1, ilast-1 ), ldh,
760 $ t( ilast-1, ilast-1 ), ldt, safmin*safety, s1,
761 $ s2, wr, wr2, wi )
762*
763 IF ( abs( (wr/s1)*t( ilast, ilast ) - h( ilast, ilast ) )
764 $ .GT. abs( (wr2/s2)*t( ilast, ilast )
765 $ - h( ilast, ilast ) ) ) THEN
766 temp = wr
767 wr = wr2
768 wr2 = temp
769 temp = s1
770 s1 = s2
771 s2 = temp
772 END IF
773 temp = max( s1, safmin*max( one, abs( wr ), abs( wi ) ) )
774 IF( wi.NE.zero )
775 $ GO TO 200
776 END IF
777*
778* Fiddle with shift to avoid overflow
779*
780 temp = min( ascale, one )*( half*safmax )
781 IF( s1.GT.temp ) THEN
782 scale = temp / s1
783 ELSE
784 scale = one
785 END IF
786*
787 temp = min( bscale, one )*( half*safmax )
788 IF( abs( wr ).GT.temp )
789 $ scale = min( scale, temp / abs( wr ) )
790 s1 = scale*s1
791 wr = scale*wr
792*
793* Now check for two consecutive small subdiagonals.
794*
795 DO 120 j = ilast - 1, ifirst + 1, -1
796 istart = j
797 temp = abs( s1*h( j, j-1 ) )
798 temp2 = abs( s1*h( j, j )-wr*t( j, j ) )
799 tempr = max( temp, temp2 )
800 IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
801 temp = temp / tempr
802 temp2 = temp2 / tempr
803 END IF
804 IF( abs( ( ascale*h( j+1, j ) )*temp ).LE.( ascale*atol )*
805 $ temp2 )GO TO 130
806 120 CONTINUE
807*
808 istart = ifirst
809 130 CONTINUE
810*
811* Do an implicit single-shift QZ sweep.
812*
813* Initial Q
814*
815 temp = s1*h( istart, istart ) - wr*t( istart, istart )
816 temp2 = s1*h( istart+1, istart )
817 CALL dlartg( temp, temp2, c, s, tempr )
818*
819* Sweep
820*
821 DO 190 j = istart, ilast - 1
822 IF( j.GT.istart ) THEN
823 temp = h( j, j-1 )
824 CALL dlartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
825 h( j+1, j-1 ) = zero
826 END IF
827*
828 DO 140 jc = j, ilastm
829 temp = c*h( j, jc ) + s*h( j+1, jc )
830 h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
831 h( j, jc ) = temp
832 temp2 = c*t( j, jc ) + s*t( j+1, jc )
833 t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
834 t( j, jc ) = temp2
835 140 CONTINUE
836 IF( ilq ) THEN
837 DO 150 jr = 1, n
838 temp = c*q( jr, j ) + s*q( jr, j+1 )
839 q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
840 q( jr, j ) = temp
841 150 CONTINUE
842 END IF
843*
844 temp = t( j+1, j+1 )
845 CALL dlartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
846 t( j+1, j ) = zero
847*
848 DO 160 jr = ifrstm, min( j+2, ilast )
849 temp = c*h( jr, j+1 ) + s*h( jr, j )
850 h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
851 h( jr, j+1 ) = temp
852 160 CONTINUE
853 DO 170 jr = ifrstm, j
854 temp = c*t( jr, j+1 ) + s*t( jr, j )
855 t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
856 t( jr, j+1 ) = temp
857 170 CONTINUE
858 IF( ilz ) THEN
859 DO 180 jr = 1, n
860 temp = c*z( jr, j+1 ) + s*z( jr, j )
861 z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
862 z( jr, j+1 ) = temp
863 180 CONTINUE
864 END IF
865 190 CONTINUE
866*
867 GO TO 350
868*
869* Use Francis double-shift
870*
871* Note: the Francis double-shift should work with real shifts,
872* but only if the block is at least 3x3.
873* This code may break if this point is reached with
874* a 2x2 block with real eigenvalues.
875*
876 200 CONTINUE
877 IF( ifirst+1.EQ.ilast ) THEN
878*
879* Special case -- 2x2 block with complex eigenvectors
880*
881* Step 1: Standardize, that is, rotate so that
882*
883* ( B11 0 )
884* B = ( ) with B11 non-negative.
885* ( 0 B22 )
886*
887 CALL dlasv2( t( ilast-1, ilast-1 ), t( ilast-1, ilast ),
888 $ t( ilast, ilast ), b22, b11, sr, cr, sl, cl )
889*
890 IF( b11.LT.zero ) THEN
891 cr = -cr
892 sr = -sr
893 b11 = -b11
894 b22 = -b22
895 END IF
896*
897 CALL drot( ilastm+1-ifirst, h( ilast-1, ilast-1 ), ldh,
898 $ h( ilast, ilast-1 ), ldh, cl, sl )
899 CALL drot( ilast+1-ifrstm, h( ifrstm, ilast-1 ), 1,
900 $ h( ifrstm, ilast ), 1, cr, sr )
901*
902 IF( ilast.LT.ilastm )
903 $ CALL drot( ilastm-ilast, t( ilast-1, ilast+1 ), ldt,
904 $ t( ilast, ilast+1 ), ldt, cl, sl )
905 IF( ifrstm.LT.ilast-1 )
906 $ CALL drot( ifirst-ifrstm, t( ifrstm, ilast-1 ), 1,
907 $ t( ifrstm, ilast ), 1, cr, sr )
908*
909 IF( ilq )
910 $ CALL drot( n, q( 1, ilast-1 ), 1, q( 1, ilast ), 1, cl,
911 $ sl )
912 IF( ilz )
913 $ CALL drot( n, z( 1, ilast-1 ), 1, z( 1, ilast ), 1, cr,
914 $ sr )
915*
916 t( ilast-1, ilast-1 ) = b11
917 t( ilast-1, ilast ) = zero
918 t( ilast, ilast-1 ) = zero
919 t( ilast, ilast ) = b22
920*
921* If B22 is negative, negate column ILAST
922*
923 IF( b22.LT.zero ) THEN
924 DO 210 j = ifrstm, ilast
925 h( j, ilast ) = -h( j, ilast )
926 t( j, ilast ) = -t( j, ilast )
927 210 CONTINUE
928*
929 IF( ilz ) THEN
930 DO 220 j = 1, n
931 z( j, ilast ) = -z( j, ilast )
932 220 CONTINUE
933 END IF
934 b22 = -b22
935 END IF
936*
937* Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
938*
939* Recompute shift
940*
941 CALL dlag2( h( ilast-1, ilast-1 ), ldh,
942 $ t( ilast-1, ilast-1 ), ldt, safmin*safety, s1,
943 $ temp, wr, temp2, wi )
944*
945* If standardization has perturbed the shift onto real line,
946* do another (real single-shift) QR step.
947*
948 IF( wi.EQ.zero )
949 $ GO TO 350
950 s1inv = one / s1
951*
952* Do EISPACK (QZVAL) computation of alpha and beta
953*
954 a11 = h( ilast-1, ilast-1 )
955 a21 = h( ilast, ilast-1 )
956 a12 = h( ilast-1, ilast )
957 a22 = h( ilast, ilast )
958*
959* Compute complex Givens rotation on right
960* (Assume some element of C = (sA - wB) > unfl )
961* __
962* (sA - wB) ( CZ -SZ )
963* ( SZ CZ )
964*
965 c11r = s1*a11 - wr*b11
966 c11i = -wi*b11
967 c12 = s1*a12
968 c21 = s1*a21
969 c22r = s1*a22 - wr*b22
970 c22i = -wi*b22
971*
972 IF( abs( c11r )+abs( c11i )+abs( c12 ).GT.abs( c21 )+
973 $ abs( c22r )+abs( c22i ) ) THEN
974 t1 = dlapy3( c12, c11r, c11i )
975 cz = c12 / t1
976 szr = -c11r / t1
977 szi = -c11i / t1
978 ELSE
979 cz = dlapy2( c22r, c22i )
980 IF( cz.LE.safmin ) THEN
981 cz = zero
982 szr = one
983 szi = zero
984 ELSE
985 tempr = c22r / cz
986 tempi = c22i / cz
987 t1 = dlapy2( cz, c21 )
988 cz = cz / t1
989 szr = -c21*tempr / t1
990 szi = c21*tempi / t1
991 END IF
992 END IF
993*
994* Compute Givens rotation on left
995*
996* ( CQ SQ )
997* ( __ ) A or B
998* ( -SQ CQ )
999*
1000 an = abs( a11 ) + abs( a12 ) + abs( a21 ) + abs( a22 )
1001 bn = abs( b11 ) + abs( b22 )
1002 wabs = abs( wr ) + abs( wi )
1003 IF( s1*an.GT.wabs*bn ) THEN
1004 cq = cz*b11
1005 sqr = szr*b22
1006 sqi = -szi*b22
1007 ELSE
1008 a1r = cz*a11 + szr*a12
1009 a1i = szi*a12
1010 a2r = cz*a21 + szr*a22
1011 a2i = szi*a22
1012 cq = dlapy2( a1r, a1i )
1013 IF( cq.LE.safmin ) THEN
1014 cq = zero
1015 sqr = one
1016 sqi = zero
1017 ELSE
1018 tempr = a1r / cq
1019 tempi = a1i / cq
1020 sqr = tempr*a2r + tempi*a2i
1021 sqi = tempi*a2r - tempr*a2i
1022 END IF
1023 END IF
1024 t1 = dlapy3( cq, sqr, sqi )
1025 cq = cq / t1
1026 sqr = sqr / t1
1027 sqi = sqi / t1
1028*
1029* Compute diagonal elements of QBZ
1030*
1031 tempr = sqr*szr - sqi*szi
1032 tempi = sqr*szi + sqi*szr
1033 b1r = cq*cz*b11 + tempr*b22
1034 b1i = tempi*b22
1035 b1a = dlapy2( b1r, b1i )
1036 b2r = cq*cz*b22 + tempr*b11
1037 b2i = -tempi*b11
1038 b2a = dlapy2( b2r, b2i )
1039*
1040* Normalize so beta > 0, and Im( alpha1 ) > 0
1041*
1042 beta( ilast-1 ) = b1a
1043 beta( ilast ) = b2a
1044 alphar( ilast-1 ) = ( wr*b1a )*s1inv
1045 alphai( ilast-1 ) = ( wi*b1a )*s1inv
1046 alphar( ilast ) = ( wr*b2a )*s1inv
1047 alphai( ilast ) = -( wi*b2a )*s1inv
1048*
1049* Step 3: Go to next block -- exit if finished.
1050*
1051 ilast = ifirst - 1
1052 IF( ilast.LT.ilo )
1053 $ GO TO 380
1054*
1055* Reset counters
1056*
1057 iiter = 0
1058 eshift = zero
1059 IF( .NOT.ilschr ) THEN
1060 ilastm = ilast
1061 IF( ifrstm.GT.ilast )
1062 $ ifrstm = ilo
1063 END IF
1064 GO TO 350
1065 ELSE
1066*
1067* Usual case: 3x3 or larger block, using Francis implicit
1068* double-shift
1069*
1070* 2
1071* Eigenvalue equation is w - c w + d = 0,
1072*
1073* -1 2 -1
1074* so compute 1st column of (A B ) - c A B + d
1075* using the formula in QZIT (from EISPACK)
1076*
1077* We assume that the block is at least 3x3
1078*
1079 ad11 = ( ascale*h( ilast-1, ilast-1 ) ) /
1080 $ ( bscale*t( ilast-1, ilast-1 ) )
1081 ad21 = ( ascale*h( ilast, ilast-1 ) ) /
1082 $ ( bscale*t( ilast-1, ilast-1 ) )
1083 ad12 = ( ascale*h( ilast-1, ilast ) ) /
1084 $ ( bscale*t( ilast, ilast ) )
1085 ad22 = ( ascale*h( ilast, ilast ) ) /
1086 $ ( bscale*t( ilast, ilast ) )
1087 u12 = t( ilast-1, ilast ) / t( ilast, ilast )
1088 ad11l = ( ascale*h( ifirst, ifirst ) ) /
1089 $ ( bscale*t( ifirst, ifirst ) )
1090 ad21l = ( ascale*h( ifirst+1, ifirst ) ) /
1091 $ ( bscale*t( ifirst, ifirst ) )
1092 ad12l = ( ascale*h( ifirst, ifirst+1 ) ) /
1093 $ ( bscale*t( ifirst+1, ifirst+1 ) )
1094 ad22l = ( ascale*h( ifirst+1, ifirst+1 ) ) /
1095 $ ( bscale*t( ifirst+1, ifirst+1 ) )
1096 ad32l = ( ascale*h( ifirst+2, ifirst+1 ) ) /
1097 $ ( bscale*t( ifirst+1, ifirst+1 ) )
1098 u12l = t( ifirst, ifirst+1 ) / t( ifirst+1, ifirst+1 )
1099*
1100 v( 1 ) = ( ad11-ad11l )*( ad22-ad11l ) - ad12*ad21 +
1101 $ ad21*u12*ad11l + ( ad12l-ad11l*u12l )*ad21l
1102 v( 2 ) = ( ( ad22l-ad11l )-ad21l*u12l-( ad11-ad11l )-
1103 $ ( ad22-ad11l )+ad21*u12 )*ad21l
1104 v( 3 ) = ad32l*ad21l
1105*
1106 istart = ifirst
1107*
1108 CALL dlarfg( 3, v( 1 ), v( 2 ), 1, tau )
1109 v( 1 ) = one
1110*
1111* Sweep
1112*
1113 DO 290 j = istart, ilast - 2
1114*
1115* All but last elements: use 3x3 Householder transforms.
1116*
1117* Zero (j-1)st column of A
1118*
1119 IF( j.GT.istart ) THEN
1120 v( 1 ) = h( j, j-1 )
1121 v( 2 ) = h( j+1, j-1 )
1122 v( 3 ) = h( j+2, j-1 )
1123*
1124 CALL dlarfg( 3, h( j, j-1 ), v( 2 ), 1, tau )
1125 v( 1 ) = one
1126 h( j+1, j-1 ) = zero
1127 h( j+2, j-1 ) = zero
1128 END IF
1129*
1130 t2 = tau*v( 2 )
1131 t3 = tau*v( 3 )
1132 DO 230 jc = j, ilastm
1133 temp = h( j, jc )+v( 2 )*h( j+1, jc )+v( 3 )*
1134 $ h( j+2, jc )
1135 h( j, jc ) = h( j, jc ) - temp*tau
1136 h( j+1, jc ) = h( j+1, jc ) - temp*t2
1137 h( j+2, jc ) = h( j+2, jc ) - temp*t3
1138 temp2 = t( j, jc )+v( 2 )*t( j+1, jc )+v( 3 )*
1139 $ t( j+2, jc )
1140 t( j, jc ) = t( j, jc ) - temp2*tau
1141 t( j+1, jc ) = t( j+1, jc ) - temp2*t2
1142 t( j+2, jc ) = t( j+2, jc ) - temp2*t3
1143 230 CONTINUE
1144 IF( ilq ) THEN
1145 DO 240 jr = 1, n
1146 temp = q( jr, j )+v( 2 )*q( jr, j+1 )+v( 3 )*
1147 $ q( jr, j+2 )
1148 q( jr, j ) = q( jr, j ) - temp*tau
1149 q( jr, j+1 ) = q( jr, j+1 ) - temp*t2
1150 q( jr, j+2 ) = q( jr, j+2 ) - temp*t3
1151 240 CONTINUE
1152 END IF
1153*
1154* Zero j-th column of B (see DLAGBC for details)
1155*
1156* Swap rows to pivot
1157*
1158 ilpivt = .false.
1159 temp = max( abs( t( j+1, j+1 ) ), abs( t( j+1, j+2 ) ) )
1160 temp2 = max( abs( t( j+2, j+1 ) ), abs( t( j+2, j+2 ) ) )
1161 IF( max( temp, temp2 ).LT.safmin ) THEN
1162 scale = zero
1163 u1 = one
1164 u2 = zero
1165 GO TO 250
1166 ELSE IF( temp.GE.temp2 ) THEN
1167 w11 = t( j+1, j+1 )
1168 w21 = t( j+2, j+1 )
1169 w12 = t( j+1, j+2 )
1170 w22 = t( j+2, j+2 )
1171 u1 = t( j+1, j )
1172 u2 = t( j+2, j )
1173 ELSE
1174 w21 = t( j+1, j+1 )
1175 w11 = t( j+2, j+1 )
1176 w22 = t( j+1, j+2 )
1177 w12 = t( j+2, j+2 )
1178 u2 = t( j+1, j )
1179 u1 = t( j+2, j )
1180 END IF
1181*
1182* Swap columns if nec.
1183*
1184 IF( abs( w12 ).GT.abs( w11 ) ) THEN
1185 ilpivt = .true.
1186 temp = w12
1187 temp2 = w22
1188 w12 = w11
1189 w22 = w21
1190 w11 = temp
1191 w21 = temp2
1192 END IF
1193*
1194* LU-factor
1195*
1196 temp = w21 / w11
1197 u2 = u2 - temp*u1
1198 w22 = w22 - temp*w12
1199 w21 = zero
1200*
1201* Compute SCALE
1202*
1203 scale = one
1204 IF( abs( w22 ).LT.safmin ) THEN
1205 scale = zero
1206 u2 = one
1207 u1 = -w12 / w11
1208 GO TO 250
1209 END IF
1210 IF( abs( w22 ).LT.abs( u2 ) )
1211 $ scale = abs( w22 / u2 )
1212 IF( abs( w11 ).LT.abs( u1 ) )
1213 $ scale = min( scale, abs( w11 / u1 ) )
1214*
1215* Solve
1216*
1217 u2 = ( scale*u2 ) / w22
1218 u1 = ( scale*u1-w12*u2 ) / w11
1219*
1220 250 CONTINUE
1221 IF( ilpivt ) THEN
1222 temp = u2
1223 u2 = u1
1224 u1 = temp
1225 END IF
1226*
1227* Compute Householder Vector
1228*
1229 t1 = sqrt( scale**2+u1**2+u2**2 )
1230 tau = one + scale / t1
1231 vs = -one / ( scale+t1 )
1232 v( 1 ) = one
1233 v( 2 ) = vs*u1
1234 v( 3 ) = vs*u2
1235*
1236* Apply transformations from the right.
1237*
1238 t2 = tau*v(2)
1239 t3 = tau*v(3)
1240 DO 260 jr = ifrstm, min( j+3, ilast )
1241 temp = h( jr, j )+v( 2 )*h( jr, j+1 )+v( 3 )*
1242 $ h( jr, j+2 )
1243 h( jr, j ) = h( jr, j ) - temp*tau
1244 h( jr, j+1 ) = h( jr, j+1 ) - temp*t2
1245 h( jr, j+2 ) = h( jr, j+2 ) - temp*t3
1246 260 CONTINUE
1247 DO 270 jr = ifrstm, j + 2
1248 temp = t( jr, j )+v( 2 )*t( jr, j+1 )+v( 3 )*
1249 $ t( jr, j+2 )
1250 t( jr, j ) = t( jr, j ) - temp*tau
1251 t( jr, j+1 ) = t( jr, j+1 ) - temp*t2
1252 t( jr, j+2 ) = t( jr, j+2 ) - temp*t3
1253 270 CONTINUE
1254 IF( ilz ) THEN
1255 DO 280 jr = 1, n
1256 temp = z( jr, j )+v( 2 )*z( jr, j+1 )+v( 3 )*
1257 $ z( jr, j+2 )
1258 z( jr, j ) = z( jr, j ) - temp*tau
1259 z( jr, j+1 ) = z( jr, j+1 ) - temp*t2
1260 z( jr, j+2 ) = z( jr, j+2 ) - temp*t3
1261 280 CONTINUE
1262 END IF
1263 t( j+1, j ) = zero
1264 t( j+2, j ) = zero
1265 290 CONTINUE
1266*
1267* Last elements: Use Givens rotations
1268*
1269* Rotations from the left
1270*
1271 j = ilast - 1
1272 temp = h( j, j-1 )
1273 CALL dlartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
1274 h( j+1, j-1 ) = zero
1275*
1276 DO 300 jc = j, ilastm
1277 temp = c*h( j, jc ) + s*h( j+1, jc )
1278 h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
1279 h( j, jc ) = temp
1280 temp2 = c*t( j, jc ) + s*t( j+1, jc )
1281 t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
1282 t( j, jc ) = temp2
1283 300 CONTINUE
1284 IF( ilq ) THEN
1285 DO 310 jr = 1, n
1286 temp = c*q( jr, j ) + s*q( jr, j+1 )
1287 q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
1288 q( jr, j ) = temp
1289 310 CONTINUE
1290 END IF
1291*
1292* Rotations from the right.
1293*
1294 temp = t( j+1, j+1 )
1295 CALL dlartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
1296 t( j+1, j ) = zero
1297*
1298 DO 320 jr = ifrstm, ilast
1299 temp = c*h( jr, j+1 ) + s*h( jr, j )
1300 h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
1301 h( jr, j+1 ) = temp
1302 320 CONTINUE
1303 DO 330 jr = ifrstm, ilast - 1
1304 temp = c*t( jr, j+1 ) + s*t( jr, j )
1305 t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
1306 t( jr, j+1 ) = temp
1307 330 CONTINUE
1308 IF( ilz ) THEN
1309 DO 340 jr = 1, n
1310 temp = c*z( jr, j+1 ) + s*z( jr, j )
1311 z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
1312 z( jr, j+1 ) = temp
1313 340 CONTINUE
1314 END IF
1315*
1316* End of Double-Shift code
1317*
1318 END IF
1319*
1320 GO TO 350
1321*
1322* End of iteration loop
1323*
1324 350 CONTINUE
1325 360 CONTINUE
1326*
1327* Drop-through = non-convergence
1328*
1329 info = ilast
1330 GO TO 420
1331*
1332* Successful completion of all QZ steps
1333*
1334 380 CONTINUE
1335*
1336* Set Eigenvalues 1:ILO-1
1337*
1338 DO 410 j = 1, ilo - 1
1339 IF( t( j, j ).LT.zero ) THEN
1340 IF( ilschr ) THEN
1341 DO 390 jr = 1, j
1342 h( jr, j ) = -h( jr, j )
1343 t( jr, j ) = -t( jr, j )
1344 390 CONTINUE
1345 ELSE
1346 h( j, j ) = -h( j, j )
1347 t( j, j ) = -t( j, j )
1348 END IF
1349 IF( ilz ) THEN
1350 DO 400 jr = 1, n
1351 z( jr, j ) = -z( jr, j )
1352 400 CONTINUE
1353 END IF
1354 END IF
1355 alphar( j ) = h( j, j )
1356 alphai( j ) = zero
1357 beta( j ) = t( j, j )
1358 410 CONTINUE
1359*
1360* Normal Termination
1361*
1362 info = 0
1363*
1364* Exit (other than argument error) -- return optimal workspace size
1365*
1366 420 CONTINUE
1367 work( 1 ) = dble( n )
1368 RETURN
1369*
1370* End of DHGEQZ
1371*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dlag2
subroutine dlag2(a, lda, b, ldb, safmin, scale1, scale2, wr1, wr2, wi)
DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
Definition dlag2.f:156
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlanhs
double precision function dlanhs(norm, n, a, lda, work)
DLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition dlanhs.f:108
dlapy2
double precision function dlapy2(x, y)
DLAPY2 returns sqrt(x2+y2).
Definition dlapy2.f:63
dlapy3
double precision function dlapy3(x, y, z)
DLAPY3 returns sqrt(x2+y2+z2).
Definition dlapy3.f:68
dlarfg
subroutine dlarfg(n, alpha, x, incx, tau)
DLARFG generates an elementary reflector (Householder matrix).
Definition dlarfg.f:106
dlartg
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
dlasv2
subroutine dlasv2(f, g, h, ssmin, ssmax, snr, csr, snl, csl)
DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
Definition dlasv2.f:136
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
Here is the call graph for this function:
Here is the caller graph for this function: