dlaqz1.f

Go to the documentation of this file.

*> \brief \b DLAQZ1
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DLAQZ1 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqz1.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqz1.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqz1.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*      SUBROUTINE DLAQZ1( A, LDA, B, LDB, SR1, SR2, SI, BETA1, BETA2,
*     $    V )
*      IMPLICIT NONE
*
*      Arguments
*      INTEGER, INTENT( IN ) :: LDA, LDB
*      DOUBLE PRECISION, INTENT( IN ) :: A( LDA, * ), B( LDB, * ), SR1,
*     $                  SR2, SI, BETA1, BETA2
*      DOUBLE PRECISION, INTENT( OUT ) :: V( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*>      Given a 3-by-3 matrix pencil (A,B), DLAQZ1 sets v to a
*>      scalar multiple of the first column of the product
*>
*>      (*)  K = (A - (beta2*sr2 - i*si)*B)*B^(-1)*(beta1*A - (sr2 + i*si2)*B)*B^(-1).
*>
*>      It is assumed that either
*>
*>              1) sr1 = sr2
*>          or
*>              2) si = 0.
*>
*>      This is useful for starting double implicit shift bulges
*>      in the QZ algorithm.
*> \endverbatim
*
*
*  Arguments:
*  ==========
*
*> \param[in] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (LDA,N)
*>              The 3-by-3 matrix A in (*).
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>              The leading dimension of A as declared in
*>              the calling procedure.
*> \endverbatim
*
*> \param[in] B
*> \verbatim
*>          B is DOUBLE PRECISION array, dimension (LDB,N)
*>              The 3-by-3 matrix B in (*).
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>              The leading dimension of B as declared in
*>              the calling procedure.
*> \endverbatim
*>
*> \param[in] SR1
*> \verbatim
*>          SR1 is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in] SR2
*> \verbatim
*>          SR2 is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in] SI
*> \verbatim
*>          SI is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in] BETA1
*> \verbatim
*>          BETA1 is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in] BETA2
*> \verbatim
*>          BETA2 is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[out] V
*> \verbatim
*>          V is DOUBLE PRECISION array, dimension (N)
*>              A scalar multiple of the first column of the
*>              matrix K in (*).
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Thijs Steel, KU Leuven
*
*> \date May 2020
*
*> \ingroup laqz1
*>
*  =====================================================================
      SUBROUTINE dlaqz1( A, LDA, B, LDB, SR1, SR2, SI, BETA1, BETA2,
     $                   V )
      IMPLICIT NONE
*
*     Arguments
      INTEGER, INTENT( IN ) :: LDA, LDB
      DOUBLE PRECISION, INTENT( IN ) :: A( LDA, * ), B( LDB, * ), SR1,
     $                  sr2, si, beta1, beta2
      DOUBLE PRECISION, INTENT( OUT ) :: V( * )
*
*     Parameters
      DOUBLE PRECISION :: ZERO, ONE, HALF
      parameter( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
*
*     Local scalars
      DOUBLE PRECISION :: W( 2 ), SAFMIN, SAFMAX, SCALE1, SCALE2
*
*     External Functions
      DOUBLE PRECISION, EXTERNAL :: DLAMCH
      LOGICAL, EXTERNAL :: DISNAN
*
      safmin = dlamch( 'SAFE MINIMUM' )
      safmax = one/safmin
*
*     Calculate first shifted vector
*
      w( 1 ) = beta1*a( 1, 1 )-sr1*b( 1, 1 )
      w( 2 ) = beta1*a( 2, 1 )-sr1*b( 2, 1 )
      scale1 = sqrt( abs( w( 1 ) ) ) * sqrt( abs( w( 2 ) ) )
      IF( scale1 .GE. safmin .AND. scale1 .LE. safmax ) THEN
         w( 1 ) = w( 1 )/scale1
         w( 2 ) = w( 2 )/scale1
      END IF
*
*     Solve linear system
*
      w( 2 ) = w( 2 )/b( 2, 2 )
      w( 1 ) = ( w( 1 )-b( 1, 2 )*w( 2 ) )/b( 1, 1 )
      scale2 = sqrt( abs( w( 1 ) ) ) * sqrt( abs( w( 2 ) ) )
      IF( scale2 .GE. safmin .AND. scale2 .LE. safmax ) THEN
         w( 1 ) = w( 1 )/scale2
         w( 2 ) = w( 2 )/scale2
      END IF
*
*     Apply second shift
*
      v( 1 ) = beta2*( a( 1, 1 )*w( 1 )+a( 1, 2 )*w( 2 ) )-sr2*( b( 1,
     $   1 )*w( 1 )+b( 1, 2 )*w( 2 ) )
      v( 2 ) = beta2*( a( 2, 1 )*w( 1 )+a( 2, 2 )*w( 2 ) )-sr2*( b( 2,
     $   1 )*w( 1 )+b( 2, 2 )*w( 2 ) )
      v( 3 ) = beta2*( a( 3, 1 )*w( 1 )+a( 3, 2 )*w( 2 ) )-sr2*( b( 3,
     $   1 )*w( 1 )+b( 3, 2 )*w( 2 ) )
*
*     Account for imaginary part
*
      v( 1 ) = v( 1 )+si*si*b( 1, 1 )/scale1/scale2
*
*     Check for overflow
*
      IF( abs( v( 1 ) ).GT.safmax .OR. abs( v( 2 ) ) .GT. safmax .OR.
     $   abs( v( 3 ) ).GT.safmax .OR. disnan( v( 1 ) ) .OR.
     $   disnan( v( 2 ) ) .OR. disnan( v( 3 ) ) ) THEN
         v( 1 ) = zero
         v( 2 ) = zero
         v( 3 ) = zero
      END IF
*
*     End of DLAQZ1
*
      END SUBROUTINE