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

◆ dlaqz0()

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

DLAQZ0

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

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

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

 as computed by DGGHRD.

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

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

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

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

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

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

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

 To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently,
 of (A,B)) are computed as a pair of values (alpha,beta), where alpha is
 complex and beta real.
 If beta is nonzero, lambda = alpha / beta is an eigenvalue of the
 generalized nonsymmetric eigenvalue problem (GNEP)
    A*x = lambda*B*x
 and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the
 alternate form of the GNEP
    mu*A*y = B*y.
 Real eigenvalues can be read directly from the generalized Schur
 form:
   alpha = S(i,i), beta = P(i,i).

 Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
      Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
      pp. 241--256.

 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
          The real parts of each scalar alpha defining an eigenvalue
          of GNEP.
[out]ALPHAI
          ALPHAI is DOUBLE PRECISION array, dimension (N)
          The imaginary parts of each scalar alpha defining an
          eigenvalue of GNEP.
          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
          positive, then the j-th and (j+1)-st eigenvalues are a
          complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j).
[out]BETA
          BETA is DOUBLE PRECISION array, dimension (N)
          The scalars beta that define the eigenvalues of GNEP.
          Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
          beta = BETA(j) represent the j-th eigenvalue of the matrix
          pair (A,B), in one of the forms lambda = alpha/beta or
          mu = beta/alpha.  Since either lambda or mu may overflow,
          they should not, in general, be computed.
[in,out]Q
          Q is DOUBLE PRECISION array, dimension (LDQ, N)
          On entry, if COMPQ = 'V', the orthogonal matrix Q1 used in
          the reduction of (A,B) to generalized Hessenberg form.
          On exit, if COMPQ = 'I', the orthogonal matrix of left Schur
          vectors of (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 DOUBLE PRECISION array, dimension (LDZ, N)
          On entry, if COMPZ = 'V', the orthogonal matrix Z1 used in
          the reduction of (A,B) to generalized Hessenberg form.
          On exit, if COMPZ = 'I', the orthogonal matrix of
          right Schur vectors of (H,T), and if COMPZ = 'V', the
          orthogonal matrix of right Schur vectors of (A,B).
          Not referenced if COMPZ = 'N'.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z.  LDZ >= 1.
          If COMPZ='V' or 'I', then LDZ >= N.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
          On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.  LWORK >= max(1,N).

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[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 302 of file dlaqz0.f.

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