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

◆ shgeqz()

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

SHGEQZ

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

Purpose:
 SHGEQZ 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 SGGHRD.

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