LAPACK 3.12.1
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, , 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 299 of file shgeqz.f.

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