LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlagv2.f
Go to the documentation of this file.
00001 *> \brief \b DLAGV2
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLAGV2 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlagv2.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlagv2.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlagv2.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLAGV2( A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL,
00022 *                          CSR, SNR )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       INTEGER            LDA, LDB
00026 *       DOUBLE PRECISION   CSL, CSR, SNL, SNR
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       DOUBLE PRECISION   A( LDA, * ), ALPHAI( 2 ), ALPHAR( 2 ),
00030 *      $                   B( LDB, * ), BETA( 2 )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> DLAGV2 computes the Generalized Schur factorization of a real 2-by-2
00040 *> matrix pencil (A,B) where B is upper triangular. This routine
00041 *> computes orthogonal (rotation) matrices given by CSL, SNL and CSR,
00042 *> SNR such that
00043 *>
00044 *> 1) if the pencil (A,B) has two real eigenvalues (include 0/0 or 1/0
00045 *>    types), then
00046 *>
00047 *>    [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]
00048 *>    [  0  a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]
00049 *>
00050 *>    [ b11 b12 ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]
00051 *>    [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ],
00052 *>
00053 *> 2) if the pencil (A,B) has a pair of complex conjugate eigenvalues,
00054 *>    then
00055 *>
00056 *>    [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]
00057 *>    [ a21 a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]
00058 *>
00059 *>    [ b11  0  ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]
00060 *>    [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ]
00061 *>
00062 *>    where b11 >= b22 > 0.
00063 *>
00064 *> \endverbatim
00065 *
00066 *  Arguments:
00067 *  ==========
00068 *
00069 *> \param[in,out] A
00070 *> \verbatim
00071 *>          A is DOUBLE PRECISION array, dimension (LDA, 2)
00072 *>          On entry, the 2 x 2 matrix A.
00073 *>          On exit, A is overwritten by the ``A-part'' of the
00074 *>          generalized Schur form.
00075 *> \endverbatim
00076 *>
00077 *> \param[in] LDA
00078 *> \verbatim
00079 *>          LDA is INTEGER
00080 *>          THe leading dimension of the array A.  LDA >= 2.
00081 *> \endverbatim
00082 *>
00083 *> \param[in,out] B
00084 *> \verbatim
00085 *>          B is DOUBLE PRECISION array, dimension (LDB, 2)
00086 *>          On entry, the upper triangular 2 x 2 matrix B.
00087 *>          On exit, B is overwritten by the ``B-part'' of the
00088 *>          generalized Schur form.
00089 *> \endverbatim
00090 *>
00091 *> \param[in] LDB
00092 *> \verbatim
00093 *>          LDB is INTEGER
00094 *>          THe leading dimension of the array B.  LDB >= 2.
00095 *> \endverbatim
00096 *>
00097 *> \param[out] ALPHAR
00098 *> \verbatim
00099 *>          ALPHAR is DOUBLE PRECISION array, dimension (2)
00100 *> \endverbatim
00101 *>
00102 *> \param[out] ALPHAI
00103 *> \verbatim
00104 *>          ALPHAI is DOUBLE PRECISION array, dimension (2)
00105 *> \endverbatim
00106 *>
00107 *> \param[out] BETA
00108 *> \verbatim
00109 *>          BETA is DOUBLE PRECISION array, dimension (2)
00110 *>          (ALPHAR(k)+i*ALPHAI(k))/BETA(k) are the eigenvalues of the
00111 *>          pencil (A,B), k=1,2, i = sqrt(-1).  Note that BETA(k) may
00112 *>          be zero.
00113 *> \endverbatim
00114 *>
00115 *> \param[out] CSL
00116 *> \verbatim
00117 *>          CSL is DOUBLE PRECISION
00118 *>          The cosine of the left rotation matrix.
00119 *> \endverbatim
00120 *>
00121 *> \param[out] SNL
00122 *> \verbatim
00123 *>          SNL is DOUBLE PRECISION
00124 *>          The sine of the left rotation matrix.
00125 *> \endverbatim
00126 *>
00127 *> \param[out] CSR
00128 *> \verbatim
00129 *>          CSR is DOUBLE PRECISION
00130 *>          The cosine of the right rotation matrix.
00131 *> \endverbatim
00132 *>
00133 *> \param[out] SNR
00134 *> \verbatim
00135 *>          SNR is DOUBLE PRECISION
00136 *>          The sine of the right rotation matrix.
00137 *> \endverbatim
00138 *
00139 *  Authors:
00140 *  ========
00141 *
00142 *> \author Univ. of Tennessee 
00143 *> \author Univ. of California Berkeley 
00144 *> \author Univ. of Colorado Denver 
00145 *> \author NAG Ltd. 
00146 *
00147 *> \date November 2011
00148 *
00149 *> \ingroup doubleOTHERauxiliary
00150 *
00151 *> \par Contributors:
00152 *  ==================
00153 *>
00154 *>     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
00155 *
00156 *  =====================================================================
00157       SUBROUTINE DLAGV2( A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL,
00158      $                   CSR, SNR )
00159 *
00160 *  -- LAPACK auxiliary routine (version 3.4.0) --
00161 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00162 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00163 *     November 2011
00164 *
00165 *     .. Scalar Arguments ..
00166       INTEGER            LDA, LDB
00167       DOUBLE PRECISION   CSL, CSR, SNL, SNR
00168 *     ..
00169 *     .. Array Arguments ..
00170       DOUBLE PRECISION   A( LDA, * ), ALPHAI( 2 ), ALPHAR( 2 ),
00171      $                   B( LDB, * ), BETA( 2 )
00172 *     ..
00173 *
00174 *  =====================================================================
00175 *
00176 *     .. Parameters ..
00177       DOUBLE PRECISION   ZERO, ONE
00178       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00179 *     ..
00180 *     .. Local Scalars ..
00181       DOUBLE PRECISION   ANORM, ASCALE, BNORM, BSCALE, H1, H2, H3, QQ,
00182      $                   R, RR, SAFMIN, SCALE1, SCALE2, T, ULP, WI, WR1,
00183      $                   WR2
00184 *     ..
00185 *     .. External Subroutines ..
00186       EXTERNAL           DLAG2, DLARTG, DLASV2, DROT
00187 *     ..
00188 *     .. External Functions ..
00189       DOUBLE PRECISION   DLAMCH, DLAPY2
00190       EXTERNAL           DLAMCH, DLAPY2
00191 *     ..
00192 *     .. Intrinsic Functions ..
00193       INTRINSIC          ABS, MAX
00194 *     ..
00195 *     .. Executable Statements ..
00196 *
00197       SAFMIN = DLAMCH( 'S' )
00198       ULP = DLAMCH( 'P' )
00199 *
00200 *     Scale A
00201 *
00202       ANORM = MAX( ABS( A( 1, 1 ) )+ABS( A( 2, 1 ) ),
00203      $        ABS( A( 1, 2 ) )+ABS( A( 2, 2 ) ), SAFMIN )
00204       ASCALE = ONE / ANORM
00205       A( 1, 1 ) = ASCALE*A( 1, 1 )
00206       A( 1, 2 ) = ASCALE*A( 1, 2 )
00207       A( 2, 1 ) = ASCALE*A( 2, 1 )
00208       A( 2, 2 ) = ASCALE*A( 2, 2 )
00209 *
00210 *     Scale B
00211 *
00212       BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 1, 2 ) )+ABS( B( 2, 2 ) ),
00213      $        SAFMIN )
00214       BSCALE = ONE / BNORM
00215       B( 1, 1 ) = BSCALE*B( 1, 1 )
00216       B( 1, 2 ) = BSCALE*B( 1, 2 )
00217       B( 2, 2 ) = BSCALE*B( 2, 2 )
00218 *
00219 *     Check if A can be deflated
00220 *
00221       IF( ABS( A( 2, 1 ) ).LE.ULP ) THEN
00222          CSL = ONE
00223          SNL = ZERO
00224          CSR = ONE
00225          SNR = ZERO
00226          A( 2, 1 ) = ZERO
00227          B( 2, 1 ) = ZERO
00228          WI = ZERO
00229 *
00230 *     Check if B is singular
00231 *
00232       ELSE IF( ABS( B( 1, 1 ) ).LE.ULP ) THEN
00233          CALL DLARTG( A( 1, 1 ), A( 2, 1 ), CSL, SNL, R )
00234          CSR = ONE
00235          SNR = ZERO
00236          CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL )
00237          CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL )
00238          A( 2, 1 ) = ZERO
00239          B( 1, 1 ) = ZERO
00240          B( 2, 1 ) = ZERO
00241          WI = ZERO
00242 *
00243       ELSE IF( ABS( B( 2, 2 ) ).LE.ULP ) THEN
00244          CALL DLARTG( A( 2, 2 ), A( 2, 1 ), CSR, SNR, T )
00245          SNR = -SNR
00246          CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR )
00247          CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR )
00248          CSL = ONE
00249          SNL = ZERO
00250          A( 2, 1 ) = ZERO
00251          B( 2, 1 ) = ZERO
00252          B( 2, 2 ) = ZERO
00253          WI = ZERO
00254 *
00255       ELSE
00256 *
00257 *        B is nonsingular, first compute the eigenvalues of (A,B)
00258 *
00259          CALL DLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2,
00260      $               WI )
00261 *
00262          IF( WI.EQ.ZERO ) THEN
00263 *
00264 *           two real eigenvalues, compute s*A-w*B
00265 *
00266             H1 = SCALE1*A( 1, 1 ) - WR1*B( 1, 1 )
00267             H2 = SCALE1*A( 1, 2 ) - WR1*B( 1, 2 )
00268             H3 = SCALE1*A( 2, 2 ) - WR1*B( 2, 2 )
00269 *
00270             RR = DLAPY2( H1, H2 )
00271             QQ = DLAPY2( SCALE1*A( 2, 1 ), H3 )
00272 *
00273             IF( RR.GT.QQ ) THEN
00274 *
00275 *              find right rotation matrix to zero 1,1 element of
00276 *              (sA - wB)
00277 *
00278                CALL DLARTG( H2, H1, CSR, SNR, T )
00279 *
00280             ELSE
00281 *
00282 *              find right rotation matrix to zero 2,1 element of
00283 *              (sA - wB)
00284 *
00285                CALL DLARTG( H3, SCALE1*A( 2, 1 ), CSR, SNR, T )
00286 *
00287             END IF
00288 *
00289             SNR = -SNR
00290             CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR )
00291             CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR )
00292 *
00293 *           compute inf norms of A and B
00294 *
00295             H1 = MAX( ABS( A( 1, 1 ) )+ABS( A( 1, 2 ) ),
00296      $           ABS( A( 2, 1 ) )+ABS( A( 2, 2 ) ) )
00297             H2 = MAX( ABS( B( 1, 1 ) )+ABS( B( 1, 2 ) ),
00298      $           ABS( B( 2, 1 ) )+ABS( B( 2, 2 ) ) )
00299 *
00300             IF( ( SCALE1*H1 ).GE.ABS( WR1 )*H2 ) THEN
00301 *
00302 *              find left rotation matrix Q to zero out B(2,1)
00303 *
00304                CALL DLARTG( B( 1, 1 ), B( 2, 1 ), CSL, SNL, R )
00305 *
00306             ELSE
00307 *
00308 *              find left rotation matrix Q to zero out A(2,1)
00309 *
00310                CALL DLARTG( A( 1, 1 ), A( 2, 1 ), CSL, SNL, R )
00311 *
00312             END IF
00313 *
00314             CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL )
00315             CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL )
00316 *
00317             A( 2, 1 ) = ZERO
00318             B( 2, 1 ) = ZERO
00319 *
00320          ELSE
00321 *
00322 *           a pair of complex conjugate eigenvalues
00323 *           first compute the SVD of the matrix B
00324 *
00325             CALL DLASV2( B( 1, 1 ), B( 1, 2 ), B( 2, 2 ), R, T, SNR,
00326      $                   CSR, SNL, CSL )
00327 *
00328 *           Form (A,B) := Q(A,B)Z**T where Q is left rotation matrix and
00329 *           Z is right rotation matrix computed from DLASV2
00330 *
00331             CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL )
00332             CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL )
00333             CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR )
00334             CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR )
00335 *
00336             B( 2, 1 ) = ZERO
00337             B( 1, 2 ) = ZERO
00338 *
00339          END IF
00340 *
00341       END IF
00342 *
00343 *     Unscaling
00344 *
00345       A( 1, 1 ) = ANORM*A( 1, 1 )
00346       A( 2, 1 ) = ANORM*A( 2, 1 )
00347       A( 1, 2 ) = ANORM*A( 1, 2 )
00348       A( 2, 2 ) = ANORM*A( 2, 2 )
00349       B( 1, 1 ) = BNORM*B( 1, 1 )
00350       B( 2, 1 ) = BNORM*B( 2, 1 )
00351       B( 1, 2 ) = BNORM*B( 1, 2 )
00352       B( 2, 2 ) = BNORM*B( 2, 2 )
00353 *
00354       IF( WI.EQ.ZERO ) THEN
00355          ALPHAR( 1 ) = A( 1, 1 )
00356          ALPHAR( 2 ) = A( 2, 2 )
00357          ALPHAI( 1 ) = ZERO
00358          ALPHAI( 2 ) = ZERO
00359          BETA( 1 ) = B( 1, 1 )
00360          BETA( 2 ) = B( 2, 2 )
00361       ELSE
00362          ALPHAR( 1 ) = ANORM*WR1 / SCALE1 / BNORM
00363          ALPHAI( 1 ) = ANORM*WI / SCALE1 / BNORM
00364          ALPHAR( 2 ) = ALPHAR( 1 )
00365          ALPHAI( 2 ) = -ALPHAI( 1 )
00366          BETA( 1 ) = ONE
00367          BETA( 2 ) = ONE
00368       END IF
00369 *
00370       RETURN
00371 *
00372 *     End of DLAGV2
00373 *
00374       END
 All Files Functions