zlaqz2.f

*> \brief \b ZLAQZ2
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download ZLAQZ2 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ZLAQZ2.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ZLAQZ2.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ZLAQZ2.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*      SUBROUTINE ZLAQZ2( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW, A, LDA, B,
*     $    LDB, Q, LDQ, Z, LDZ, NS, ND, ALPHA, BETA, QC, LDQC, ZC, LDZC,
*     $    WORK, LWORK, RWORK, 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
*
*      COMPLEX*16, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ,
*     $    * ), Z( LDZ, * ), ALPHA( * ), BETA( * )
*      INTEGER, INTENT( OUT ) :: NS, ND, INFO
*      COMPLEX*16 :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
*      DOUBLE PRECISION :: RWORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> ZLAQZ2 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 array, dimension (LDQ, N)
*> \endverbatim
*>
*> \param[in] LDQ
*> \verbatim
*>          LDQ is INTEGER
*> \endverbatim
*>
*> \param[in,out] Z
*> \verbatim
*>          Z is COMPLEX*16 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] ALPHA
*> \verbatim
*>          ALPHA is COMPLEX*16 array, dimension (N)
*>          Each scalar alpha defining an eigenvalue
*>          of GNEP.
*> \endverbatim
*>
*> \param[out] BETA
*> \verbatim
*>          BETA is COMPLEX*16 array, dimension (N)
*>          The scalars beta that define the eigenvalues of GNEP.
*>          Together, the quantities alpha = ALPHA(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 COMPLEX*16 array, dimension (LDQC, NW)
*> \endverbatim
*>
*> \param[in] LDQC
*> \verbatim
*>          LDQC is INTEGER
*> \endverbatim
*>
*> \param[in,out] ZC
*> \verbatim
*>          ZC is COMPLEX*16 array, dimension (LDZC, NW)
*> \endverbatim
*>
*> \param[in] LDZC
*> \verbatim
*>          LDZ is INTEGER
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is COMPLEX*16 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[out] RWORK
*> \verbatim
*>          RWORK is DOUBLE PRECISION array, dimension (N)
*> \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 laqz2
*>
*  =====================================================================
      RECURSIVE SUBROUTINE zlaqz2( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW,
     $                             A, LDA, B, LDB, Q, LDQ, Z, LDZ, NS,
     $                             ND, ALPHA, BETA, QC, LDQC, ZC, LDZC,
     $                             WORK, LWORK, RWORK, 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

      COMPLEX*16, INTENT( INOUT ) :: A( lda, * ), B( ldb, * ), Q( ldq,
     $   * ), z( ldz, * ), alpha( * ), beta( * )
      INTEGER, INTENT( OUT ) :: NS, ND, INFO
      COMPLEX*16 :: QC( ldqc, * ), ZC( ldzc, * ), WORK( * )
      DOUBLE PRECISION :: RWORK( * )

*     Parameters
      COMPLEX*16         CZERO, CONE
      parameter( czero = ( 0.0d+0, 0.0d+0 ), cone = ( 1.0d+0,
     $                     0.0d+0 ) )
      DOUBLE PRECISION :: ZERO, ONE, HALF
      parameter( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )

*     Local Scalars
      INTEGER :: JW, KWTOP, KWBOT, ISTOPM, ISTARTM, K, K2, ZTGEXC_INFO,
     $           ifst, ilst, lworkreq, qz_small_info
      DOUBLE PRECISION ::SMLNUM, ULP, SAFMIN, SAFMAX, C1, TEMPR
      COMPLEX*16 :: S, S1, TEMP

*     External Functions
      EXTERNAL :: xerbla, zlaqz0, zlaqz1, dlabad, zlacpy, zlaset, zgemm,
     $            ztgexc, zlartg, zrot
      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 = czero
      ELSE
         s = a( kwtop, kwtop-1 )
      END IF

*     Determine required workspace
      ifst = 1
      ilst = jw
      CALL zlaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
     $             b( kwtop, kwtop ), ldb, alpha, beta, qc, ldqc, zc,
     $             ldzc, work, -1, rwork, rec+1, qz_small_info )
      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( 'ZLAQZ2', -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
         alpha( kwtop ) = a( kwtop, kwtop )
         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 ) = czero
            END IF
         END IF
      END IF


*     Store window in case of convergence failure
      CALL zlacpy( 'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
      CALL zlacpy( 'ALL', jw, jw, b( kwtop, kwtop ), ldb, work( jw**2+
     $             1 ), jw )

*     Transform window to real schur form
      CALL zlaset( 'FULL', jw, jw, czero, cone, qc, ldqc )
      CALL zlaset( 'FULL', jw, jw, czero, cone, zc, ldzc )
      CALL zlaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
     $             b( kwtop, kwtop ), ldb, alpha, beta, qc, ldqc, zc,
     $             ldzc, work( 2*jw**2+1 ), lwork-2*jw**2, rwork,
     $             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 zlacpy( 'ALL', jw, jw, work, jw, a( kwtop, kwtop ), lda )
         CALL zlacpy( '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. czero ) THEN
         kwbot = kwtop-1
      ELSE
         kwbot = ihi
         k = 1
         k2 = 1
         DO WHILE ( k .LE. jw )
*              Try to deflate eigenvalue
               tempr = abs( a( kwbot, kwbot ) )
               IF( tempr .EQ. zero ) THEN
                  tempr = abs( s )
               END IF
               IF ( ( abs( s*qc( 1, kwbot-kwtop+1 ) ) ) .LE. max( ulp*
     $            tempr, smlnum ) ) THEN
*                 Deflatable
                  kwbot = kwbot-1
               ELSE
*                 Not deflatable, move out of the way
                  ifst = kwbot-kwtop+1
                  ilst = k2
                  CALL ztgexc( .true., .true., jw, a( kwtop, kwtop ),
     $                         lda, b( kwtop, kwtop ), ldb, qc, ldqc,
     $                         zc, ldzc, ifst, ilst, ztgexc_info )
                  k2 = k2+1
               END IF

               k = k+1
         END DO
      END IF

*     Store eigenvalues
      nd = ihi-kwbot
      ns = jw-nd
      k = kwtop
      DO WHILE ( k .LE. ihi )
         alpha( k ) = a( k, k )
         beta( k ) = b( k, k )
         k = k+1
      END DO

      IF ( kwtop .NE. ilo .AND. s .NE. czero ) THEN
*        Reflect spike back, this will create optimally packed bulges
         a( kwtop:kwbot, kwtop-1 ) = a( kwtop, kwtop-1 ) *dconjg( qc( 1,
     $      1:jw-nd ) )
         DO k = kwbot-1, kwtop, -1
            CALL zlartg( a( k, kwtop-1 ), a( k+1, kwtop-1 ), c1, s1,
     $                   temp )
            a( k, kwtop-1 ) = temp
            a( k+1, kwtop-1 ) = czero
            k2 = max( kwtop, k-1 )
            CALL zrot( ihi-k2+1, a( k, k2 ), lda, a( k+1, k2 ), lda, c1,
     $                 s1 )
            CALL zrot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1, k-1 ),
     $                 ldb, c1, s1 )
            CALL zrot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1, k+1-kwtop+1 ),
     $                 1, c1, dconjg( s1 ) )
         END DO

*        Chase bulges down
         istartm = kwtop
         istopm = ihi
         k = kwbot-1
         DO WHILE ( k .GE. kwtop )

*           Move bulge down and remove it
            DO k2 = k, kwbot-1
               CALL zlaqz1( .true., .true., k2, kwtop, kwtop+jw-1,
     $                      kwbot, a, lda, b, ldb, jw, kwtop, qc, ldqc,
     $                      jw, kwtop, zc, ldzc )
            END DO

            k = k-1
         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 zgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
     $               a( kwtop, ihi+1 ), lda, czero, work, jw )
         CALL zlacpy( 'ALL', jw, istopm-ihi, work, jw, a( kwtop,
     $                ihi+1 ), lda )
         CALL zgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
     $               b( kwtop, ihi+1 ), ldb, czero, work, jw )
         CALL zlacpy( 'ALL', jw, istopm-ihi, work, jw, b( kwtop,
     $                ihi+1 ), ldb )
      END IF
      IF ( ilq ) THEN
         CALL zgemm( 'N', 'N', n, jw, jw, cone, q( 1, kwtop ), ldq, qc,
     $               ldqc, czero, work, n )
         CALL zlacpy( 'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
      END IF

      IF ( kwtop-1-istartm+1 > 0 ) THEN
         CALL zgemm( 'N', 'N', kwtop-istartm, jw, jw, cone, a( istartm,
     $               kwtop ), lda, zc, ldzc, czero, work,
     $               kwtop-istartm )
        CALL zlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
     $               a( istartm, kwtop ), lda )
         CALL zgemm( 'N', 'N', kwtop-istartm, jw, jw, cone, b( istartm,
     $               kwtop ), ldb, zc, ldzc, czero, work,
     $               kwtop-istartm )
        CALL zlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
     $               b( istartm, kwtop ), ldb )
      END IF
      IF ( ilz ) THEN
         CALL zgemm( 'N', 'N', n, jw, jw, cone, z( 1, kwtop ), ldz, zc,
     $               ldzc, czero, work, n )
         CALL zlacpy( 'ALL', n, jw, work, n, z( 1, kwtop ), ldz )
      END IF

      END SUBROUTINE