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: *> \ingroup doubleOTHERauxiliary
148: *
149: *> \par Contributors:
150: * ==================
151: *>
152: *> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
153: *
154: * =====================================================================
155: SUBROUTINE DLAGV2( A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL,
156: $ CSR, SNR )
157: *
158: * -- LAPACK auxiliary routine --
159: * -- LAPACK is a software package provided by Univ. of Tennessee, --
160: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
161: *
162: * .. Scalar Arguments ..
163: INTEGER LDA, LDB
164: DOUBLE PRECISION CSL, CSR, SNL, SNR
165: * ..
166: * .. Array Arguments ..
167: DOUBLE PRECISION A( LDA, * ), ALPHAI( 2 ), ALPHAR( 2 ),
168: $ B( LDB, * ), BETA( 2 )
169: * ..
170: *
171: * =====================================================================
172: *
173: * .. Parameters ..
174: DOUBLE PRECISION ZERO, ONE
175: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
176: * ..
177: * .. Local Scalars ..
178: DOUBLE PRECISION ANORM, ASCALE, BNORM, BSCALE, H1, H2, H3, QQ,
179: $ R, RR, SAFMIN, SCALE1, SCALE2, T, ULP, WI, WR1,
180: $ WR2
181: * ..
182: * .. External Subroutines ..
183: EXTERNAL DLAG2, DLARTG, DLASV2, DROT
184: * ..
185: * .. External Functions ..
186: DOUBLE PRECISION DLAMCH, DLAPY2
187: EXTERNAL DLAMCH, DLAPY2
188: * ..
189: * .. Intrinsic Functions ..
190: INTRINSIC ABS, MAX
191: * ..
192: * .. Executable Statements ..
193: *
194: SAFMIN = DLAMCH( 'S' )
195: ULP = DLAMCH( 'P' )
196: *
197: * Scale A
198: *
199: ANORM = MAX( ABS( A( 1, 1 ) )+ABS( A( 2, 1 ) ),
200: $ ABS( A( 1, 2 ) )+ABS( A( 2, 2 ) ), SAFMIN )
201: ASCALE = ONE / ANORM
202: A( 1, 1 ) = ASCALE*A( 1, 1 )
203: A( 1, 2 ) = ASCALE*A( 1, 2 )
204: A( 2, 1 ) = ASCALE*A( 2, 1 )
205: A( 2, 2 ) = ASCALE*A( 2, 2 )
206: *
207: * Scale B
208: *
209: BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 1, 2 ) )+ABS( B( 2, 2 ) ),
210: $ SAFMIN )
211: BSCALE = ONE / BNORM
212: B( 1, 1 ) = BSCALE*B( 1, 1 )
213: B( 1, 2 ) = BSCALE*B( 1, 2 )
214: B( 2, 2 ) = BSCALE*B( 2, 2 )
215: *
216: * Check if A can be deflated
217: *
218: IF( ABS( A( 2, 1 ) ).LE.ULP ) THEN
219: CSL = ONE
220: SNL = ZERO
221: CSR = ONE
222: SNR = ZERO
223: A( 2, 1 ) = ZERO
224: B( 2, 1 ) = ZERO
225: WI = ZERO
226: *
227: * Check if B is singular
228: *
229: ELSE IF( ABS( B( 1, 1 ) ).LE.ULP ) THEN
230: CALL DLARTG( A( 1, 1 ), A( 2, 1 ), CSL, SNL, R )
231: CSR = ONE
232: SNR = ZERO
233: CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL )
234: CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL )
235: A( 2, 1 ) = ZERO
236: B( 1, 1 ) = ZERO
237: B( 2, 1 ) = ZERO
238: WI = ZERO
239: *
240: ELSE IF( ABS( B( 2, 2 ) ).LE.ULP ) THEN
241: CALL DLARTG( A( 2, 2 ), A( 2, 1 ), CSR, SNR, T )
242: SNR = -SNR
243: CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR )
244: CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR )
245: CSL = ONE
246: SNL = ZERO
247: A( 2, 1 ) = ZERO
248: B( 2, 1 ) = ZERO
249: B( 2, 2 ) = ZERO
250: WI = ZERO
251: *
252: ELSE
253: *
254: * B is nonsingular, first compute the eigenvalues of (A,B)
255: *
256: CALL DLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2,
257: $ WI )
258: *
259: IF( WI.EQ.ZERO ) THEN
260: *
261: * two real eigenvalues, compute s*A-w*B
262: *
263: H1 = SCALE1*A( 1, 1 ) - WR1*B( 1, 1 )
264: H2 = SCALE1*A( 1, 2 ) - WR1*B( 1, 2 )
265: H3 = SCALE1*A( 2, 2 ) - WR1*B( 2, 2 )
266: *
267: RR = DLAPY2( H1, H2 )
268: QQ = DLAPY2( SCALE1*A( 2, 1 ), H3 )
269: *
270: IF( RR.GT.QQ ) THEN
271: *
272: * find right rotation matrix to zero 1,1 element of
273: * (sA - wB)
274: *
275: CALL DLARTG( H2, H1, CSR, SNR, T )
276: *
277: ELSE
278: *
279: * find right rotation matrix to zero 2,1 element of
280: * (sA - wB)
281: *
282: CALL DLARTG( H3, SCALE1*A( 2, 1 ), CSR, SNR, T )
283: *
284: END IF
285: *
286: SNR = -SNR
287: CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR )
288: CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR )
289: *
290: * compute inf norms of A and B
291: *
292: H1 = MAX( ABS( A( 1, 1 ) )+ABS( A( 1, 2 ) ),
293: $ ABS( A( 2, 1 ) )+ABS( A( 2, 2 ) ) )
294: H2 = MAX( ABS( B( 1, 1 ) )+ABS( B( 1, 2 ) ),
295: $ ABS( B( 2, 1 ) )+ABS( B( 2, 2 ) ) )
296: *
297: IF( ( SCALE1*H1 ).GE.ABS( WR1 )*H2 ) THEN
298: *
299: * find left rotation matrix Q to zero out B(2,1)
300: *
301: CALL DLARTG( B( 1, 1 ), B( 2, 1 ), CSL, SNL, R )
302: *
303: ELSE
304: *
305: * find left rotation matrix Q to zero out A(2,1)
306: *
307: CALL DLARTG( A( 1, 1 ), A( 2, 1 ), CSL, SNL, R )
308: *
309: END IF
310: *
311: CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL )
312: CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL )
313: *
314: A( 2, 1 ) = ZERO
315: B( 2, 1 ) = ZERO
316: *
317: ELSE
318: *
319: * a pair of complex conjugate eigenvalues
320: * first compute the SVD of the matrix B
321: *
322: CALL DLASV2( B( 1, 1 ), B( 1, 2 ), B( 2, 2 ), R, T, SNR,
323: $ CSR, SNL, CSL )
324: *
325: * Form (A,B) := Q(A,B)Z**T where Q is left rotation matrix and
326: * Z is right rotation matrix computed from DLASV2
327: *
328: CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL )
329: CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL )
330: CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR )
331: CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR )
332: *
333: B( 2, 1 ) = ZERO
334: B( 1, 2 ) = ZERO
335: *
336: END IF
337: *
338: END IF
339: *
340: * Unscaling
341: *
342: A( 1, 1 ) = ANORM*A( 1, 1 )
343: A( 2, 1 ) = ANORM*A( 2, 1 )
344: A( 1, 2 ) = ANORM*A( 1, 2 )
345: A( 2, 2 ) = ANORM*A( 2, 2 )
346: B( 1, 1 ) = BNORM*B( 1, 1 )
347: B( 2, 1 ) = BNORM*B( 2, 1 )
348: B( 1, 2 ) = BNORM*B( 1, 2 )
349: B( 2, 2 ) = BNORM*B( 2, 2 )
350: *
351: IF( WI.EQ.ZERO ) THEN
352: ALPHAR( 1 ) = A( 1, 1 )
353: ALPHAR( 2 ) = A( 2, 2 )
354: ALPHAI( 1 ) = ZERO
355: ALPHAI( 2 ) = ZERO
356: BETA( 1 ) = B( 1, 1 )
357: BETA( 2 ) = B( 2, 2 )
358: ELSE
359: ALPHAR( 1 ) = ANORM*WR1 / SCALE1 / BNORM
360: ALPHAI( 1 ) = ANORM*WI / SCALE1 / BNORM
361: ALPHAR( 2 ) = ALPHAR( 1 )
362: ALPHAI( 2 ) = -ALPHAI( 1 )
363: BETA( 1 ) = ONE
364: BETA( 2 ) = ONE
365: END IF
366: *
367: RETURN
368: *
369: * End of DLAGV2
370: *
371: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlagv2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack