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 233 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, 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
      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: