SRC/dlagv2.f

   1 *> \brief \b DLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A,B) where B is upper triangular.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLAGV2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlagv2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlagv2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlagv2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLAGV2( A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL,
  22 *                          CSR, SNR )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            LDA, LDB
  26 *       DOUBLE PRECISION   CSL, CSR, SNL, SNR
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   A( LDA, * ), ALPHAI( 2 ), ALPHAR( 2 ),
  30 *      $                   B( LDB, * ), BETA( 2 )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> DLAGV2 computes the Generalized Schur factorization of a real 2-by-2
  40 *> matrix pencil (A,B) where B is upper triangular. This routine
  41 *> computes orthogonal (rotation) matrices given by CSL, SNL and CSR,
  42 *> SNR such that
  43 *>
  44 *> 1) if the pencil (A,B) has two real eigenvalues (include 0/0 or 1/0
  45 *>    types), then
  46 *>
  47 *>    [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]
  48 *>    [  0  a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]
  49 *>
  50 *>    [ b11 b12 ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]
  51 *>    [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ],
  52 *>
  53 *> 2) if the pencil (A,B) has a pair of complex conjugate eigenvalues,
  54 *>    then
  55 *>
  56 *>    [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]
  57 *>    [ a21 a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]
  58 *>
  59 *>    [ b11  0  ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]
  60 *>    [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ]
  61 *>
  62 *>    where b11 >= b22 > 0.
  63 *>
  64 *> \endverbatim
  65 *
  66 *  Arguments:
  67 *  ==========
  68 *
  69 *> \param[in,out] A
  70 *> \verbatim
  71 *>          A is DOUBLE PRECISION array, dimension (LDA, 2)
  72 *>          On entry, the 2 x 2 matrix A.
  73 *>          On exit, A is overwritten by the ``A-part'' of the
  74 *>          generalized Schur form.
  75 *> \endverbatim
  76 *>
  77 *> \param[in] LDA
  78 *> \verbatim
  79 *>          LDA is INTEGER
  80 *>          THe leading dimension of the array A.  LDA >= 2.
  81 *> \endverbatim
  82 *>
  83 *> \param[in,out] B
  84 *> \verbatim
  85 *>          B is DOUBLE PRECISION array, dimension (LDB, 2)
  86 *>          On entry, the upper triangular 2 x 2 matrix B.
  87 *>          On exit, B is overwritten by the ``B-part'' of the
  88 *>          generalized Schur form.
  89 *> \endverbatim
  90 *>
  91 *> \param[in] LDB
  92 *> \verbatim
  93 *>          LDB is INTEGER
  94 *>          THe leading dimension of the array B.  LDB >= 2.
  95 *> \endverbatim
  96 *>
  97 *> \param[out] ALPHAR
  98 *> \verbatim
  99 *>          ALPHAR is DOUBLE PRECISION array, dimension (2)
 100 *> \endverbatim
 101 *>
 102 *> \param[out] ALPHAI
 103 *> \verbatim
 104 *>          ALPHAI is DOUBLE PRECISION array, dimension (2)
 105 *> \endverbatim
 106 *>
 107 *> \param[out] BETA
 108 *> \verbatim
 109 *>          BETA is DOUBLE PRECISION array, dimension (2)
 110 *>          (ALPHAR(k)+i*ALPHAI(k))/BETA(k) are the eigenvalues of the
 111 *>          pencil (A,B), k=1,2, i = sqrt(-1).  Note that BETA(k) may
 112 *>          be zero.
 113 *> \endverbatim
 114 *>
 115 *> \param[out] CSL
 116 *> \verbatim
 117 *>          CSL is DOUBLE PRECISION
 118 *>          The cosine of the left rotation matrix.
 119 *> \endverbatim
 120 *>
 121 *> \param[out] SNL
 122 *> \verbatim
 123 *>          SNL is DOUBLE PRECISION
 124 *>          The sine of the left rotation matrix.
 125 *> \endverbatim
 126 *>
 127 *> \param[out] CSR
 128 *> \verbatim
 129 *>          CSR is DOUBLE PRECISION
 130 *>          The cosine of the right rotation matrix.
 131 *> \endverbatim
 132 *>
 133 *> \param[out] SNR
 134 *> \verbatim
 135 *>          SNR is DOUBLE PRECISION
 136 *>          The sine of the right rotation matrix.
 137 *> \endverbatim
 138 *
 139 *  Authors:
 140 *  ========
 141 *
 142 *> \author Univ. of Tennessee
 143 *> \author Univ. of California Berkeley
 144 *> \author Univ. of Colorado Denver
 145 *> \author NAG Ltd.
 146 *
 147 *> \date September 2012
 148 *
 149 *> \ingroup doubleOTHERauxiliary
 150 *
 151 *> \par Contributors:
 152 *  ==================
 153 *>
 154 *>     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
 155 *
 156 *  =====================================================================
 157       SUBROUTINE DLAGV2( A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL,
 158      $                   CSR, SNR )
 159 *
 160 *  -- LAPACK auxiliary routine (version 3.4.2) --
 161 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 162 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 163 *     September 2012
 164 *
 165 *     .. Scalar Arguments ..
 166       INTEGER            LDA, LDB
 167       DOUBLE PRECISION   CSL, CSR, SNL, SNR
 168 *     ..
 169 *     .. Array Arguments ..
 170       DOUBLE PRECISION   A( LDA, * ), ALPHAI( 2 ), ALPHAR( 2 ),
 171      $                   B( LDB, * ), BETA( 2 )
 172 *     ..
 173 *
 174 *  =====================================================================
 175 *
 176 *     .. Parameters ..
 177       DOUBLE PRECISION   ZERO, ONE
 178       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 179 *     ..
 180 *     .. Local Scalars ..
 181       DOUBLE PRECISION   ANORM, ASCALE, BNORM, BSCALE, H1, H2, H3, QQ,
 182      $                   R, RR, SAFMIN, SCALE1, SCALE2, T, ULP, WI, WR1,
 183      $                   WR2
 184 *     ..
 185 *     .. External Subroutines ..
 186       EXTERNAL           DLAG2, DLARTG, DLASV2, DROT
 187 *     ..
 188 *     .. External Functions ..
 189       DOUBLE PRECISION   DLAMCH, DLAPY2
 190       EXTERNAL           DLAMCH, DLAPY2
 191 *     ..
 192 *     .. Intrinsic Functions ..
 193       INTRINSIC          ABS, MAX
 194 *     ..
 195 *     .. Executable Statements ..
 196 *
 197       SAFMIN = DLAMCH( 'S' )
 198       ULP = DLAMCH( 'P' )
 199 *
 200 *     Scale A
 201 *
 202       ANORM = MAX( ABS( A( 1, 1 ) )+ABS( A( 2, 1 ) ),
 203      $        ABS( A( 1, 2 ) )+ABS( A( 2, 2 ) ), SAFMIN )
 204       ASCALE = ONE / ANORM
 205       A( 1, 1 ) = ASCALE*A( 1, 1 )
 206       A( 1, 2 ) = ASCALE*A( 1, 2 )
 207       A( 2, 1 ) = ASCALE*A( 2, 1 )
 208       A( 2, 2 ) = ASCALE*A( 2, 2 )
 209 *
 210 *     Scale B
 211 *
 212       BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 1, 2 ) )+ABS( B( 2, 2 ) ),
 213      $        SAFMIN )
 214       BSCALE = ONE / BNORM
 215       B( 1, 1 ) = BSCALE*B( 1, 1 )
 216       B( 1, 2 ) = BSCALE*B( 1, 2 )
 217       B( 2, 2 ) = BSCALE*B( 2, 2 )
 218 *
 219 *     Check if A can be deflated
 220 *
 221       IF( ABS( A( 2, 1 ) ).LE.ULP ) THEN
 222          CSL = ONE
 223          SNL = ZERO
 224          CSR = ONE
 225          SNR = ZERO
 226          A( 2, 1 ) = ZERO
 227          B( 2, 1 ) = ZERO
 228          WI = ZERO
 229 *
 230 *     Check if B is singular
 231 *
 232       ELSE IF( ABS( B( 1, 1 ) ).LE.ULP ) THEN
 233          CALL DLARTG( A( 1, 1 ), A( 2, 1 ), CSL, SNL, R )
 234          CSR = ONE
 235          SNR = ZERO
 236          CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL )
 237          CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL )
 238          A( 2, 1 ) = ZERO
 239          B( 1, 1 ) = ZERO
 240          B( 2, 1 ) = ZERO
 241          WI = ZERO
 242 *
 243       ELSE IF( ABS( B( 2, 2 ) ).LE.ULP ) THEN
 244          CALL DLARTG( A( 2, 2 ), A( 2, 1 ), CSR, SNR, T )
 245          SNR = -SNR
 246          CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR )
 247          CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR )
 248          CSL = ONE
 249          SNL = ZERO
 250          A( 2, 1 ) = ZERO
 251          B( 2, 1 ) = ZERO
 252          B( 2, 2 ) = ZERO
 253          WI = ZERO
 254 *
 255       ELSE
 256 *
 257 *        B is nonsingular, first compute the eigenvalues of (A,B)
 258 *
 259          CALL DLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2,
 260      $               WI )
 261 *
 262          IF( WI.EQ.ZERO ) THEN
 263 *
 264 *           two real eigenvalues, compute s*A-w*B
 265 *
 266             H1 = SCALE1*A( 1, 1 ) - WR1*B( 1, 1 )
 267             H2 = SCALE1*A( 1, 2 ) - WR1*B( 1, 2 )
 268             H3 = SCALE1*A( 2, 2 ) - WR1*B( 2, 2 )
 269 *
 270             RR = DLAPY2( H1, H2 )
 271             QQ = DLAPY2( SCALE1*A( 2, 1 ), H3 )
 272 *
 273             IF( RR.GT.QQ ) THEN
 274 *
 275 *              find right rotation matrix to zero 1,1 element of
 276 *              (sA - wB)
 277 *
 278                CALL DLARTG( H2, H1, CSR, SNR, T )
 279 *
 280             ELSE
 281 *
 282 *              find right rotation matrix to zero 2,1 element of
 283 *              (sA - wB)
 284 *
 285                CALL DLARTG( H3, SCALE1*A( 2, 1 ), CSR, SNR, T )
 286 *
 287             END IF
 288 *
 289             SNR = -SNR
 290             CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR )
 291             CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR )
 292 *
 293 *           compute inf norms of A and B
 294 *
 295             H1 = MAX( ABS( A( 1, 1 ) )+ABS( A( 1, 2 ) ),
 296      $           ABS( A( 2, 1 ) )+ABS( A( 2, 2 ) ) )
 297             H2 = MAX( ABS( B( 1, 1 ) )+ABS( B( 1, 2 ) ),
 298      $           ABS( B( 2, 1 ) )+ABS( B( 2, 2 ) ) )
 299 *
 300             IF( ( SCALE1*H1 ).GE.ABS( WR1 )*H2 ) THEN
 301 *
 302 *              find left rotation matrix Q to zero out B(2,1)
 303 *
 304                CALL DLARTG( B( 1, 1 ), B( 2, 1 ), CSL, SNL, R )
 305 *
 306             ELSE
 307 *
 308 *              find left rotation matrix Q to zero out A(2,1)
 309 *
 310                CALL DLARTG( A( 1, 1 ), A( 2, 1 ), CSL, SNL, R )
 311 *
 312             END IF
 313 *
 314             CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL )
 315             CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL )
 316 *
 317             A( 2, 1 ) = ZERO
 318             B( 2, 1 ) = ZERO
 319 *
 320          ELSE
 321 *
 322 *           a pair of complex conjugate eigenvalues
 323 *           first compute the SVD of the matrix B
 324 *
 325             CALL DLASV2( B( 1, 1 ), B( 1, 2 ), B( 2, 2 ), R, T, SNR,
 326      $                   CSR, SNL, CSL )
 327 *
 328 *           Form (A,B) := Q(A,B)Z**T where Q is left rotation matrix and
 329 *           Z is right rotation matrix computed from DLASV2
 330 *
 331             CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL )
 332             CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL )
 333             CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR )
 334             CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR )
 335 *
 336             B( 2, 1 ) = ZERO
 337             B( 1, 2 ) = ZERO
 338 *
 339          END IF
 340 *
 341       END IF
 342 *
 343 *     Unscaling
 344 *
 345       A( 1, 1 ) = ANORM*A( 1, 1 )
 346       A( 2, 1 ) = ANORM*A( 2, 1 )
 347       A( 1, 2 ) = ANORM*A( 1, 2 )
 348       A( 2, 2 ) = ANORM*A( 2, 2 )
 349       B( 1, 1 ) = BNORM*B( 1, 1 )
 350       B( 2, 1 ) = BNORM*B( 2, 1 )
 351       B( 1, 2 ) = BNORM*B( 1, 2 )
 352       B( 2, 2 ) = BNORM*B( 2, 2 )
 353 *
 354       IF( WI.EQ.ZERO ) THEN
 355          ALPHAR( 1 ) = A( 1, 1 )
 356          ALPHAR( 2 ) = A( 2, 2 )
 357          ALPHAI( 1 ) = ZERO
 358          ALPHAI( 2 ) = ZERO
 359          BETA( 1 ) = B( 1, 1 )
 360          BETA( 2 ) = B( 2, 2 )
 361       ELSE
 362          ALPHAR( 1 ) = ANORM*WR1 / SCALE1 / BNORM
 363          ALPHAI( 1 ) = ANORM*WI / SCALE1 / BNORM
 364          ALPHAR( 2 ) = ALPHAR( 1 )
 365          ALPHAI( 2 ) = -ALPHAI( 1 )
 366          BETA( 1 ) = ONE
 367          BETA( 2 ) = ONE
 368       END IF
 369 *
 370       RETURN
 371 *
 372 *     End of DLAGV2
 373 *
 374       END