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