LAPACK 3.12.1
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, , SIAM J. Numer. Anal., 10(1973),
!>      pp. 241--256.
!>
!> Ref: B. Kagstrom, D. Kressner, , SIAM J. Numer.
!>      Anal., 29(2006), pp. 199--227.
!>
!> Ref: T. Steel, D. Camps, K. Meerbergen, R. Vandebril 
!>
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 298 of file slaqz0.f.

303 IMPLICIT NONE
304
305* Arguments
306 CHARACTER, INTENT( IN ) :: WANTS, WANTQ, WANTZ
307 INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
308 $ REC
309
310 INTEGER, INTENT( OUT ) :: INFO
311
312 REAL, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
313 $ Z( LDZ, * ), ALPHAR( * ), ALPHAI( * ), BETA( * ), WORK( * )
314
315* Parameters
316 REAL :: ZERO, ONE, HALF
317 parameter( zero = 0.0, one = 1.0, half = 0.5 )
318
319* Local scalars
320 REAL :: SMLNUM, ULP, ESHIFT, SAFMIN, SAFMAX, C1, S1, TEMP, SWAP,
321 $ BNORM, BTOL
322 INTEGER :: ISTART, ISTOP, IITER, MAXIT, ISTART2, K, LD, NSHIFTS,
323 $ NBLOCK, NW, NMIN, NIBBLE, N_UNDEFLATED, N_DEFLATED,
324 $ NS, SWEEP_INFO, SHIFTPOS, LWORKREQ, K2, ISTARTM,
325 $ ISTOPM, IWANTS, IWANTQ, IWANTZ, NORM_INFO, AED_INFO,
326 $ NWR, NBR, NSR, ITEMP1, ITEMP2, RCOST, I
327 LOGICAL :: ILSCHUR, ILQ, ILZ
328 CHARACTER :: JBCMPZ*3
329
330* External Functions
331 EXTERNAL :: xerbla, shgeqz, slaqz3, slaqz4, slaset,
332 $ slartg, srot
333 REAL, EXTERNAL :: SLAMCH, SLANHS, SROUNDUP_LWORK
334 LOGICAL, EXTERNAL :: LSAME
335 INTEGER, EXTERNAL :: ILAENV
336
337*
338* Decode wantS,wantQ,wantZ
339*
340 IF( lsame( wants, 'E' ) ) THEN
341 ilschur = .false.
342 iwants = 1
343 ELSE IF( lsame( wants, 'S' ) ) THEN
344 ilschur = .true.
345 iwants = 2
346 ELSE
347 iwants = 0
348 END IF
349
350 IF( lsame( wantq, 'N' ) ) THEN
351 ilq = .false.
352 iwantq = 1
353 ELSE IF( lsame( wantq, 'V' ) ) THEN
354 ilq = .true.
355 iwantq = 2
356 ELSE IF( lsame( wantq, 'I' ) ) THEN
357 ilq = .true.
358 iwantq = 3
359 ELSE
360 iwantq = 0
361 END IF
362
363 IF( lsame( wantz, 'N' ) ) THEN
364 ilz = .false.
365 iwantz = 1
366 ELSE IF( lsame( wantz, 'V' ) ) THEN
367 ilz = .true.
368 iwantz = 2
369 ELSE IF( lsame( wantz, 'I' ) ) THEN
370 ilz = .true.
371 iwantz = 3
372 ELSE
373 iwantz = 0
374 END IF
375*
376* Check Argument Values
377*
378 info = 0
379 IF( iwants.EQ.0 ) THEN
380 info = -1
381 ELSE IF( iwantq.EQ.0 ) THEN
382 info = -2
383 ELSE IF( iwantz.EQ.0 ) THEN
384 info = -3
385 ELSE IF( n.LT.0 ) THEN
386 info = -4
387 ELSE IF( ilo.LT.1 ) THEN
388 info = -5
389 ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 ) THEN
390 info = -6
391 ELSE IF( lda.LT.n ) THEN
392 info = -8
393 ELSE IF( ldb.LT.n ) THEN
394 info = -10
395 ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) ) THEN
396 info = -15
397 ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) ) THEN
398 info = -17
399 END IF
400 IF( info.NE.0 ) THEN
401 CALL xerbla( 'SLAQZ0', -info )
402 RETURN
403 END IF
404
405*
406* Quick return if possible
407*
408 IF( n.LE.0 ) THEN
409 work( 1 ) = real( 1 )
410 RETURN
411 END IF
412
413*
414* Get the parameters
415*
416 jbcmpz( 1:1 ) = wants
417 jbcmpz( 2:2 ) = wantq
418 jbcmpz( 3:3 ) = wantz
419
420 nmin = ilaenv( 12, 'SLAQZ0', jbcmpz, n, ilo, ihi, lwork )
421
422 nwr = ilaenv( 13, 'SLAQZ0', jbcmpz, n, ilo, ihi, lwork )
423 nwr = max( 2, nwr )
424 nwr = min( ihi-ilo+1, ( n-1 ) / 3, nwr )
425
426 nibble = ilaenv( 14, 'SLAQZ0', jbcmpz, n, ilo, ihi, lwork )
427
428 nsr = ilaenv( 15, 'SLAQZ0', jbcmpz, n, ilo, ihi, lwork )
429 nsr = min( nsr, ( n+6 ) / 9, ihi-ilo )
430 nsr = max( 2, nsr-mod( nsr, 2 ) )
431
432 rcost = ilaenv( 17, 'SLAQZ0', jbcmpz, n, ilo, ihi, lwork )
433 itemp1 = int( real( nsr )/sqrt( 1+2*real( nsr )/
434 $ ( real( rcost )/100*real( 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,
440 $ ldb,
441 $ alphar, alphai, beta, q, ldq, z, ldz, work,
442 $ lwork, info )
443 RETURN
444 END IF
445
446*
447* Find out required workspace
448*
449
450* Workspace query to slaqz3
451 nw = max( nwr, nmin )
452 CALL slaqz3( ilschur, ilq, ilz, n, ilo, ihi, nw, a, lda, b,
453 $ ldb,
454 $ q, ldq, z, ldz, n_undeflated, n_deflated, alphar,
455 $ alphai, beta, work, nw, work, nw, work, -1, rec,
456 $ aed_info )
457 itemp1 = int( work( 1 ) )
458* Workspace query to slaqz4
459 CALL slaqz4( ilschur, ilq, ilz, n, ilo, ihi, nsr, nbr, alphar,
460 $ alphai, beta, a, lda, b, ldb, q, ldq, z, ldz, work,
461 $ nbr, work, nbr, work, -1, sweep_info )
462 itemp2 = int( work( 1 ) )
463
464 lworkreq = max( itemp1+2*nw**2, itemp2+2*nbr**2 )
465 IF ( lwork .EQ.-1 ) THEN
466 work( 1 ) = sroundup_lwork( lworkreq )
467 RETURN
468 ELSE IF ( lwork .LT. lworkreq ) THEN
469 info = -19
470 END IF
471 IF( info.NE.0 ) THEN
472 CALL xerbla( 'SLAQZ0', info )
473 RETURN
474 END IF
475*
476* Initialize Q and Z
477*
478 IF( iwantq.EQ.3 ) CALL slaset( 'FULL', n, n, zero, one, q,
479 $ ldq )
480 IF( iwantz.EQ.3 ) CALL slaset( 'FULL', n, n, zero, one, z,
481 $ ldz )
482
483* Get machine constants
484 safmin = slamch( 'SAFE MINIMUM' )
485 safmax = one/safmin
486 ulp = slamch( 'PRECISION' )
487 smlnum = safmin*( real( n )/ulp )
488
489 bnorm = slanhs( '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 negligible, move it
571* to the top and deflate it
572
573 DO k2 = k, istart2+1, -1
574 CALL slartg( b( k2-1, k2 ), b( k2-1, k2-1 ), c1,
575 $ s1,
576 $ temp )
577 b( k2-1, k2 ) = temp
578 b( k2-1, k2-1 ) = zero
579
580 CALL srot( k2-2-istartm+1, b( istartm, k2 ), 1,
581 $ b( istartm, k2-1 ), 1, c1, s1 )
582 CALL srot( min( k2+1, istop )-istartm+1,
583 $ a( istartm,
584 $ k2 ), 1, a( istartm, k2-1 ), 1, c1, s1 )
585 IF ( ilz ) THEN
586 CALL srot( n, z( 1, k2 ), 1, z( 1, k2-1 ), 1,
587 $ c1,
588 $ s1 )
589 END IF
590
591 IF( k2.LT.istop ) THEN
592 CALL slartg( a( k2, k2-1 ), a( k2+1, k2-1 ), c1,
593 $ s1, temp )
594 a( k2, k2-1 ) = temp
595 a( k2+1, k2-1 ) = zero
596
597 CALL srot( istopm-k2+1, a( k2, k2 ), lda,
598 $ a( k2+1,
599 $ k2 ), lda, c1, s1 )
600 CALL srot( istopm-k2+1, b( k2, k2 ), ldb,
601 $ b( k2+1,
602 $ k2 ), ldb, c1, s1 )
603 IF( ilq ) THEN
604 CALL srot( n, q( 1, k2 ), 1, q( 1, k2+1 ), 1,
605 $ c1, s1 )
606 END IF
607 END IF
608
609 END DO
610
611 IF( istart2.LT.istop )THEN
612 CALL slartg( a( istart2, istart2 ), a( istart2+1,
613 $ istart2 ), c1, s1, temp )
614 a( istart2, istart2 ) = temp
615 a( istart2+1, istart2 ) = zero
616
617 CALL srot( istopm-( istart2+1 )+1, a( istart2,
618 $ istart2+1 ), lda, a( istart2+1,
619 $ istart2+1 ), lda, c1, s1 )
620 CALL srot( istopm-( istart2+1 )+1, b( istart2,
621 $ istart2+1 ), ldb, b( istart2+1,
622 $ istart2+1 ), ldb, c1, s1 )
623 IF( ilq ) THEN
624 CALL srot( n, q( 1, istart2 ), 1, q( 1,
625 $ istart2+1 ), 1, c1, s1 )
626 END IF
627 END IF
628
629 istart2 = istart2+1
630
631 END IF
632 k = k-1
633 END DO
634
635* istart2 now points to the top of the bottom right
636* unreduced Hessenberg block
637 IF ( istart2 .GE. istop ) THEN
638 istop = istart2-1
639 ld = 0
640 eshift = zero
641 cycle
642 END IF
643
644 nw = nwr
645 nshifts = nsr
646 nblock = nbr
647
648 IF ( istop-istart2+1 .LT. nmin ) THEN
649* Setting nw to the size of the subblock will make AED deflate
650* all the eigenvalues. This is slightly more efficient than just
651* using qz_small because the off diagonal part gets updated via BLAS.
652 IF ( istop-istart+1 .LT. nmin ) THEN
653 nw = istop-istart+1
654 istart2 = istart
655 ELSE
656 nw = istop-istart2+1
657 END IF
658 END IF
659
660*
661* Time for AED
662*
663 CALL slaqz3( ilschur, ilq, ilz, n, istart2, istop, nw, a,
664 $ lda,
665 $ b, ldb, q, ldq, z, ldz, n_undeflated, n_deflated,
666 $ alphar, alphai, beta, work, nw, work( nw**2+1 ),
667 $ nw, work( 2*nw**2+1 ), lwork-2*nw**2, rec,
668 $ aed_info )
669
670 IF ( n_deflated > 0 ) THEN
671 istop = istop-n_deflated
672 ld = 0
673 eshift = zero
674 END IF
675
676 IF ( 100*n_deflated > nibble*( n_deflated+n_undeflated ) .OR.
677 $ istop-istart2+1 .LT. nmin ) THEN
678* AED has uncovered many eigenvalues. Skip a QZ sweep and run
679* AED again.
680 cycle
681 END IF
682
683 ld = ld+1
684
685 ns = min( nshifts, istop-istart2 )
686 ns = min( ns, n_undeflated )
687 shiftpos = istop-n_undeflated+1
688*
689* Shuffle shifts to put double shifts in front
690* This ensures that we don't split up a double shift
691*
692 DO i = shiftpos, shiftpos+n_undeflated-1, 2
693 IF( alphai( i ).NE.-alphai( i+1 ) ) THEN
694*
695 swap = alphar( i )
696 alphar( i ) = alphar( i+1 )
697 alphar( i+1 ) = alphar( i+2 )
698 alphar( i+2 ) = swap
699
700 swap = alphai( i )
701 alphai( i ) = alphai( i+1 )
702 alphai( i+1 ) = alphai( i+2 )
703 alphai( i+2 ) = swap
704
705 swap = beta( i )
706 beta( i ) = beta( i+1 )
707 beta( i+1 ) = beta( i+2 )
708 beta( i+2 ) = swap
709 END IF
710 END DO
711
712 IF ( mod( ld, 6 ) .EQ. 0 ) THEN
713*
714* Exceptional shift. Chosen for no particularly good reason.
715*
716 IF( ( real( maxit )*safmin )*abs( a( istop,
717 $ istop-1 ) ).LT.abs( a( istop-1, istop-1 ) ) ) THEN
718 eshift = a( istop, istop-1 )/b( istop-1, istop-1 )
719 ELSE
720 eshift = eshift+one/( safmin*real( maxit ) )
721 END IF
722 alphar( shiftpos ) = one
723 alphar( shiftpos+1 ) = zero
724 alphai( shiftpos ) = zero
725 alphai( shiftpos+1 ) = zero
726 beta( shiftpos ) = eshift
727 beta( shiftpos+1 ) = eshift
728 ns = 2
729 END IF
730
731*
732* Time for a QZ sweep
733*
734 CALL slaqz4( ilschur, ilq, ilz, n, istart2, istop, ns,
735 $ nblock,
736 $ alphar( shiftpos ), alphai( shiftpos ),
737 $ beta( shiftpos ), a, lda, b, ldb, q, ldq, z, ldz,
738 $ work, nblock, work( nblock**2+1 ), nblock,
739 $ work( 2*nblock**2+1 ), lwork-2*nblock**2,
740 $ sweep_info )
741
742 END DO
743
744*
745* Call SHGEQZ to normalize the eigenvalue blocks and set the eigenvalues
746* If all the eigenvalues have been found, SHGEQZ will not do any iterations
747* and only normalize the blocks. In case of a rare convergence failure,
748* the single shift might perform better.
749*
750 80 CALL shgeqz( wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb,
751 $ alphar, alphai, beta, q, ldq, z, ldz, work, lwork,
752 $ norm_info )
753
754 info = norm_info
755
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
shgeqz
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:303
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
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
slaqz3
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:237
slaqz4
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:212
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
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: