slaqz3.f

*> \brief \b SLAQZ3
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download SLAQZ3 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaqz3.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaqz3.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaqz3.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*      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 )
*      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( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> SLAQZ3 performs AED
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] ILSCHUR
*> \verbatim
*>          ILSCHUR is LOGICAL
*>              Determines whether or not to update the full Schur form
*> \endverbatim
*> \param[in] ILQ
*> \verbatim
*>          ILQ is LOGICAL
*>              Determines whether or not to update the matrix Q
*> \endverbatim
*>
*> \param[in] ILZ
*> \verbatim
*>          ILZ is LOGICAL
*>              Determines whether or not to update the matrix Z
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrices A, B, Q, and Z.  N >= 0.
*> \endverbatim
*>
*> \param[in] ILO
*> \verbatim
*>          ILO is INTEGER
*> \endverbatim
*>
*> \param[in] IHI
*> \verbatim
*>          IHI is INTEGER
*>          ILO and IHI mark the rows and columns of (A,B) which
*>          are to be normalized
*> \endverbatim
*>
*> \param[in] NW
*> \verbatim
*>          NW is INTEGER
*>          The desired size of the deflation window.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is REAL array, dimension (LDA, N)
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A.  LDA >= max( 1, N ).
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*>          B is REAL array, dimension (LDB, N)
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>          The leading dimension of the array B.  LDB >= max( 1, N ).
*> \endverbatim
*>
*> \param[in,out] Q
*> \verbatim
*>          Q is REAL array, dimension (LDQ, N)
*> \endverbatim
*>
*> \param[in] LDQ
*> \verbatim
*>          LDQ is INTEGER
*> \endverbatim
*>
*> \param[in,out] Z
*> \verbatim
*>          Z is REAL array, dimension (LDZ, N)
*> \endverbatim
*>
*> \param[in] LDZ
*> \verbatim
*>          LDZ is INTEGER
*> \endverbatim
*>
*> \param[out] NS
*> \verbatim
*>          NS is INTEGER
*>          The number of unconverged eigenvalues available to
*>          use as shifts.
*> \endverbatim
*>
*> \param[out] ND
*> \verbatim
*>          ND is INTEGER
*>          The number of converged eigenvalues found.
*> \endverbatim
*>
*> \param[out] ALPHAR
*> \verbatim
*>          ALPHAR is REAL array, dimension (N)
*>          The real parts of each scalar alpha defining an eigenvalue
*>          of GNEP.
*> \endverbatim
*>
*> \param[out] ALPHAI
*> \verbatim
*>          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).
*> \endverbatim
*>
*> \param[out] BETA
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in,out] QC
*> \verbatim
*>          QC is REAL array, dimension (LDQC, NW)
*> \endverbatim
*>
*> \param[in] LDQC
*> \verbatim
*>          LDQC is INTEGER
*> \endverbatim
*>
*> \param[in,out] ZC
*> \verbatim
*>          ZC is REAL array, dimension (LDZC, NW)
*> \endverbatim
*>
*> \param[in] LDZC
*> \verbatim
*>          LDZ is INTEGER
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is REAL array, dimension (MAX(1,LWORK))
*>          On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] REC
*> \verbatim
*>          REC is INTEGER
*>             REC indicates the current recursion level. Should be set
*>             to 0 on first call.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0: successful exit
*>          < 0: if INFO = -i, the i-th argument had an illegal value
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Thijs Steel, KU Leuven
*
*> \date May 2020
*
*> \ingroup laqz3
*>
*  =====================================================================
      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 )
      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, slabad, slaqz0, slacpy, slaset,
     $            slaqz2, srot, slartg, slag2, sgemm
      REAL, EXTERNAL :: SLAMCH

      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 ) = 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
      CALL slabad( safmin, safmax )
      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

      END SUBROUTINE