slaqz3.f

Go to the documentation of this file.

*> \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 doubleGEcomputational
*>
*  =====================================================================
      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
.EQ..OR..EQ.      IF ( KWTOP  ILO  S  ZERO ) THEN
         KWBOT = KWTOP-1
      ELSE
         KWBOT = IHI
         K = 1
         K2 = 1
.LE.         DO WHILE ( K  JW )
            BULGE = .FALSE.
.GE.            IF ( KWBOT-KWTOP+1  2 ) THEN
.NE.               BULGE = A( KWBOT, KWBOT-1 )  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 ) ) )
.EQ.               IF( TEMP  ZERO )THEN
                  TEMP = ABS( S )
               END IF
               IF ( MAX( ABS( S*QC( 1, KWBOT-KWTOP ) ), ABS( S*QC( 1,
.LE.     $            KWBOT-KWTOP+1 ) ) )  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 ) )
.EQ.               IF( TEMP  ZERO ) THEN
                  TEMP = ABS( S )
               END IF
.LE.               IF ( ( ABS( S*QC( 1, KWBOT-KWTOP+1 ) ) )  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
.LE.      DO WHILE ( K  IHI )
         BULGE = .FALSE.
.LT.         IF ( K  IHI ) THEN
.NE.            IF ( A( K+1, K )  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
 
.NE..AND..NE.      IF ( KWTOP  ILO  S  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
.GE.         DO WHILE ( K  KWTOP )
.GE..AND..NE.            IF ( ( K  KWTOP+1 )  A( K+1, K-1 )  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