LAPACK 3.11.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 $ TAU, TEMP, TEMP2, TEMPI, TEMPR, U1, U12, U12L,
341 $ U2, ULP, VS, W11, W12, W21, W22, WABS, WI, WR,
342 $ WR2
343* ..
344* .. Local Arrays ..
345 REAL V( 3 )
346* ..
347* .. External Functions ..
348 LOGICAL LSAME
349 REAL SLAMCH, SLANHS, SLAPY2, SLAPY3
350 EXTERNAL lsame, slamch, slanhs, slapy2, slapy3
351* ..
352* .. External Subroutines ..
353 EXTERNAL slag2, slarfg, slartg, slaset, slasv2, srot,
354 $ xerbla
355* ..
356* .. Intrinsic Functions ..
357 INTRINSIC abs, max, min, real, 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( 'SHGEQZ', -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 ) = real( 1 )
438 RETURN
439 END IF
440*
441* Initialize Q and Z
442*
443 IF( icompq.EQ.3 )
444 $ CALL slaset( 'Full', n, n, zero, one, q, ldq )
445 IF( icompz.EQ.3 )
446 $ CALL slaset( 'Full', n, n, zero, one, z, ldz )
447*
448* Machine Constants
449*
450 in = ihi + 1 - ilo
451 safmin = slamch( 'S' )
452 safmax = one / safmin
453 ulp = slamch( 'E' )*slamch( 'B' )
454 anorm = slanhs( 'F', in, h( ilo, ilo ), ldh, work )
455 bnorm = slanhs( '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 slartg( temp, h( jch+1, jch ), c, s,
593 $ h( jch, jch ) )
594 h( jch+1, jch ) = zero
595 CALL srot( ilastm-jch, h( jch, jch+1 ), ldh,
596 $ h( jch+1, jch+1 ), ldh, c, s )
597 CALL srot( ilastm-jch, t( jch, jch+1 ), ldt,
598 $ t( jch+1, jch+1 ), ldt, c, s )
599 IF( ilq )
600 $ CALL srot( 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 slartg( 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 srot( ilastm-jch-1, t( jch, jch+2 ), ldt,
628 $ t( jch+1, jch+2 ), ldt, c, s )
629 CALL srot( ilastm-jch+2, h( jch, jch-1 ), ldh,
630 $ h( jch+1, jch-1 ), ldh, c, s )
631 IF( ilq )
632 $ CALL srot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
633 $ c, s )
634 temp = h( jch+1, jch )
635 CALL slartg( temp, h( jch+1, jch-1 ), c, s,
636 $ h( jch+1, jch ) )
637 h( jch+1, jch-1 ) = zero
638 CALL srot( jch+1-ifrstm, h( ifrstm, jch ), 1,
639 $ h( ifrstm, jch-1 ), 1, c, s )
640 CALL srot( jch-ifrstm, t( ifrstm, jch ), 1,
641 $ t( ifrstm, jch-1 ), 1, c, s )
642 IF( ilz )
643 $ CALL srot( 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 slartg( temp, h( ilast, ilast-1 ), c, s,
671 $ h( ilast, ilast ) )
672 h( ilast, ilast-1 ) = zero
673 CALL srot( ilast-ifrstm, h( ifrstm, ilast ), 1,
674 $ h( ifrstm, ilast-1 ), 1, c, s )
675 CALL srot( ilast-ifrstm, t( ifrstm, ilast ), 1,
676 $ t( ifrstm, ilast-1 ), 1, c, s )
677 IF( ilz )
678 $ CALL srot( 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( ( real( 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*real( 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 SLAG2 is the Wilkinson shift (AEP p.512),
758*
759 CALL slag2( 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 slartg( 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 slartg( 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 slartg( 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 slasv2( 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 srot( ilastm+1-ifirst, h( ilast-1, ilast-1 ), ldh,
898 $ h( ilast, ilast-1 ), ldh, cl, sl )
899 CALL srot( ilast+1-ifrstm, h( ifrstm, ilast-1 ), 1,
900 $ h( ifrstm, ilast ), 1, cr, sr )
901*
902 IF( ilast.LT.ilastm )
903 $ CALL srot( ilastm-ilast, t( ilast-1, ilast+1 ), ldt,
904 $ t( ilast, ilast+1 ), ldt, cl, sl )
905 IF( ifrstm.LT.ilast-1 )
906 $ CALL srot( ifirst-ifrstm, t( ifrstm, ilast-1 ), 1,
907 $ t( ifrstm, ilast ), 1, cr, sr )
908*
909 IF( ilq )
910 $ CALL srot( n, q( 1, ilast-1 ), 1, q( 1, ilast ), 1, cl,
911 $ sl )
912 IF( ilz )
913 $ CALL srot( 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 slag2( 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 = slapy3( c12, c11r, c11i )
975 cz = c12 / t1
976 szr = -c11r / t1
977 szi = -c11i / t1
978 ELSE
979 cz = slapy2( 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 = slapy2( 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 = slapy2( 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 = slapy3( 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 = slapy2( b1r, b1i )
1036 b2r = cq*cz*b22 + tempr*b11
1037 b2i = -tempi*b11
1038 b2a = slapy2( 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 slarfg( 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 slarfg( 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 DO 230 jc = j, ilastm
1131 temp = tau*( h( j, jc )+v( 2 )*h( j+1, jc )+v( 3 )*
1132 $ h( j+2, jc ) )
1133 h( j, jc ) = h( j, jc ) - temp
1134 h( j+1, jc ) = h( j+1, jc ) - temp*v( 2 )
1135 h( j+2, jc ) = h( j+2, jc ) - temp*v( 3 )
1136 temp2 = tau*( t( j, jc )+v( 2 )*t( j+1, jc )+v( 3 )*
1137 $ t( j+2, jc ) )
1138 t( j, jc ) = t( j, jc ) - temp2
1139 t( j+1, jc ) = t( j+1, jc ) - temp2*v( 2 )
1140 t( j+2, jc ) = t( j+2, jc ) - temp2*v( 3 )
1141 230 CONTINUE
1142 IF( ilq ) THEN
1143 DO 240 jr = 1, n
1144 temp = tau*( q( jr, j )+v( 2 )*q( jr, j+1 )+v( 3 )*
1145 $ q( jr, j+2 ) )
1146 q( jr, j ) = q( jr, j ) - temp
1147 q( jr, j+1 ) = q( jr, j+1 ) - temp*v( 2 )
1148 q( jr, j+2 ) = q( jr, j+2 ) - temp*v( 3 )
1149 240 CONTINUE
1150 END IF
1151*
1152* Zero j-th column of B (see SLAGBC for details)
1153*
1154* Swap rows to pivot
1155*
1156 ilpivt = .false.
1157 temp = max( abs( t( j+1, j+1 ) ), abs( t( j+1, j+2 ) ) )
1158 temp2 = max( abs( t( j+2, j+1 ) ), abs( t( j+2, j+2 ) ) )
1159 IF( max( temp, temp2 ).LT.safmin ) THEN
1160 scale = zero
1161 u1 = one
1162 u2 = zero
1163 GO TO 250
1164 ELSE IF( temp.GE.temp2 ) THEN
1165 w11 = t( j+1, j+1 )
1166 w21 = t( j+2, j+1 )
1167 w12 = t( j+1, j+2 )
1168 w22 = t( j+2, j+2 )
1169 u1 = t( j+1, j )
1170 u2 = t( j+2, j )
1171 ELSE
1172 w21 = t( j+1, j+1 )
1173 w11 = t( j+2, j+1 )
1174 w22 = t( j+1, j+2 )
1175 w12 = t( j+2, j+2 )
1176 u2 = t( j+1, j )
1177 u1 = t( j+2, j )
1178 END IF
1179*
1180* Swap columns if nec.
1181*
1182 IF( abs( w12 ).GT.abs( w11 ) ) THEN
1183 ilpivt = .true.
1184 temp = w12
1185 temp2 = w22
1186 w12 = w11
1187 w22 = w21
1188 w11 = temp
1189 w21 = temp2
1190 END IF
1191*
1192* LU-factor
1193*
1194 temp = w21 / w11
1195 u2 = u2 - temp*u1
1196 w22 = w22 - temp*w12
1197 w21 = zero
1198*
1199* Compute SCALE
1200*
1201 scale = one
1202 IF( abs( w22 ).LT.safmin ) THEN
1203 scale = zero
1204 u2 = one
1205 u1 = -w12 / w11
1206 GO TO 250
1207 END IF
1208 IF( abs( w22 ).LT.abs( u2 ) )
1209 $ scale = abs( w22 / u2 )
1210 IF( abs( w11 ).LT.abs( u1 ) )
1211 $ scale = min( scale, abs( w11 / u1 ) )
1212*
1213* Solve
1214*
1215 u2 = ( scale*u2 ) / w22
1216 u1 = ( scale*u1-w12*u2 ) / w11
1217*
1218 250 CONTINUE
1219 IF( ilpivt ) THEN
1220 temp = u2
1221 u2 = u1
1222 u1 = temp
1223 END IF
1224*
1225* Compute Householder Vector
1226*
1227 t1 = sqrt( scale**2+u1**2+u2**2 )
1228 tau = one + scale / t1
1229 vs = -one / ( scale+t1 )
1230 v( 1 ) = one
1231 v( 2 ) = vs*u1
1232 v( 3 ) = vs*u2
1233*
1234* Apply transformations from the right.
1235*
1236 DO 260 jr = ifrstm, min( j+3, ilast )
1237 temp = tau*( h( jr, j )+v( 2 )*h( jr, j+1 )+v( 3 )*
1238 $ h( jr, j+2 ) )
1239 h( jr, j ) = h( jr, j ) - temp
1240 h( jr, j+1 ) = h( jr, j+1 ) - temp*v( 2 )
1241 h( jr, j+2 ) = h( jr, j+2 ) - temp*v( 3 )
1242 260 CONTINUE
1243 DO 270 jr = ifrstm, j + 2
1244 temp = tau*( t( jr, j )+v( 2 )*t( jr, j+1 )+v( 3 )*
1245 $ t( jr, j+2 ) )
1246 t( jr, j ) = t( jr, j ) - temp
1247 t( jr, j+1 ) = t( jr, j+1 ) - temp*v( 2 )
1248 t( jr, j+2 ) = t( jr, j+2 ) - temp*v( 3 )
1249 270 CONTINUE
1250 IF( ilz ) THEN
1251 DO 280 jr = 1, n
1252 temp = tau*( z( jr, j )+v( 2 )*z( jr, j+1 )+v( 3 )*
1253 $ z( jr, j+2 ) )
1254 z( jr, j ) = z( jr, j ) - temp
1255 z( jr, j+1 ) = z( jr, j+1 ) - temp*v( 2 )
1256 z( jr, j+2 ) = z( jr, j+2 ) - temp*v( 3 )
1257 280 CONTINUE
1258 END IF
1259 t( j+1, j ) = zero
1260 t( j+2, j ) = zero
1261 290 CONTINUE
1262*
1263* Last elements: Use Givens rotations
1264*
1265* Rotations from the left
1266*
1267 j = ilast - 1
1268 temp = h( j, j-1 )
1269 CALL slartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
1270 h( j+1, j-1 ) = zero
1271*
1272 DO 300 jc = j, ilastm
1273 temp = c*h( j, jc ) + s*h( j+1, jc )
1274 h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
1275 h( j, jc ) = temp
1276 temp2 = c*t( j, jc ) + s*t( j+1, jc )
1277 t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
1278 t( j, jc ) = temp2
1279 300 CONTINUE
1280 IF( ilq ) THEN
1281 DO 310 jr = 1, n
1282 temp = c*q( jr, j ) + s*q( jr, j+1 )
1283 q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
1284 q( jr, j ) = temp
1285 310 CONTINUE
1286 END IF
1287*
1288* Rotations from the right.
1289*
1290 temp = t( j+1, j+1 )
1291 CALL slartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
1292 t( j+1, j ) = zero
1293*
1294 DO 320 jr = ifrstm, ilast
1295 temp = c*h( jr, j+1 ) + s*h( jr, j )
1296 h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
1297 h( jr, j+1 ) = temp
1298 320 CONTINUE
1299 DO 330 jr = ifrstm, ilast - 1
1300 temp = c*t( jr, j+1 ) + s*t( jr, j )
1301 t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
1302 t( jr, j+1 ) = temp
1303 330 CONTINUE
1304 IF( ilz ) THEN
1305 DO 340 jr = 1, n
1306 temp = c*z( jr, j+1 ) + s*z( jr, j )
1307 z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
1308 z( jr, j+1 ) = temp
1309 340 CONTINUE
1310 END IF
1311*
1312* End of Double-Shift code
1313*
1314 END IF
1315*
1316 GO TO 350
1317*
1318* End of iteration loop
1319*
1320 350 CONTINUE
1321 360 CONTINUE
1322*
1323* Drop-through = non-convergence
1324*
1325 info = ilast
1326 GO TO 420
1327*
1328* Successful completion of all QZ steps
1329*
1330 380 CONTINUE
1331*
1332* Set Eigenvalues 1:ILO-1
1333*
1334 DO 410 j = 1, ilo - 1
1335 IF( t( j, j ).LT.zero ) THEN
1336 IF( ilschr ) THEN
1337 DO 390 jr = 1, j
1338 h( jr, j ) = -h( jr, j )
1339 t( jr, j ) = -t( jr, j )
1340 390 CONTINUE
1341 ELSE
1342 h( j, j ) = -h( j, j )
1343 t( j, j ) = -t( j, j )
1344 END IF
1345 IF( ilz ) THEN
1346 DO 400 jr = 1, n
1347 z( jr, j ) = -z( jr, j )
1348 400 CONTINUE
1349 END IF
1350 END IF
1351 alphar( j ) = h( j, j )
1352 alphai( j ) = zero
1353 beta( j ) = t( j, j )
1354 410 CONTINUE
1355*
1356* Normal Termination
1357*
1358 info = 0
1359*
1360* Exit (other than argument error) -- return optimal workspace size
1361*
1362 420 CONTINUE
1363 work( 1 ) = real( n )
1364 RETURN
1365*
1366* End of SHGEQZ
1367*
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:138
slartg
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f90:111
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
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
slarfg
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
Definition: slarfg.f:106
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
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
srot
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:92
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: