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

◆ slaqz0()

recursive subroutine slaqz0 ( character, intent(in)  WANTS,
character, intent(in)  WANTQ,
character, intent(in)  WANTZ,
integer, intent(in)  N,
integer, intent(in)  ILO,
integer, intent(in)  IHI,
real, dimension( lda, * ), intent(inout)  A,
integer, intent(in)  LDA,
real, dimension( ldb, * ), intent(inout)  B,
integer, intent(in)  LDB,
real, dimension( * ), intent(inout)  ALPHAR,
real, dimension( * ), intent(inout)  ALPHAI,
real, dimension( * ), intent(inout)  BETA,
real, dimension( ldq, * ), intent(inout)  Q,
integer, intent(in)  LDQ,
real, dimension( ldz, * ), intent(inout)  Z,
integer, intent(in)  LDZ,
real, dimension( * ), intent(inout)  WORK,
integer, intent(in)  LWORK,
integer, intent(in)  REC,
integer, intent(out)  INFO 
)

SLAQZ0

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

Purpose:
 SLAQZ0 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.

 Ref: B. Kagstrom, D. Kressner, "Multishift Variants of the QZ
      Algorithm with Aggressive Early Deflation", SIAM J. Numer.
      Anal., 29(2006), pp. 199--227.

 Ref: T. Steel, D. Camps, K. Meerbergen, R. Vandebril "A multishift,
      multipole rational QZ method with agressive early deflation"
Parameters
[in]WANTS
          WANTS is CHARACTER*1
          = 'E': Compute eigenvalues only;
          = 'S': Compute eigenvalues and the Schur form.
[in]WANTQ
          WANTQ 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 (A,B) is returned;
          = 'V': Q must contain an orthogonal matrix Q1 on entry and
                 the product Q1*Q is returned.
[in]WANTZ
          WANTZ 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 (A,B) 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 A, B, Q, and Z.  N >= 0.
[in]ILO
          ILO is INTEGER
[in]IHI
          IHI is INTEGER
          ILO and IHI mark the rows and columns of A 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]A
          A is REAL array, dimension (LDA, N)
          On entry, the N-by-N upper Hessenberg matrix A.
          On exit, if JOB = 'S', A contains the upper quasi-triangular
          matrix S from the generalized Schur factorization.
          If JOB = 'E', the diagonal blocks of A match those of S, but
          the rest of A is unspecified.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max( 1, N ).
[in,out]B
          B is REAL array, dimension (LDB, N)
          On entry, the N-by-N upper triangular matrix B.
          On exit, if JOB = 'S', B 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 A(j+1,j) is
          non-zero, then B(j+1,j) = B(j,j+1) = 0, B(j,j) > 0, and
          B(j+1,j+1) > 0.
          If JOB = 'E', the diagonal blocks of B match those of P, but
          the rest of B is unspecified.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= 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 (A,B), 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.
[in]REC
          REC is INTEGER
             REC indicates the current recursion level. Should be set
             to 0 on first call.
[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.  (A,B) is not
                     in Schur form, but ALPHAR(i), ALPHAI(i), and
                     BETA(i), i=INFO+1,...,N should be correct.
Author
Thijs Steel, KU Leuven
Date
May 2020

Definition at line 300 of file slaqz0.f.

304 IMPLICIT NONE
305
306* Arguments
307 CHARACTER, INTENT( IN ) :: WANTS, WANTQ, WANTZ
308 INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
309 $ REC
310
311 INTEGER, INTENT( OUT ) :: INFO
312
313 REAL, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
314 $ Z( LDZ, * ), ALPHAR( * ), ALPHAI( * ), BETA( * ), WORK( * )
315
316* Parameters
317 REAL :: ZERO, ONE, HALF
318 parameter( zero = 0.0, one = 1.0, half = 0.5 )
319
320* Local scalars
321 REAL :: SMLNUM, ULP, ESHIFT, SAFMIN, SAFMAX, C1, S1, TEMP, SWAP,
322 $ BNORM, BTOL
323 INTEGER :: ISTART, ISTOP, IITER, MAXIT, ISTART2, K, LD, NSHIFTS,
324 $ NBLOCK, NW, NMIN, NIBBLE, N_UNDEFLATED, N_DEFLATED,
325 $ NS, SWEEP_INFO, SHIFTPOS, LWORKREQ, K2, ISTARTM,
326 $ ISTOPM, IWANTS, IWANTQ, IWANTZ, NORM_INFO, AED_INFO,
327 $ NWR, NBR, NSR, ITEMP1, ITEMP2, RCOST, I
328 LOGICAL :: ILSCHUR, ILQ, ILZ
329 CHARACTER :: JBCMPZ*3
330
331* External Functions
332 EXTERNAL :: xerbla, shgeqz, slaqz3, slaqz4, slaset, slabad,
333 $ slartg, srot
334 REAL, EXTERNAL :: SLAMCH, SLANHS
335 LOGICAL, EXTERNAL :: LSAME
336 INTEGER, EXTERNAL :: ILAENV
337
338*
339* Decode wantS,wantQ,wantZ
340*
341 IF( lsame( wants, 'E' ) ) THEN
342 ilschur = .false.
343 iwants = 1
344 ELSE IF( lsame( wants, 'S' ) ) THEN
345 ilschur = .true.
346 iwants = 2
347 ELSE
348 iwants = 0
349 END IF
350
351 IF( lsame( wantq, 'N' ) ) THEN
352 ilq = .false.
353 iwantq = 1
354 ELSE IF( lsame( wantq, 'V' ) ) THEN
355 ilq = .true.
356 iwantq = 2
357 ELSE IF( lsame( wantq, 'I' ) ) THEN
358 ilq = .true.
359 iwantq = 3
360 ELSE
361 iwantq = 0
362 END IF
363
364 IF( lsame( wantz, 'N' ) ) THEN
365 ilz = .false.
366 iwantz = 1
367 ELSE IF( lsame( wantz, 'V' ) ) THEN
368 ilz = .true.
369 iwantz = 2
370 ELSE IF( lsame( wantz, 'I' ) ) THEN
371 ilz = .true.
372 iwantz = 3
373 ELSE
374 iwantz = 0
375 END IF
376*
377* Check Argument Values
378*
379 info = 0
380 IF( iwants.EQ.0 ) THEN
381 info = -1
382 ELSE IF( iwantq.EQ.0 ) THEN
383 info = -2
384 ELSE IF( iwantz.EQ.0 ) THEN
385 info = -3
386 ELSE IF( n.LT.0 ) THEN
387 info = -4
388 ELSE IF( ilo.LT.1 ) THEN
389 info = -5
390 ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 ) THEN
391 info = -6
392 ELSE IF( lda.LT.n ) THEN
393 info = -8
394 ELSE IF( ldb.LT.n ) THEN
395 info = -10
396 ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) ) THEN
397 info = -15
398 ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) ) THEN
399 info = -17
400 END IF
401 IF( info.NE.0 ) THEN
402 CALL xerbla( 'SLAQZ0', -info )
403 RETURN
404 END IF
405
406*
407* Quick return if possible
408*
409 IF( n.LE.0 ) THEN
410 work( 1 ) = real( 1 )
411 RETURN
412 END IF
413
414*
415* Get the parameters
416*
417 jbcmpz( 1:1 ) = wants
418 jbcmpz( 2:2 ) = wantq
419 jbcmpz( 3:3 ) = wantz
420
421 nmin = ilaenv( 12, 'SLAQZ0', jbcmpz, n, ilo, ihi, lwork )
422
423 nwr = ilaenv( 13, 'SLAQZ0', jbcmpz, n, ilo, ihi, lwork )
424 nwr = max( 2, nwr )
425 nwr = min( ihi-ilo+1, ( n-1 ) / 3, nwr )
426
427 nibble = ilaenv( 14, 'SLAQZ0', jbcmpz, n, ilo, ihi, lwork )
428
429 nsr = ilaenv( 15, 'SLAQZ0', jbcmpz, n, ilo, ihi, lwork )
430 nsr = min( nsr, ( n+6 ) / 9, ihi-ilo )
431 nsr = max( 2, nsr-mod( nsr, 2 ) )
432
433 rcost = ilaenv( 17, 'SLAQZ0', jbcmpz, n, ilo, ihi, lwork )
434 itemp1 = int( nsr/sqrt( 1+2*nsr/( real( rcost )/100*n ) ) )
435 itemp1 = ( ( itemp1-1 )/4 )*4+4
436 nbr = nsr+itemp1
437
438 IF( n .LT. nmin .OR. rec .GE. 2 ) THEN
439 CALL shgeqz( wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb,
440 $ alphar, alphai, beta, q, ldq, z, ldz, work,
441 $ lwork, info )
442 RETURN
443 END IF
444
445*
446* Find out required workspace
447*
448
449* Workspace query to slaqz3
450 nw = max( nwr, nmin )
451 CALL slaqz3( ilschur, ilq, ilz, n, ilo, ihi, nw, a, lda, b, ldb,
452 $ q, ldq, z, ldz, n_undeflated, n_deflated, alphar,
453 $ alphai, beta, work, nw, work, nw, work, -1, rec,
454 $ aed_info )
455 itemp1 = int( work( 1 ) )
456* Workspace query to slaqz4
457 CALL slaqz4( ilschur, ilq, ilz, n, ilo, ihi, nsr, nbr, alphar,
458 $ alphai, beta, a, lda, b, ldb, q, ldq, z, ldz, work,
459 $ nbr, work, nbr, work, -1, sweep_info )
460 itemp2 = int( work( 1 ) )
461
462 lworkreq = max( itemp1+2*nw**2, itemp2+2*nbr**2 )
463 IF ( lwork .EQ.-1 ) THEN
464 work( 1 ) = real( lworkreq )
465 RETURN
466 ELSE IF ( lwork .LT. lworkreq ) THEN
467 info = -19
468 END IF
469 IF( info.NE.0 ) THEN
470 CALL xerbla( 'SLAQZ0', info )
471 RETURN
472 END IF
473*
474* Initialize Q and Z
475*
476 IF( iwantq.EQ.3 ) CALL slaset( 'FULL', n, n, zero, one, q, ldq )
477 IF( iwantz.EQ.3 ) CALL slaset( 'FULL', n, n, zero, one, z, ldz )
478
479* Get machine constants
480 safmin = slamch( 'SAFE MINIMUM' )
481 safmax = one/safmin
482 CALL slabad( safmin, safmax )
483 ulp = slamch( 'PRECISION' )
484 smlnum = safmin*( real( n )/ulp )
485
486 bnorm = slanhs( 'F', ihi-ilo+1, b( ilo, ilo ), ldb, work )
487 btol = max( safmin, ulp*bnorm )
488
489 istart = ilo
490 istop = ihi
491 maxit = 3*( ihi-ilo+1 )
492 ld = 0
493
494 DO iiter = 1, maxit
495 IF( iiter .GE. maxit ) THEN
496 info = istop+1
497 GOTO 80
498 END IF
499 IF ( istart+1 .GE. istop ) THEN
500 istop = istart
501 EXIT
502 END IF
503
504* Check deflations at the end
505 IF ( abs( a( istop-1, istop-2 ) ) .LE. max( smlnum,
506 $ ulp*( abs( a( istop-1, istop-1 ) )+abs( a( istop-2,
507 $ istop-2 ) ) ) ) ) THEN
508 a( istop-1, istop-2 ) = zero
509 istop = istop-2
510 ld = 0
511 eshift = zero
512 ELSE IF ( abs( a( istop, istop-1 ) ) .LE. max( smlnum,
513 $ ulp*( abs( a( istop, istop ) )+abs( a( istop-1,
514 $ istop-1 ) ) ) ) ) THEN
515 a( istop, istop-1 ) = zero
516 istop = istop-1
517 ld = 0
518 eshift = zero
519 END IF
520* Check deflations at the start
521 IF ( abs( a( istart+2, istart+1 ) ) .LE. max( smlnum,
522 $ ulp*( abs( a( istart+1, istart+1 ) )+abs( a( istart+2,
523 $ istart+2 ) ) ) ) ) THEN
524 a( istart+2, istart+1 ) = zero
525 istart = istart+2
526 ld = 0
527 eshift = zero
528 ELSE IF ( abs( a( istart+1, istart ) ) .LE. max( smlnum,
529 $ ulp*( abs( a( istart, istart ) )+abs( a( istart+1,
530 $ istart+1 ) ) ) ) ) THEN
531 a( istart+1, istart ) = zero
532 istart = istart+1
533 ld = 0
534 eshift = zero
535 END IF
536
537 IF ( istart+1 .GE. istop ) THEN
538 EXIT
539 END IF
540
541* Check interior deflations
542 istart2 = istart
543 DO k = istop, istart+1, -1
544 IF ( abs( a( k, k-1 ) ) .LE. max( smlnum, ulp*( abs( a( k,
545 $ k ) )+abs( a( k-1, k-1 ) ) ) ) ) THEN
546 a( k, k-1 ) = zero
547 istart2 = k
548 EXIT
549 END IF
550 END DO
551
552* Get range to apply rotations to
553 IF ( ilschur ) THEN
554 istartm = 1
555 istopm = n
556 ELSE
557 istartm = istart2
558 istopm = istop
559 END IF
560
561* Check infinite eigenvalues, this is done without blocking so might
562* slow down the method when many infinite eigenvalues are present
563 k = istop
564 DO WHILE ( k.GE.istart2 )
565
566 IF( abs( b( k, k ) ) .LT. btol ) THEN
567* A diagonal element of B is negligable, move it
568* to the top and deflate it
569
570 DO k2 = k, istart2+1, -1
571 CALL slartg( b( k2-1, k2 ), b( k2-1, k2-1 ), c1, s1,
572 $ temp )
573 b( k2-1, k2 ) = temp
574 b( k2-1, k2-1 ) = zero
575
576 CALL srot( k2-2-istartm+1, b( istartm, k2 ), 1,
577 $ b( istartm, k2-1 ), 1, c1, s1 )
578 CALL srot( min( k2+1, istop )-istartm+1, a( istartm,
579 $ k2 ), 1, a( istartm, k2-1 ), 1, c1, s1 )
580 IF ( ilz ) THEN
581 CALL srot( n, z( 1, k2 ), 1, z( 1, k2-1 ), 1, c1,
582 $ s1 )
583 END IF
584
585 IF( k2.LT.istop ) THEN
586 CALL slartg( a( k2, k2-1 ), a( k2+1, k2-1 ), c1,
587 $ s1, temp )
588 a( k2, k2-1 ) = temp
589 a( k2+1, k2-1 ) = zero
590
591 CALL srot( istopm-k2+1, a( k2, k2 ), lda, a( k2+1,
592 $ k2 ), lda, c1, s1 )
593 CALL srot( istopm-k2+1, b( k2, k2 ), ldb, b( k2+1,
594 $ k2 ), ldb, c1, s1 )
595 IF( ilq ) THEN
596 CALL srot( n, q( 1, k2 ), 1, q( 1, k2+1 ), 1,
597 $ c1, s1 )
598 END IF
599 END IF
600
601 END DO
602
603 IF( istart2.LT.istop )THEN
604 CALL slartg( a( istart2, istart2 ), a( istart2+1,
605 $ istart2 ), c1, s1, temp )
606 a( istart2, istart2 ) = temp
607 a( istart2+1, istart2 ) = zero
608
609 CALL srot( istopm-( istart2+1 )+1, a( istart2,
610 $ istart2+1 ), lda, a( istart2+1,
611 $ istart2+1 ), lda, c1, s1 )
612 CALL srot( istopm-( istart2+1 )+1, b( istart2,
613 $ istart2+1 ), ldb, b( istart2+1,
614 $ istart2+1 ), ldb, c1, s1 )
615 IF( ilq ) THEN
616 CALL srot( n, q( 1, istart2 ), 1, q( 1,
617 $ istart2+1 ), 1, c1, s1 )
618 END IF
619 END IF
620
621 istart2 = istart2+1
622
623 END IF
624 k = k-1
625 END DO
626
627* istart2 now points to the top of the bottom right
628* unreduced Hessenberg block
629 IF ( istart2 .GE. istop ) THEN
630 istop = istart2-1
631 ld = 0
632 eshift = zero
633 cycle
634 END IF
635
636 nw = nwr
637 nshifts = nsr
638 nblock = nbr
639
640 IF ( istop-istart2+1 .LT. nmin ) THEN
641* Setting nw to the size of the subblock will make AED deflate
642* all the eigenvalues. This is slightly more efficient than just
643* using qz_small because the off diagonal part gets updated via BLAS.
644 IF ( istop-istart+1 .LT. nmin ) THEN
645 nw = istop-istart+1
646 istart2 = istart
647 ELSE
648 nw = istop-istart2+1
649 END IF
650 END IF
651
652*
653* Time for AED
654*
655 CALL slaqz3( ilschur, ilq, ilz, n, istart2, istop, nw, a, lda,
656 $ b, ldb, q, ldq, z, ldz, n_undeflated, n_deflated,
657 $ alphar, alphai, beta, work, nw, work( nw**2+1 ),
658 $ nw, work( 2*nw**2+1 ), lwork-2*nw**2, rec,
659 $ aed_info )
660
661 IF ( n_deflated > 0 ) THEN
662 istop = istop-n_deflated
663 ld = 0
664 eshift = zero
665 END IF
666
667 IF ( 100*n_deflated > nibble*( n_deflated+n_undeflated ) .OR.
668 $ istop-istart2+1 .LT. nmin ) THEN
669* AED has uncovered many eigenvalues. Skip a QZ sweep and run
670* AED again.
671 cycle
672 END IF
673
674 ld = ld+1
675
676 ns = min( nshifts, istop-istart2 )
677 ns = min( ns, n_undeflated )
678 shiftpos = istop-n_deflated-n_undeflated+1
679*
680* Shuffle shifts to put double shifts in front
681* This ensures that we don't split up a double shift
682*
683 DO i = shiftpos, shiftpos+n_undeflated-1, 2
684 IF( alphai( i ).NE.-alphai( i+1 ) ) THEN
685*
686 swap = alphar( i )
687 alphar( i ) = alphar( i+1 )
688 alphar( i+1 ) = alphar( i+2 )
689 alphar( i+2 ) = swap
690
691 swap = alphai( i )
692 alphai( i ) = alphai( i+1 )
693 alphai( i+1 ) = alphai( i+2 )
694 alphai( i+2 ) = swap
695
696 swap = beta( i )
697 beta( i ) = beta( i+1 )
698 beta( i+1 ) = beta( i+2 )
699 beta( i+2 ) = swap
700 END IF
701 END DO
702
703 IF ( mod( ld, 6 ) .EQ. 0 ) THEN
704*
705* Exceptional shift. Chosen for no particularly good reason.
706*
707 IF( ( real( maxit )*safmin )*abs( a( istop,
708 $ istop-1 ) ).LT.abs( a( istop-1, istop-1 ) ) ) THEN
709 eshift = a( istop, istop-1 )/b( istop-1, istop-1 )
710 ELSE
711 eshift = eshift+one/( safmin*real( maxit ) )
712 END IF
713 alphar( shiftpos ) = one
714 alphar( shiftpos+1 ) = zero
715 alphai( shiftpos ) = zero
716 alphai( shiftpos+1 ) = zero
717 beta( shiftpos ) = eshift
718 beta( shiftpos+1 ) = eshift
719 ns = 2
720 END IF
721
722*
723* Time for a QZ sweep
724*
725 CALL slaqz4( ilschur, ilq, ilz, n, istart2, istop, ns, nblock,
726 $ alphar( shiftpos ), alphai( shiftpos ),
727 $ beta( shiftpos ), a, lda, b, ldb, q, ldq, z, ldz,
728 $ work, nblock, work( nblock**2+1 ), nblock,
729 $ work( 2*nblock**2+1 ), lwork-2*nblock**2,
730 $ sweep_info )
731
732 END DO
733
734*
735* Call SHGEQZ to normalize the eigenvalue blocks and set the eigenvalues
736* If all the eigenvalues have been found, SHGEQZ will not do any iterations
737* and only normalize the blocks. In case of a rare convergence failure,
738* the single shift might perform better.
739*
740 80 CALL shgeqz( wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb,
741 $ alphar, alphai, beta, q, ldq, z, ldz, work, lwork,
742 $ norm_info )
743
744 info = norm_info
745
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
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
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f90:111
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine slaqz4(ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS, NBLOCK_DESIRED, SR, SI, SS, A, LDA, B, LDB, Q, LDQ, Z, LDZ, QC, LDQC, ZC, LDZC, WORK, LWORK, INFO)
SLAQZ4
Definition: slaqz4.f:214
recursive subroutine slaqz3(ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW, A, LDA, B, LDB, Q, LDQ, Z, LDZ, NS, ND, ALPHAR, ALPHAI, BETA, QC, LDQC, ZC, LDZC, WORK, LWORK, REC, INFO)
SLAQZ3
Definition: slaqz3.f:238
subroutine shgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
SHGEQZ
Definition: shgeqz.f:304
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
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:92
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: