LAPACK: SRC/dlagv2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLAGV2( A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL,
00002      $                   CSR, SNR )
00003 *
00004 *  -- LAPACK auxiliary routine (version 3.2.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     June 2010
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            LDA, LDB
00011       DOUBLE PRECISION   CSL, CSR, SNL, SNR
00012 *     ..
00013 *     .. Array Arguments ..
00014       DOUBLE PRECISION   A( LDA, * ), ALPHAI( 2 ), ALPHAR( 2 ),
00015      $                   B( LDB, * ), BETA( 2 )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  DLAGV2 computes the Generalized Schur factorization of a real 2-by-2
00022 *  matrix pencil (A,B) where B is upper triangular. This routine
00023 *  computes orthogonal (rotation) matrices given by CSL, SNL and CSR,
00024 *  SNR such that
00025 *
00026 *  1) if the pencil (A,B) has two real eigenvalues (include 0/0 or 1/0
00027 *     types), then
00028 *
00029 *     [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]
00030 *     [  0  a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]
00031 *
00032 *     [ b11 b12 ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]
00033 *     [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ],
00034 *
00035 *  2) if the pencil (A,B) has a pair of complex conjugate eigenvalues,
00036 *     then
00037 *
00038 *     [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]
00039 *     [ a21 a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]
00040 *
00041 *     [ b11  0  ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]
00042 *     [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ]
00043 *
00044 *     where b11 >= b22 > 0.
00045 *
00046 *
00047 *  Arguments
00048 *  =========
00049 *
00050 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA, 2)
00051 *          On entry, the 2 x 2 matrix A.
00052 *          On exit, A is overwritten by the ``A-part'' of the
00053 *          generalized Schur form.
00054 *
00055 *  LDA     (input) INTEGER
00056 *          THe leading dimension of the array A.  LDA >= 2.
00057 *
00058 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB, 2)
00059 *          On entry, the upper triangular 2 x 2 matrix B.
00060 *          On exit, B is overwritten by the ``B-part'' of the
00061 *          generalized Schur form.
00062 *
00063 *  LDB     (input) INTEGER
00064 *          THe leading dimension of the array B.  LDB >= 2.
00065 *
00066 *  ALPHAR  (output) DOUBLE PRECISION array, dimension (2)
00067 *  ALPHAI  (output) DOUBLE PRECISION array, dimension (2)
00068 *  BETA    (output) DOUBLE PRECISION array, dimension (2)
00069 *          (ALPHAR(k)+i*ALPHAI(k))/BETA(k) are the eigenvalues of the
00070 *          pencil (A,B), k=1,2, i = sqrt(-1).  Note that BETA(k) may
00071 *          be zero.
00072 *
00073 *  CSL     (output) DOUBLE PRECISION
00074 *          The cosine of the left rotation matrix.
00075 *
00076 *  SNL     (output) DOUBLE PRECISION
00077 *          The sine of the left rotation matrix.
00078 *
00079 *  CSR     (output) DOUBLE PRECISION
00080 *          The cosine of the right rotation matrix.
00081 *
00082 *  SNR     (output) DOUBLE PRECISION
00083 *          The sine of the right rotation matrix.
00084 *
00085 *  Further Details
00086 *  ===============
00087 *
00088 *  Based on contributions by
00089 *     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
00090 *
00091 *  =====================================================================
00092 *
00093 *     .. Parameters ..
00094       DOUBLE PRECISION   ZERO, ONE
00095       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00096 *     ..
00097 *     .. Local Scalars ..
00098       DOUBLE PRECISION   ANORM, ASCALE, BNORM, BSCALE, H1, H2, H3, QQ,
00099      $                   R, RR, SAFMIN, SCALE1, SCALE2, T, ULP, WI, WR1,
00100      $                   WR2
00101 *     ..
00102 *     .. External Subroutines ..
00103       EXTERNAL           DLAG2, DLARTG, DLASV2, DROT
00104 *     ..
00105 *     .. External Functions ..
00106       DOUBLE PRECISION   DLAMCH, DLAPY2
00107       EXTERNAL           DLAMCH, DLAPY2
00108 *     ..
00109 *     .. Intrinsic Functions ..
00110       INTRINSIC          ABS, MAX
00111 *     ..
00112 *     .. Executable Statements ..
00113 *
00114       SAFMIN = DLAMCH( 'S' )
00115       ULP = DLAMCH( 'P' )
00116 *
00117 *     Scale A
00118 *
00119       ANORM = MAX( ABS( A( 1, 1 ) )+ABS( A( 2, 1 ) ),
00120      $        ABS( A( 1, 2 ) )+ABS( A( 2, 2 ) ), SAFMIN )
00121       ASCALE = ONE / ANORM
00122       A( 1, 1 ) = ASCALE*A( 1, 1 )
00123       A( 1, 2 ) = ASCALE*A( 1, 2 )
00124       A( 2, 1 ) = ASCALE*A( 2, 1 )
00125       A( 2, 2 ) = ASCALE*A( 2, 2 )
00126 *
00127 *     Scale B
00128 *
00129       BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 1, 2 ) )+ABS( B( 2, 2 ) ),
00130      $        SAFMIN )
00131       BSCALE = ONE / BNORM
00132       B( 1, 1 ) = BSCALE*B( 1, 1 )
00133       B( 1, 2 ) = BSCALE*B( 1, 2 )
00134       B( 2, 2 ) = BSCALE*B( 2, 2 )
00135 *
00136 *     Check if A can be deflated
00137 *
00138       IF( ABS( A( 2, 1 ) ).LE.ULP ) THEN
00139          CSL = ONE
00140          SNL = ZERO
00141          CSR = ONE
00142          SNR = ZERO
00143          A( 2, 1 ) = ZERO
00144          B( 2, 1 ) = ZERO
00145          WI = ZERO
00146 *
00147 *     Check if B is singular
00148 *
00149       ELSE IF( ABS( B( 1, 1 ) ).LE.ULP ) THEN
00150          CALL DLARTG( A( 1, 1 ), A( 2, 1 ), CSL, SNL, R )
00151          CSR = ONE
00152          SNR = ZERO
00153          CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL )
00154          CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL )
00155          A( 2, 1 ) = ZERO
00156          B( 1, 1 ) = ZERO
00157          B( 2, 1 ) = ZERO
00158          WI = ZERO
00159 *
00160       ELSE IF( ABS( B( 2, 2 ) ).LE.ULP ) THEN
00161          CALL DLARTG( A( 2, 2 ), A( 2, 1 ), CSR, SNR, T )
00162          SNR = -SNR
00163          CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR )
00164          CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR )
00165          CSL = ONE
00166          SNL = ZERO
00167          A( 2, 1 ) = ZERO
00168          B( 2, 1 ) = ZERO
00169          B( 2, 2 ) = ZERO
00170          WI = ZERO
00171 *
00172       ELSE
00173 *
00174 *        B is nonsingular, first compute the eigenvalues of (A,B)
00175 *
00176          CALL DLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2,
00177      $               WI )
00178 *
00179          IF( WI.EQ.ZERO ) THEN
00180 *
00181 *           two real eigenvalues, compute s*A-w*B
00182 *
00183             H1 = SCALE1*A( 1, 1 ) - WR1*B( 1, 1 )
00184             H2 = SCALE1*A( 1, 2 ) - WR1*B( 1, 2 )
00185             H3 = SCALE1*A( 2, 2 ) - WR1*B( 2, 2 )
00186 *
00187             RR = DLAPY2( H1, H2 )
00188             QQ = DLAPY2( SCALE1*A( 2, 1 ), H3 )
00189 *
00190             IF( RR.GT.QQ ) THEN
00191 *
00192 *              find right rotation matrix to zero 1,1 element of
00193 *              (sA - wB)
00194 *
00195                CALL DLARTG( H2, H1, CSR, SNR, T )
00196 *
00197             ELSE
00198 *
00199 *              find right rotation matrix to zero 2,1 element of
00200 *              (sA - wB)
00201 *
00202                CALL DLARTG( H3, SCALE1*A( 2, 1 ), CSR, SNR, T )
00203 *
00204             END IF
00205 *
00206             SNR = -SNR
00207             CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR )
00208             CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR )
00209 *
00210 *           compute inf norms of A and B
00211 *
00212             H1 = MAX( ABS( A( 1, 1 ) )+ABS( A( 1, 2 ) ),
00213      $           ABS( A( 2, 1 ) )+ABS( A( 2, 2 ) ) )
00214             H2 = MAX( ABS( B( 1, 1 ) )+ABS( B( 1, 2 ) ),
00215      $           ABS( B( 2, 1 ) )+ABS( B( 2, 2 ) ) )
00216 *
00217             IF( ( SCALE1*H1 ).GE.ABS( WR1 )*H2 ) THEN
00218 *
00219 *              find left rotation matrix Q to zero out B(2,1)
00220 *
00221                CALL DLARTG( B( 1, 1 ), B( 2, 1 ), CSL, SNL, R )
00222 *
00223             ELSE
00224 *
00225 *              find left rotation matrix Q to zero out A(2,1)
00226 *
00227                CALL DLARTG( A( 1, 1 ), A( 2, 1 ), CSL, SNL, R )
00228 *
00229             END IF
00230 *
00231             CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL )
00232             CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL )
00233 *
00234             A( 2, 1 ) = ZERO
00235             B( 2, 1 ) = ZERO
00236 *
00237          ELSE
00238 *
00239 *           a pair of complex conjugate eigenvalues
00240 *           first compute the SVD of the matrix B
00241 *
00242             CALL DLASV2( B( 1, 1 ), B( 1, 2 ), B( 2, 2 ), R, T, SNR,
00243      $                   CSR, SNL, CSL )
00244 *
00245 *           Form (A,B) := Q(A,B)Z' where Q is left rotation matrix and
00246 *           Z is right rotation matrix computed from DLASV2
00247 *
00248             CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL )
00249             CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL )
00250             CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR )
00251             CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR )
00252 *
00253             B( 2, 1 ) = ZERO
00254             B( 1, 2 ) = ZERO
00255 *
00256          END IF
00257 *
00258       END IF
00259 *
00260 *     Unscaling
00261 *
00262       A( 1, 1 ) = ANORM*A( 1, 1 )
00263       A( 2, 1 ) = ANORM*A( 2, 1 )
00264       A( 1, 2 ) = ANORM*A( 1, 2 )
00265       A( 2, 2 ) = ANORM*A( 2, 2 )
00266       B( 1, 1 ) = BNORM*B( 1, 1 )
00267       B( 2, 1 ) = BNORM*B( 2, 1 )
00268       B( 1, 2 ) = BNORM*B( 1, 2 )
00269       B( 2, 2 ) = BNORM*B( 2, 2 )
00270 *
00271       IF( WI.EQ.ZERO ) THEN
00272          ALPHAR( 1 ) = A( 1, 1 )
00273          ALPHAR( 2 ) = A( 2, 2 )
00274          ALPHAI( 1 ) = ZERO
00275          ALPHAI( 2 ) = ZERO
00276          BETA( 1 ) = B( 1, 1 )
00277          BETA( 2 ) = B( 2, 2 )
00278       ELSE
00279          ALPHAR( 1 ) = ANORM*WR1 / SCALE1 / BNORM
00280          ALPHAI( 1 ) = ANORM*WI / SCALE1 / BNORM
00281          ALPHAR( 2 ) = ALPHAR( 1 )
00282          ALPHAI( 2 ) = -ALPHAI( 1 )
00283          BETA( 1 ) = ONE
00284          BETA( 2 ) = ONE
00285       END IF
00286 *
00287       RETURN
00288 *
00289 *     End of DLAGV2
00290 *
00291       END