dlaqz3()

recursive subroutine dlaqz3	(	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,
		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( ldq, * ), intent(inout)	Q,
		integer, intent(in)	LDQ,
		double precision, dimension( ldz, * ), intent(inout)	Z,
		integer, intent(in)	LDZ,
		integer, intent(out)	NS,
		integer, intent(out)	ND,
		double precision, dimension( * ), intent(inout)	ALPHAR,
		double precision, dimension( * ), intent(inout)	ALPHAI,
		double precision, dimension( * ), intent(inout)	BETA,
		double precision, dimension( ldqc, * )	QC,
		integer, intent(in)	LDQC,
		double precision, dimension( ldzc, * )	ZC,
		integer, intent(in)	LDZC,
		double precision, dimension( * )	WORK,
		integer, intent(in)	LWORK,
		integer, intent(in)	REC,
		integer, intent(out)	INFO
	)

DLAQZ3

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

Purpose:

 DLAQZ3 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDQ, N)
[in]	LDQ	LDQ is INTEGER
[in,out]	Z	Z is DOUBLE PRECISION 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 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]	QC	QC is DOUBLE PRECISION array, dimension (LDQC, NW)
[in]	LDQC	LDQC is INTEGER
[in,out]	ZC	ZC is DOUBLE PRECISION array, dimension (LDZC, NW)
[in]	LDZC	LDZ is INTEGER
[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

Author

Thijs Steel, KU Leuven

Date

May 2020

Definition at line 235 of file dlaqz3.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
 
      DOUBLE PRECISION, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ),
     $                  Q( LDQ, * ), Z( LDZ, * ), ALPHAR( * ),
     $                  ALPHAI( * ), BETA( * )
      INTEGER, INTENT( OUT ) :: NS, ND, INFO
      DOUBLE PRECISION :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
 
*     Parameters
      DOUBLE PRECISION :: ZERO, ONE, HALF
      parameter( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
 
*     Local Scalars
      LOGICAL :: BULGE
      INTEGER :: JW, KWTOP, KWBOT, ISTOPM, ISTARTM, K, K2, DTGEXC_INFO,
     $           IFST, ILST, LWORKREQ, QZ_SMALL_INFO
      DOUBLE PRECISION :: S, SMLNUM, ULP, SAFMIN, SAFMAX, C1, S1, TEMP
 
*     External Functions
      EXTERNAL :: xerbla, dtgexc, dlabad, dlaqz0, dlacpy, dlaset,
     $            dlaqz2, drot, dlartg, dlag2, dgemm
      DOUBLE PRECISION, EXTERNAL :: DLAMCH
 
      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 dtgexc( .true., .true., jw, a, lda, b, ldb, qc, ldqc, zc,
     $             ldzc, ifst, ilst, work, -1, dtgexc_info )
      lworkreq = int( work( 1 ) )
      CALL dlaqz0( '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 ) = lworkreq
         RETURN
      ELSE IF ( lwork .LT. lworkreq ) THEN
         info = -26
      END IF
 
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DLAQZ3', -info )
         RETURN
      END IF
 
*     Get machine constants
      safmin = dlamch( 'SAFE MINIMUM' )
      safmax = one/safmin
      CALL dlabad( safmin, safmax )
      ulp = dlamch( 'PRECISION' )
      smlnum = safmin*( dble( 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 dlacpy( 'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
      CALL dlacpy( 'ALL', jw, jw, b( kwtop, kwtop ), ldb, work( jw**2+
     $             1 ), jw )
 
*     Transform window to real schur form
      CALL dlaset( 'FULL', jw, jw, zero, one, qc, ldqc )
      CALL dlaset( 'FULL', jw, jw, zero, one, zc, ldzc )
      CALL dlaqz0( '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 dlacpy( 'ALL', jw, jw, work, jw, a( kwtop, kwtop ), lda )
         CALL dlacpy( '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 dtgexc( .true., .true., jw, a( kwtop, kwtop ),
     $                         lda, b( kwtop, kwtop ), ldb, qc, ldqc,
     $                         zc, ldzc, ifst, ilst, work, lwork,
     $                         dtgexc_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 dtgexc( .true., .true., jw, a( kwtop, kwtop ),
     $                         lda, b( kwtop, kwtop ), ldb, qc, ldqc,
     $                         zc, ldzc, ifst, ilst, work, lwork,
     $                         dtgexc_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 dlag2( 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 dlartg( 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 drot( ihi-k2+1, a( k, k2 ), lda, a( k+1, k2 ), lda, c1,
     $                 s1 )
            CALL drot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1, k-1 ),
     $                 ldb, c1, s1 )
            CALL drot( 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 dlaqz2( .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 dlartg( b( k2+1, k2+1 ), b( k2+1, k2 ), c1, s1,
     $                         temp )
                  b( k2+1, k2+1 ) = temp
                  b( k2+1, k2 ) = zero
                  CALL drot( k2+2-istartm+1, a( istartm, k2+1 ), 1,
     $                       a( istartm, k2 ), 1, c1, s1 )
                  CALL drot( k2-istartm+1, b( istartm, k2+1 ), 1,
     $                       b( istartm, k2 ), 1, c1, s1 )
                  CALL drot( jw, zc( 1, k2+1-kwtop+1 ), 1, zc( 1,
     $                       k2-kwtop+1 ), 1, c1, s1 )
            
                  CALL dlartg( a( k2+1, k2 ), a( k2+2, k2 ), c1, s1,
     $                         temp )
                  a( k2+1, k2 ) = temp
                  a( k2+2, k2 ) = zero
                  CALL drot( istopm-k2, a( k2+1, k2+1 ), lda, a( k2+2,
     $                       k2+1 ), lda, c1, s1 )
                  CALL drot( istopm-k2, b( k2+1, k2+1 ), ldb, b( k2+2,
     $                       k2+1 ), ldb, c1, s1 )
                  CALL drot( jw, qc( 1, k2+1-kwtop+1 ), 1, qc( 1,
     $                       k2+2-kwtop+1 ), 1, c1, s1 )
 
               END DO
 
*              Remove the shift
               CALL dlartg( b( kwbot, kwbot ), b( kwbot, kwbot-1 ), c1,
     $                      s1, temp )
               b( kwbot, kwbot ) = temp
               b( kwbot, kwbot-1 ) = zero
               CALL drot( kwbot-istartm, b( istartm, kwbot ), 1,
     $                    b( istartm, kwbot-1 ), 1, c1, s1 )
               CALL drot( kwbot-istartm+1, a( istartm, kwbot ), 1,
     $                    a( istartm, kwbot-1 ), 1, c1, s1 )
               CALL drot( 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 dgemm( 'T', 'N', jw, istopm-ihi, jw, one, qc, ldqc,
     $               a( kwtop, ihi+1 ), lda, zero, work, jw )
         CALL dlacpy( 'ALL', jw, istopm-ihi, work, jw, a( kwtop,
     $                ihi+1 ), lda )
         CALL dgemm( 'T', 'N', jw, istopm-ihi, jw, one, qc, ldqc,
     $               b( kwtop, ihi+1 ), ldb, zero, work, jw )
         CALL dlacpy( 'ALL', jw, istopm-ihi, work, jw, b( kwtop,
     $                ihi+1 ), ldb )
      END IF
      IF ( ilq ) THEN
         CALL dgemm( 'N', 'N', n, jw, jw, one, q( 1, kwtop ), ldq, qc,
     $               ldqc, zero, work, n )
         CALL dlacpy( 'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
      END IF
 
      IF ( kwtop-1-istartm+1 > 0 ) THEN
         CALL dgemm( 'N', 'N', kwtop-istartm, jw, jw, one, a( istartm,
     $               kwtop ), lda, zc, ldzc, zero, work,
     $               kwtop-istartm )
         CALL dlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
     $                a( istartm, kwtop ), lda )
         CALL dgemm( 'N', 'N', kwtop-istartm, jw, jw, one, b( istartm,
     $               kwtop ), ldb, zc, ldzc, zero, work,
     $               kwtop-istartm )
         CALL dlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
     $                b( istartm, kwtop ), ldb )
      END IF
      IF ( ilz ) THEN
         CALL dgemm( 'N', 'N', n, jw, jw, one, z( 1, kwtop ), ldz, zc,
     $               ldzc, zero, work, n )
         CALL dlacpy( '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: