![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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