LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlagv2.f
Go to the documentation of this file.
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 lagv2
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
dlag2
subroutine dlag2(a, lda, b, ldb, safmin, scale1, scale2, wr1, wr2, wi)
DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
Definition dlag2.f:156
dlagv2
subroutine dlagv2(a, lda, b, ldb, alphar, alphai, beta, csl, snl, csr, snr)
DLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A,...
Definition dlagv2.f:157
dlartg
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
dlasv2
subroutine dlasv2(f, g, h, ssmin, ssmax, snr, csr, snl, csl)
DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
Definition dlasv2.f:136
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92