dlagv2.f

Go to the documentation of this file.

*> \brief \b DLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A,B) where B is upper triangular.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*> \htmlonly
*> Download DLAGV2 + dependencies 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlagv2.f"> 
*> [TGZ]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlagv2.f"> 
*> [ZIP]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlagv2.f"> 
*> [TXT]</a>
*> \endhtmlonly 
*
*  Definition:
*  ===========
*
*       SUBROUTINE DLAGV2( A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL,
*                          CSR, SNR )
* 
*       .. Scalar Arguments ..
*       INTEGER            LDA, LDB
*       DOUBLE PRECISION   CSL, CSR, SNL, SNR
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   A( LDA, * ), ALPHAI( 2 ), ALPHAR( 2 ),
*      $                   B( LDB, * ), BETA( 2 )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DLAGV2 computes the Generalized Schur factorization of a real 2-by-2
*> matrix pencil (A,B) where B is upper triangular. This routine
*> computes orthogonal (rotation) matrices given by CSL, SNL and CSR,
*> SNR such that
*>
*> 1) if the pencil (A,B) has two real eigenvalues (include 0/0 or 1/0
*>    types), then
*>
*>    [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]
*>    [  0  a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]
*>
*>    [ b11 b12 ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]
*>    [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ],
*>
*> 2) if the pencil (A,B) has a pair of complex conjugate eigenvalues,
*>    then
*>
*>    [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]
*>    [ a21 a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]
*>
*>    [ b11  0  ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]
*>    [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ]
*>
*>    where b11 >= b22 > 0.
*>
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in,out] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (LDA, 2)
*>          On entry, the 2 x 2 matrix A.
*>          On exit, A is overwritten by the ``A-part'' of the
*>          generalized Schur form.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          THe leading dimension of the array A.  LDA >= 2.
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*>          B is DOUBLE PRECISION array, dimension (LDB, 2)
*>          On entry, the upper triangular 2 x 2 matrix B.
*>          On exit, B is overwritten by the ``B-part'' of the
*>          generalized Schur form.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>          THe leading dimension of the array B.  LDB >= 2.
*> \endverbatim
*>
*> \param[out] ALPHAR
*> \verbatim
*>          ALPHAR is DOUBLE PRECISION array, dimension (2)
*> \endverbatim
*>
*> \param[out] ALPHAI
*> \verbatim
*>          ALPHAI is DOUBLE PRECISION array, dimension (2)
*> \endverbatim
*>
*> \param[out] BETA
*> \verbatim
*>          BETA is DOUBLE PRECISION array, dimension (2)
*>          (ALPHAR(k)+i*ALPHAI(k))/BETA(k) are the eigenvalues of the
*>          pencil (A,B), k=1,2, i = sqrt(-1).  Note that BETA(k) may
*>          be zero.
*> \endverbatim
*>
*> \param[out] CSL
*> \verbatim
*>          CSL is DOUBLE PRECISION
*>          The cosine of the left rotation matrix.
*> \endverbatim
*>
*> \param[out] SNL
*> \verbatim
*>          SNL is DOUBLE PRECISION
*>          The sine of the left rotation matrix.
*> \endverbatim
*>
*> \param[out] CSR
*> \verbatim
*>          CSR is DOUBLE PRECISION
*>          The cosine of the right rotation matrix.
*> \endverbatim
*>
*> \param[out] SNR
*> \verbatim
*>          SNR is DOUBLE PRECISION
*>          The sine of the right rotation matrix.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date September 2012
*
*> \ingroup doubleOTHERauxiliary
*
*> \par Contributors:
*  ==================
*>
*>     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
*
*  =====================================================================
      SUBROUTINE dlagv2( A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL,
     $                   csr, snr )
*
*  -- LAPACK auxiliary routine (version 3.4.2) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     September 2012
*
*     .. Scalar Arguments ..
      INTEGER            LDA, LDB
      DOUBLE PRECISION   CSL, CSR, SNL, SNR
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   A( lda, * ), ALPHAI( 2 ), ALPHAR( 2 ),
     $                   b( ldb, * ), beta( 2 )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      parameter                ( zero = 0.0d+0, one = 1.0d+0 )
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION   ANORM, ASCALE, BNORM, BSCALE, H1, H2, H3, QQ,
     $                   r, rr, safmin, scale1, scale2, t, ulp, wi, wr1,
     $                   wr2
*     ..
*     .. External Subroutines ..
      EXTERNAL           dlag2, dlartg, dlasv2, drot
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH, DLAPY2
      EXTERNAL           dlamch, dlapy2
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max
*     ..
*     .. Executable Statements ..
*
      safmin = dlamch( 'S' )
      ulp = dlamch( 'P' )
*
*     Scale A
*
      anorm = max( abs( a( 1, 1 ) )+abs( a( 2, 1 ) ),
     $        abs( a( 1, 2 ) )+abs( a( 2, 2 ) ), safmin )
      ascale = one / anorm
      a( 1, 1 ) = ascale*a( 1, 1 )
      a( 1, 2 ) = ascale*a( 1, 2 )
      a( 2, 1 ) = ascale*a( 2, 1 )
      a( 2, 2 ) = ascale*a( 2, 2 )
*
*     Scale B
*
      bnorm = max( abs( b( 1, 1 ) ), abs( b( 1, 2 ) )+abs( b( 2, 2 ) ),
     $        safmin )
      bscale = one / bnorm
      b( 1, 1 ) = bscale*b( 1, 1 )
      b( 1, 2 ) = bscale*b( 1, 2 )
      b( 2, 2 ) = bscale*b( 2, 2 )
*
*     Check if A can be deflated
*
      IF( abs( a( 2, 1 ) ).LE.ulp ) THEN
         csl = one
         snl = zero
         csr = one
         snr = zero
         a( 2, 1 ) = zero
         b( 2, 1 ) = zero
         wi = zero
*
*     Check if B is singular
*
      ELSE IF( abs( b( 1, 1 ) ).LE.ulp ) THEN
         CALL dlartg( a( 1, 1 ), a( 2, 1 ), csl, snl, r )
         csr = one
         snr = zero
         CALL drot( 2, a( 1, 1 ), lda, a( 2, 1 ), lda, csl, snl )
         CALL drot( 2, b( 1, 1 ), ldb, b( 2, 1 ), ldb, csl, snl )
         a( 2, 1 ) = zero
         b( 1, 1 ) = zero
         b( 2, 1 ) = zero
         wi = zero
*
      ELSE IF( abs( b( 2, 2 ) ).LE.ulp ) THEN
         CALL dlartg( a( 2, 2 ), a( 2, 1 ), csr, snr, t )
         snr = -snr
         CALL drot( 2, a( 1, 1 ), 1, a( 1, 2 ), 1, csr, snr )
         CALL drot( 2, b( 1, 1 ), 1, b( 1, 2 ), 1, csr, snr )
         csl = one
         snl = zero
         a( 2, 1 ) = zero
         b( 2, 1 ) = zero
         b( 2, 2 ) = zero
         wi = zero
*
      ELSE
*
*        B is nonsingular, first compute the eigenvalues of (A,B)
*
         CALL dlag2( a, lda, b, ldb, safmin, scale1, scale2, wr1, wr2,
     $               wi )
*
         IF( wi.EQ.zero ) THEN
*
*           two real eigenvalues, compute s*A-w*B
*
            h1 = scale1*a( 1, 1 ) - wr1*b( 1, 1 )
            h2 = scale1*a( 1, 2 ) - wr1*b( 1, 2 )
            h3 = scale1*a( 2, 2 ) - wr1*b( 2, 2 )
*
            rr = dlapy2( h1, h2 )
            qq = dlapy2( scale1*a( 2, 1 ), h3 )
*
            IF( rr.GT.qq ) THEN
*
*              find right rotation matrix to zero 1,1 element of
*              (sA - wB)
*
               CALL dlartg( h2, h1, csr, snr, t )
*
            ELSE
*
*              find right rotation matrix to zero 2,1 element of
*              (sA - wB)
*
               CALL dlartg( h3, scale1*a( 2, 1 ), csr, snr, t )
*
            END IF
*
            snr = -snr
            CALL drot( 2, a( 1, 1 ), 1, a( 1, 2 ), 1, csr, snr )
            CALL drot( 2, b( 1, 1 ), 1, b( 1, 2 ), 1, csr, snr )
*
*           compute inf norms of A and B
*
            h1 = max( abs( a( 1, 1 ) )+abs( a( 1, 2 ) ),
     $           abs( a( 2, 1 ) )+abs( a( 2, 2 ) ) )
            h2 = max( abs( b( 1, 1 ) )+abs( b( 1, 2 ) ),
     $           abs( b( 2, 1 ) )+abs( b( 2, 2 ) ) )
*
            IF( ( scale1*h1 ).GE.abs( wr1 )*h2 ) THEN
*
*              find left rotation matrix Q to zero out B(2,1)
*
               CALL dlartg( b( 1, 1 ), b( 2, 1 ), csl, snl, r )
*
            ELSE
*
*              find left rotation matrix Q to zero out A(2,1)
*
               CALL dlartg( a( 1, 1 ), a( 2, 1 ), csl, snl, r )
*
            END IF
*
            CALL drot( 2, a( 1, 1 ), lda, a( 2, 1 ), lda, csl, snl )
            CALL drot( 2, b( 1, 1 ), ldb, b( 2, 1 ), ldb, csl, snl )
*
            a( 2, 1 ) = zero
            b( 2, 1 ) = zero
*
         ELSE
*
*           a pair of complex conjugate eigenvalues
*           first compute the SVD of the matrix B
*
            CALL dlasv2( b( 1, 1 ), b( 1, 2 ), b( 2, 2 ), r, t, snr,
     $                   csr, snl, csl )
*
*           Form (A,B) := Q(A,B)Z**T where Q is left rotation matrix and
*           Z is right rotation matrix computed from DLASV2
*
            CALL drot( 2, a( 1, 1 ), lda, a( 2, 1 ), lda, csl, snl )
            CALL drot( 2, b( 1, 1 ), ldb, b( 2, 1 ), ldb, csl, snl )
            CALL drot( 2, a( 1, 1 ), 1, a( 1, 2 ), 1, csr, snr )
            CALL drot( 2, b( 1, 1 ), 1, b( 1, 2 ), 1, csr, snr )
*
            b( 2, 1 ) = zero
            b( 1, 2 ) = zero
*
         END IF
*
      END IF
*
*     Unscaling
*
      a( 1, 1 ) = anorm*a( 1, 1 )
      a( 2, 1 ) = anorm*a( 2, 1 )
      a( 1, 2 ) = anorm*a( 1, 2 )
      a( 2, 2 ) = anorm*a( 2, 2 )
      b( 1, 1 ) = bnorm*b( 1, 1 )
      b( 2, 1 ) = bnorm*b( 2, 1 )
      b( 1, 2 ) = bnorm*b( 1, 2 )
      b( 2, 2 ) = bnorm*b( 2, 2 )
*
      IF( wi.EQ.zero ) THEN
         alphar( 1 ) = a( 1, 1 )
         alphar( 2 ) = a( 2, 2 )
         alphai( 1 ) = zero
         alphai( 2 ) = zero
         beta( 1 ) = b( 1, 1 )
         beta( 2 ) = b( 2, 2 )
      ELSE
         alphar( 1 ) = anorm*wr1 / scale1 / bnorm
         alphai( 1 ) = anorm*wi / scale1 / bnorm
         alphar( 2 ) = alphar( 1 )
         alphai( 2 ) = -alphai( 1 )
         beta( 1 ) = one
         beta( 2 ) = one
      END IF
*
      RETURN
*
*     End of DLAGV2
*
      END