LAPACK 3.12.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 aggressive 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,
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 ulp = dlamch( 'PRECISION' )
486 smlnum = safmin*( dble( n )/ulp )
487
488 bnorm = dlanhs( 'F', ihi-ilo+1, b( ilo, ilo ), ldb, work )
489 btol = max( safmin, ulp*bnorm )
490
491 istart = ilo
492 istop = ihi
493 maxit = 3*( ihi-ilo+1 )
494 ld = 0
495
496 DO iiter = 1, maxit
497 IF( iiter .GE. maxit ) THEN
498 info = istop+1
499 GOTO 80
500 END IF
501 IF ( istart+1 .GE. istop ) THEN
502 istop = istart
503 EXIT
504 END IF
505
506* Check deflations at the end
507 IF ( abs( a( istop-1, istop-2 ) ) .LE. max( smlnum,
508 $ ulp*( abs( a( istop-1, istop-1 ) )+abs( a( istop-2,
509 $ istop-2 ) ) ) ) ) THEN
510 a( istop-1, istop-2 ) = zero
511 istop = istop-2
512 ld = 0
513 eshift = zero
514 ELSE IF ( abs( a( istop, istop-1 ) ) .LE. max( smlnum,
515 $ ulp*( abs( a( istop, istop ) )+abs( a( istop-1,
516 $ istop-1 ) ) ) ) ) THEN
517 a( istop, istop-1 ) = zero
518 istop = istop-1
519 ld = 0
520 eshift = zero
521 END IF
522* Check deflations at the start
523 IF ( abs( a( istart+2, istart+1 ) ) .LE. max( smlnum,
524 $ ulp*( abs( a( istart+1, istart+1 ) )+abs( a( istart+2,
525 $ istart+2 ) ) ) ) ) THEN
526 a( istart+2, istart+1 ) = zero
527 istart = istart+2
528 ld = 0
529 eshift = zero
530 ELSE IF ( abs( a( istart+1, istart ) ) .LE. max( smlnum,
531 $ ulp*( abs( a( istart, istart ) )+abs( a( istart+1,
532 $ istart+1 ) ) ) ) ) THEN
533 a( istart+1, istart ) = zero
534 istart = istart+1
535 ld = 0
536 eshift = zero
537 END IF
538
539 IF ( istart+1 .GE. istop ) THEN
540 EXIT
541 END IF
542
543* Check interior deflations
544 istart2 = istart
545 DO k = istop, istart+1, -1
546 IF ( abs( a( k, k-1 ) ) .LE. max( smlnum, ulp*( abs( a( k,
547 $ k ) )+abs( a( k-1, k-1 ) ) ) ) ) THEN
548 a( k, k-1 ) = zero
549 istart2 = k
550 EXIT
551 END IF
552 END DO
553
554* Get range to apply rotations to
555 IF ( ilschur ) THEN
556 istartm = 1
557 istopm = n
558 ELSE
559 istartm = istart2
560 istopm = istop
561 END IF
562
563* Check infinite eigenvalues, this is done without blocking so might
564* slow down the method when many infinite eigenvalues are present
565 k = istop
566 DO WHILE ( k.GE.istart2 )
567
568 IF( abs( b( k, k ) ) .LT. btol ) THEN
569* A diagonal element of B is negligible, move it
570* to the top and deflate it
571
572 DO k2 = k, istart2+1, -1
573 CALL dlartg( b( k2-1, k2 ), b( k2-1, k2-1 ), c1, s1,
574 $ temp )
575 b( k2-1, k2 ) = temp
576 b( k2-1, k2-1 ) = zero
577
578 CALL drot( k2-2-istartm+1, b( istartm, k2 ), 1,
579 $ b( istartm, k2-1 ), 1, c1, s1 )
580 CALL drot( min( k2+1, istop )-istartm+1, a( istartm,
581 $ k2 ), 1, a( istartm, k2-1 ), 1, c1, s1 )
582 IF ( ilz ) THEN
583 CALL drot( n, z( 1, k2 ), 1, z( 1, k2-1 ), 1, c1,
584 $ s1 )
585 END IF
586
587 IF( k2.LT.istop ) THEN
588 CALL dlartg( a( k2, k2-1 ), a( k2+1, k2-1 ), c1,
589 $ s1, temp )
590 a( k2, k2-1 ) = temp
591 a( k2+1, k2-1 ) = zero
592
593 CALL drot( istopm-k2+1, a( k2, k2 ), lda, a( k2+1,
594 $ k2 ), lda, c1, s1 )
595 CALL drot( istopm-k2+1, b( k2, k2 ), ldb, b( k2+1,
596 $ k2 ), ldb, c1, s1 )
597 IF( ilq ) THEN
598 CALL drot( n, q( 1, k2 ), 1, q( 1, k2+1 ), 1,
599 $ c1, s1 )
600 END IF
601 END IF
602
603 END DO
604
605 IF( istart2.LT.istop )THEN
606 CALL dlartg( a( istart2, istart2 ), a( istart2+1,
607 $ istart2 ), c1, s1, temp )
608 a( istart2, istart2 ) = temp
609 a( istart2+1, istart2 ) = zero
610
611 CALL drot( istopm-( istart2+1 )+1, a( istart2,
612 $ istart2+1 ), lda, a( istart2+1,
613 $ istart2+1 ), lda, c1, s1 )
614 CALL drot( istopm-( istart2+1 )+1, b( istart2,
615 $ istart2+1 ), ldb, b( istart2+1,
616 $ istart2+1 ), ldb, c1, s1 )
617 IF( ilq ) THEN
618 CALL drot( n, q( 1, istart2 ), 1, q( 1,
619 $ istart2+1 ), 1, c1, s1 )
620 END IF
621 END IF
622
623 istart2 = istart2+1
624
625 END IF
626 k = k-1
627 END DO
628
629* istart2 now points to the top of the bottom right
630* unreduced Hessenberg block
631 IF ( istart2 .GE. istop ) THEN
632 istop = istart2-1
633 ld = 0
634 eshift = zero
635 cycle
636 END IF
637
638 nw = nwr
639 nshifts = nsr
640 nblock = nbr
641
642 IF ( istop-istart2+1 .LT. nmin ) THEN
643* Setting nw to the size of the subblock will make AED deflate
644* all the eigenvalues. This is slightly more efficient than just
645* using DHGEQZ because the off diagonal part gets updated via BLAS.
646 IF ( istop-istart+1 .LT. nmin ) THEN
647 nw = istop-istart+1
648 istart2 = istart
649 ELSE
650 nw = istop-istart2+1
651 END IF
652 END IF
653
654*
655* Time for AED
656*
657 CALL dlaqz3( ilschur, ilq, ilz, n, istart2, istop, nw, a, lda,
658 $ b, ldb, q, ldq, z, ldz, n_undeflated, n_deflated,
659 $ alphar, alphai, beta, work, nw, work( nw**2+1 ),
660 $ nw, work( 2*nw**2+1 ), lwork-2*nw**2, rec,
661 $ aed_info )
662
663 IF ( n_deflated > 0 ) THEN
664 istop = istop-n_deflated
665 ld = 0
666 eshift = zero
667 END IF
668
669 IF ( 100*n_deflated > nibble*( n_deflated+n_undeflated ) .OR.
670 $ istop-istart2+1 .LT. nmin ) THEN
671* AED has uncovered many eigenvalues. Skip a QZ sweep and run
672* AED again.
673 cycle
674 END IF
675
676 ld = ld+1
677
678 ns = min( nshifts, istop-istart2 )
679 ns = min( ns, n_undeflated )
680 shiftpos = istop-n_undeflated+1
681*
682* Shuffle shifts to put double shifts in front
683* This ensures that we don't split up a double shift
684*
685 DO i = shiftpos, shiftpos+n_undeflated-1, 2
686 IF( alphai( i ).NE.-alphai( i+1 ) ) THEN
687*
688 swap = alphar( i )
689 alphar( i ) = alphar( i+1 )
690 alphar( i+1 ) = alphar( i+2 )
691 alphar( i+2 ) = swap
692
693 swap = alphai( i )
694 alphai( i ) = alphai( i+1 )
695 alphai( i+1 ) = alphai( i+2 )
696 alphai( i+2 ) = swap
697
698 swap = beta( i )
699 beta( i ) = beta( i+1 )
700 beta( i+1 ) = beta( i+2 )
701 beta( i+2 ) = swap
702 END IF
703 END DO
704
705 IF ( mod( ld, 6 ) .EQ. 0 ) THEN
706*
707* Exceptional shift. Chosen for no particularly good reason.
708*
709 IF( ( dble( maxit )*safmin )*abs( a( istop,
710 $ istop-1 ) ).LT.abs( a( istop-1, istop-1 ) ) ) THEN
711 eshift = a( istop, istop-1 )/b( istop-1, istop-1 )
712 ELSE
713 eshift = eshift+one/( safmin*dble( maxit ) )
714 END IF
715 alphar( shiftpos ) = one
716 alphar( shiftpos+1 ) = zero
717 alphai( shiftpos ) = zero
718 alphai( shiftpos+1 ) = zero
719 beta( shiftpos ) = eshift
720 beta( shiftpos+1 ) = eshift
721 ns = 2
722 END IF
723
724*
725* Time for a QZ sweep
726*
727 CALL dlaqz4( ilschur, ilq, ilz, n, istart2, istop, ns, nblock,
728 $ alphar( shiftpos ), alphai( shiftpos ),
729 $ beta( shiftpos ), a, lda, b, ldb, q, ldq, z, ldz,
730 $ work, nblock, work( nblock**2+1 ), nblock,
731 $ work( 2*nblock**2+1 ), lwork-2*nblock**2,
732 $ sweep_info )
733
734 END DO
735
736*
737* Call DHGEQZ to normalize the eigenvalue blocks and set the eigenvalues
738* If all the eigenvalues have been found, DHGEQZ will not do any iterations
739* and only normalize the blocks. In case of a rare convergence failure,
740* the single shift might perform better.
741*
742 80 CALL dhgeqz( wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb,
743 $ alphar, alphai, beta, q, ldq, z, ldz, work, lwork,
744 $ norm_info )
745
746 info = norm_info
747
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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
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
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
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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
Here is the call graph for this function:
Here is the caller graph for this function: