slaqz3()

recursive subroutine slaqz3	(	logical, intent(in)	ilschur,
		logical, intent(in)	ilq,
		logical, intent(in)	ilz,
		integer, intent(in)	n,
		integer, intent(in)	ilo,
		integer, intent(in)	ihi,
		integer, intent(in)	nw,
		real, dimension( lda, * ), intent(inout)	a,
		integer, intent(in)	lda,
		real, dimension( ldb, * ), intent(inout)	b,
		integer, intent(in)	ldb,
		real, dimension( ldq, * ), intent(inout)	q,
		integer, intent(in)	ldq,
		real, dimension( ldz, * ), intent(inout)	z,
		integer, intent(in)	ldz,
		integer, intent(out)	ns,
		integer, intent(out)	nd,
		real, dimension( * ), intent(inout)	alphar,
		real, dimension( * ), intent(inout)	alphai,
		real, dimension( * ), intent(inout)	beta,
		real, dimension( ldqc, * )	qc,
		integer, intent(in)	ldqc,
		real, dimension( ldzc, * )	zc,
		integer, intent(in)	ldzc,
		real, dimension( * )	work,
		integer, intent(in)	lwork,
		integer, intent(in)	rec,
		integer, intent(out)	info
	)

SLAQZ3

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

Purpose:

 SLAQZ3 performs AED

Parameters

[in]	ILSCHUR	ILSCHUR is LOGICAL Determines whether or not to update the full Schur form
[in]	ILQ	ILQ is LOGICAL Determines whether or not to update the matrix Q
[in]	ILZ	ILZ is LOGICAL Determines whether or not to update the matrix Z
[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,B) which are to be normalized
[in]	NW	NW is INTEGER The desired size of the deflation window.
[in,out]	A	A is REAL array, dimension (LDA, N)
[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)
[in]	LDB	LDB is INTEGER The leading dimension of the array B. LDB >= max( 1, N ).
[in,out]	Q	Q is REAL array, dimension (LDQ, N)
[in]	LDQ	LDQ is INTEGER
[in,out]	Z	Z is REAL array, dimension (LDZ, N)
[in]	LDZ	LDZ is INTEGER
[out]	NS	NS is INTEGER The number of unconverged eigenvalues available to use as shifts.
[out]	ND	ND is INTEGER The number of converged eigenvalues found.
[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]	QC	QC is REAL array, dimension (LDQC, NW)
[in]	LDQC	LDQC is INTEGER
[in,out]	ZC	ZC is REAL array, dimension (LDZC, NW)
[in]	LDZC	LDZ is INTEGER
[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

Author

Thijs Steel, KU Leuven

Date

May 2020

Definition at line 234 of file slaqz3.f.

      IMPLICIT NONE
 
*     Arguments
      LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
      INTEGER, INTENT( IN ) :: N, ILO, IHI, NW, LDA, LDB, LDQ, LDZ,
     $         LDQC, LDZC, LWORK, REC
 
      REAL, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
     $   Z( LDZ, * ), ALPHAR( * ), ALPHAI( * ), BETA( * )
      INTEGER, INTENT( OUT ) :: NS, ND, INFO
      REAL :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
 
*     Parameters
      REAL :: ZERO, ONE, HALF
      parameter( zero = 0.0, one = 1.0, half = 0.5 )
 
*     Local Scalars
      LOGICAL :: BULGE
      INTEGER :: JW, KWTOP, KWBOT, ISTOPM, ISTARTM, K, K2, STGEXC_INFO,
     $           IFST, ILST, LWORKREQ, QZ_SMALL_INFO
      REAL :: S, SMLNUM, ULP, SAFMIN, SAFMAX, C1, S1, TEMP
 
*     External Functions
      EXTERNAL :: xerbla, stgexc, slaqz0, slacpy, slaset,
     $            slaqz2, srot, slartg, slag2, sgemm
      REAL, EXTERNAL :: SLAMCH, SROUNDUP_LWORK
 
      info = 0
 
*     Set up deflation window
      jw = min( nw, ihi-ilo+1 )
      kwtop = ihi-jw+1
      IF ( kwtop .EQ. ilo ) THEN
         s = zero
      ELSE
         s = a( kwtop, kwtop-1 )
      END IF
 
*     Determine required workspace
      ifst = 1
      ilst = jw
      CALL stgexc( .true., .true., jw, a, lda, b, ldb, qc, ldqc, zc,
     $             ldzc, ifst, ilst, work, -1, stgexc_info )
      lworkreq = int( work( 1 ) )
      CALL slaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
     $             b( kwtop, kwtop ), ldb, alphar, alphai, beta, qc,
     $             ldqc, zc, ldzc, work, -1, rec+1, qz_small_info )
      lworkreq = max( lworkreq, int( work( 1 ) )+2*jw**2 )
      lworkreq = max( lworkreq, n*nw, 2*nw**2+n )
      IF ( lwork .EQ.-1 ) THEN
*        workspace query, quick return
         work( 1 ) = sroundup_lwork(lworkreq)
         RETURN
      ELSE IF ( lwork .LT. lworkreq ) THEN
         info = -26
      END IF
 
      IF( info.NE.0 ) THEN
         CALL xerbla( 'SLAQZ3', -info )
         RETURN
      END IF
 
*     Get machine constants
      safmin = slamch( 'SAFE MINIMUM' )
      safmax = one/safmin
      ulp = slamch( 'PRECISION' )
      smlnum = safmin*( real( n )/ulp )
 
      IF ( ihi .EQ. kwtop ) THEN
*        1 by 1 deflation window, just try a regular deflation
         alphar( kwtop ) = a( kwtop, kwtop )
         alphai( kwtop ) = zero
         beta( kwtop ) = b( kwtop, kwtop )
         ns = 1
         nd = 0
         IF ( abs( s ) .LE. max( smlnum, ulp*abs( a( kwtop,
     $      kwtop ) ) ) ) THEN
            ns = 0
            nd = 1
            IF ( kwtop .GT. ilo ) THEN
               a( kwtop, kwtop-1 ) = zero
            END IF
         END IF
      END IF
 
 
*     Store window in case of convergence failure
      CALL slacpy( 'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
      CALL slacpy( 'ALL', jw, jw, b( kwtop, kwtop ), ldb, work( jw**2+
     $             1 ), jw )
 
*     Transform window to real schur form
      CALL slaset( 'FULL', jw, jw, zero, one, qc, ldqc )
      CALL slaset( 'FULL', jw, jw, zero, one, zc, ldzc )
      CALL slaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
     $             b( kwtop, kwtop ), ldb, alphar, alphai, beta, qc,
     $             ldqc, zc, ldzc, work( 2*jw**2+1 ), lwork-2*jw**2,
     $             rec+1, qz_small_info )
 
      IF( qz_small_info .NE. 0 ) THEN
*        Convergence failure, restore the window and exit
         nd = 0
         ns = jw-qz_small_info
         CALL slacpy( 'ALL', jw, jw, work, jw, a( kwtop, kwtop ), lda )
         CALL slacpy( 'ALL', jw, jw, work( jw**2+1 ), jw, b( kwtop,
     $                kwtop ), ldb )
         RETURN
      END IF
 
*     Deflation detection loop
      IF ( kwtop .EQ. ilo .OR. s .EQ. zero ) THEN
         kwbot = kwtop-1
      ELSE
         kwbot = ihi
         k = 1
         k2 = 1
         DO WHILE ( k .LE. jw )
            bulge = .false.
            IF ( kwbot-kwtop+1 .GE. 2 ) THEN
               bulge = a( kwbot, kwbot-1 ) .NE. zero
            END IF
            IF ( bulge ) THEN
 
*              Try to deflate complex conjugate eigenvalue pair
               temp = abs( a( kwbot, kwbot ) )+sqrt( abs( a( kwbot,
     $            kwbot-1 ) ) )*sqrt( abs( a( kwbot-1, kwbot ) ) )
               IF( temp .EQ. zero )THEN
                  temp = abs( s )
               END IF
               IF ( max( abs( s*qc( 1, kwbot-kwtop ) ), abs( s*qc( 1,
     $            kwbot-kwtop+1 ) ) ) .LE. max( smlnum,
     $            ulp*temp ) ) THEN
*                 Deflatable
                  kwbot = kwbot-2
               ELSE
*                 Not deflatable, move out of the way
                  ifst = kwbot-kwtop+1
                  ilst = k2
                  CALL stgexc( .true., .true., jw, a( kwtop, kwtop ),
     $                         lda, b( kwtop, kwtop ), ldb, qc, ldqc,
     $                         zc, ldzc, ifst, ilst, work, lwork,
     $                         stgexc_info )
                  k2 = k2+2
               END IF
               k = k+2
            ELSE
 
*              Try to deflate real eigenvalue
               temp = abs( a( kwbot, kwbot ) )
               IF( temp .EQ. zero ) THEN
                  temp = abs( s )
               END IF
               IF ( ( abs( s*qc( 1, kwbot-kwtop+1 ) ) ) .LE. max( ulp*
     $            temp, smlnum ) ) THEN
*                 Deflatable
                  kwbot = kwbot-1
               ELSE
*                 Not deflatable, move out of the way
                  ifst = kwbot-kwtop+1
                  ilst = k2
                  CALL stgexc( .true., .true., jw, a( kwtop, kwtop ),
     $                         lda, b( kwtop, kwtop ), ldb, qc, ldqc,
     $                         zc, ldzc, ifst, ilst, work, lwork,
     $                         stgexc_info )
                  k2 = k2+1
               END IF
 
               k = k+1
 
            END IF
         END DO
      END IF
 
*     Store eigenvalues
      nd = ihi-kwbot
      ns = jw-nd
      k = kwtop
      DO WHILE ( k .LE. ihi )
         bulge = .false.
         IF ( k .LT. ihi ) THEN
            IF ( a( k+1, k ) .NE. zero ) THEN
               bulge = .true.
            END IF
         END IF
         IF ( bulge ) THEN
*           2x2 eigenvalue block
            CALL slag2( a( k, k ), lda, b( k, k ), ldb, safmin,
     $                  beta( k ), beta( k+1 ), alphar( k ),
     $                  alphar( k+1 ), alphai( k ) )
            alphai( k+1 ) = -alphai( k )
            k = k+2
         ELSE
*           1x1 eigenvalue block
            alphar( k ) = a( k, k )
            alphai( k ) = zero
            beta( k ) = b( k, k )
            k = k+1
         END IF
      END DO
 
      IF ( kwtop .NE. ilo .AND. s .NE. zero ) THEN
*        Reflect spike back, this will create optimally packed bulges
         a( kwtop:kwbot, kwtop-1 ) = a( kwtop, kwtop-1 )*qc( 1,
     $      1:jw-nd )
         DO k = kwbot-1, kwtop, -1
            CALL slartg( a( k, kwtop-1 ), a( k+1, kwtop-1 ), c1, s1,
     $                   temp )
            a( k, kwtop-1 ) = temp
            a( k+1, kwtop-1 ) = zero
            k2 = max( kwtop, k-1 )
            CALL srot( ihi-k2+1, a( k, k2 ), lda, a( k+1, k2 ), lda, c1,
     $                 s1 )
            CALL srot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1, k-1 ),
     $                 ldb, c1, s1 )
            CALL srot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1, k+1-kwtop+1 ),
     $                 1, c1, s1 )
         END DO
 
*        Chase bulges down
         istartm = kwtop
         istopm = ihi
         k = kwbot-1
         DO WHILE ( k .GE. kwtop )
            IF ( ( k .GE. kwtop+1 ) .AND. a( k+1, k-1 ) .NE. zero ) THEN
 
*              Move double pole block down and remove it
               DO k2 = k-1, kwbot-2
                  CALL slaqz2( .true., .true., k2, kwtop, kwtop+jw-1,
     $                         kwbot, a, lda, b, ldb, jw, kwtop, qc,
     $                         ldqc, jw, kwtop, zc, ldzc )
               END DO
 
               k = k-2
            ELSE
 
*              k points to single shift
               DO k2 = k, kwbot-2
 
*                 Move shift down
                  CALL slartg( b( k2+1, k2+1 ), b( k2+1, k2 ), c1, s1,
     $                         temp )
                  b( k2+1, k2+1 ) = temp
                  b( k2+1, k2 ) = zero
                  CALL srot( k2+2-istartm+1, a( istartm, k2+1 ), 1,
     $                       a( istartm, k2 ), 1, c1, s1 )
                  CALL srot( k2-istartm+1, b( istartm, k2+1 ), 1,
     $                       b( istartm, k2 ), 1, c1, s1 )
                  CALL srot( jw, zc( 1, k2+1-kwtop+1 ), 1, zc( 1,
     $                       k2-kwtop+1 ), 1, c1, s1 )
            
                  CALL slartg( a( k2+1, k2 ), a( k2+2, k2 ), c1, s1,
     $                         temp )
                  a( k2+1, k2 ) = temp
                  a( k2+2, k2 ) = zero
                  CALL srot( istopm-k2, a( k2+1, k2+1 ), lda, a( k2+2,
     $                       k2+1 ), lda, c1, s1 )
                  CALL srot( istopm-k2, b( k2+1, k2+1 ), ldb, b( k2+2,
     $                       k2+1 ), ldb, c1, s1 )
                  CALL srot( jw, qc( 1, k2+1-kwtop+1 ), 1, qc( 1,
     $                       k2+2-kwtop+1 ), 1, c1, s1 )
 
               END DO
 
*              Remove the shift
               CALL slartg( b( kwbot, kwbot ), b( kwbot, kwbot-1 ), c1,
     $                      s1, temp )
               b( kwbot, kwbot ) = temp
               b( kwbot, kwbot-1 ) = zero
               CALL srot( kwbot-istartm, b( istartm, kwbot ), 1,
     $                    b( istartm, kwbot-1 ), 1, c1, s1 )
               CALL srot( kwbot-istartm+1, a( istartm, kwbot ), 1,
     $                    a( istartm, kwbot-1 ), 1, c1, s1 )
               CALL srot( jw, zc( 1, kwbot-kwtop+1 ), 1, zc( 1,
     $                    kwbot-1-kwtop+1 ), 1, c1, s1 )
 
               k = k-1
            END IF
         END DO
 
      END IF
 
*     Apply Qc and Zc to rest of the matrix
      IF ( ilschur ) THEN
         istartm = 1
         istopm = n
      ELSE
         istartm = ilo
         istopm = ihi
      END IF
 
      IF ( istopm-ihi > 0 ) THEN
         CALL sgemm( 'T', 'N', jw, istopm-ihi, jw, one, qc, ldqc,
     $               a( kwtop, ihi+1 ), lda, zero, work, jw )
         CALL slacpy( 'ALL', jw, istopm-ihi, work, jw, a( kwtop,
     $                ihi+1 ), lda )
         CALL sgemm( 'T', 'N', jw, istopm-ihi, jw, one, qc, ldqc,
     $               b( kwtop, ihi+1 ), ldb, zero, work, jw )
         CALL slacpy( 'ALL', jw, istopm-ihi, work, jw, b( kwtop,
     $                ihi+1 ), ldb )
      END IF
      IF ( ilq ) THEN
         CALL sgemm( 'N', 'N', n, jw, jw, one, q( 1, kwtop ), ldq, qc,
     $               ldqc, zero, work, n )
         CALL slacpy( 'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
      END IF
 
      IF ( kwtop-1-istartm+1 > 0 ) THEN
         CALL sgemm( 'N', 'N', kwtop-istartm, jw, jw, one, a( istartm,
     $               kwtop ), lda, zc, ldzc, zero, work,
     $               kwtop-istartm )
         CALL slacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
     $                a( istartm, kwtop ), lda )
         CALL sgemm( 'N', 'N', kwtop-istartm, jw, jw, one, b( istartm,
     $               kwtop ), ldb, zc, ldzc, zero, work,
     $               kwtop-istartm )
         CALL slacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
     $                b( istartm, kwtop ), ldb )
      END IF
      IF ( ilz ) THEN
         CALL sgemm( 'N', 'N', n, jw, jw, one, z( 1, kwtop ), ldz, zc,
     $               ldzc, zero, work, n )
         CALL slacpy( 'ALL', n, jw, work, n, z( 1, kwtop ), ldz )
      END IF
 

Here is the call graph for this function:

Here is the caller graph for this function: