LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \date September 2012
148 *
149 *> \ingroup doubleOTHERauxiliary
150 *
151 *> \par Contributors:
152 * ==================
153 *>
154 *> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
155 *
156 * =====================================================================
157  SUBROUTINE dlagv2( A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL,
158  $ csr, snr )
159 *
160 * -- LAPACK auxiliary routine (version 3.4.2) --
161 * -- LAPACK is a software package provided by Univ. of Tennessee, --
162 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
163 * September 2012
164 *
165 * .. Scalar Arguments ..
166  INTEGER lda, ldb
167  DOUBLE PRECISION csl, csr, snl, snr
168 * ..
169 * .. Array Arguments ..
170  DOUBLE PRECISION a( lda, * ), alphai( 2 ), alphar( 2 ),
171  $ b( ldb, * ), beta( 2 )
172 * ..
173 *
174 * =====================================================================
175 *
176 * .. Parameters ..
177  DOUBLE PRECISION zero, one
178  parameter( zero = 0.0d+0, one = 1.0d+0 )
179 * ..
180 * .. Local Scalars ..
181  DOUBLE PRECISION anorm, ascale, bnorm, bscale, h1, h2, h3, qq,
182  $ r, rr, safmin, scale1, scale2, t, ulp, wi, wr1,
183  $ wr2
184 * ..
185 * .. External Subroutines ..
186  EXTERNAL dlag2, dlartg, dlasv2, drot
187 * ..
188 * .. External Functions ..
189  DOUBLE PRECISION dlamch, dlapy2
190  EXTERNAL dlamch, dlapy2
191 * ..
192 * .. Intrinsic Functions ..
193  INTRINSIC abs, max
194 * ..
195 * .. Executable Statements ..
196 *
197  safmin = dlamch( 'S' )
198  ulp = dlamch( 'P' )
199 *
200 * Scale A
201 *
202  anorm = max( abs( a( 1, 1 ) )+abs( a( 2, 1 ) ),
203  $ abs( a( 1, 2 ) )+abs( a( 2, 2 ) ), safmin )
204  ascale = one / anorm
205  a( 1, 1 ) = ascale*a( 1, 1 )
206  a( 1, 2 ) = ascale*a( 1, 2 )
207  a( 2, 1 ) = ascale*a( 2, 1 )
208  a( 2, 2 ) = ascale*a( 2, 2 )
209 *
210 * Scale B
211 *
212  bnorm = max( abs( b( 1, 1 ) ), abs( b( 1, 2 ) )+abs( b( 2, 2 ) ),
213  $ safmin )
214  bscale = one / bnorm
215  b( 1, 1 ) = bscale*b( 1, 1 )
216  b( 1, 2 ) = bscale*b( 1, 2 )
217  b( 2, 2 ) = bscale*b( 2, 2 )
218 *
219 * Check if A can be deflated
220 *
221  IF( abs( a( 2, 1 ) ).LE.ulp ) THEN
222  csl = one
223  snl = zero
224  csr = one
225  snr = zero
226  a( 2, 1 ) = zero
227  b( 2, 1 ) = zero
228  wi = zero
229 *
230 * Check if B is singular
231 *
232  ELSE IF( abs( b( 1, 1 ) ).LE.ulp ) THEN
233  CALL dlartg( a( 1, 1 ), a( 2, 1 ), csl, snl, r )
234  csr = one
235  snr = zero
236  CALL drot( 2, a( 1, 1 ), lda, a( 2, 1 ), lda, csl, snl )
237  CALL drot( 2, b( 1, 1 ), ldb, b( 2, 1 ), ldb, csl, snl )
238  a( 2, 1 ) = zero
239  b( 1, 1 ) = zero
240  b( 2, 1 ) = zero
241  wi = zero
242 *
243  ELSE IF( abs( b( 2, 2 ) ).LE.ulp ) THEN
244  CALL dlartg( a( 2, 2 ), a( 2, 1 ), csr, snr, t )
245  snr = -snr
246  CALL drot( 2, a( 1, 1 ), 1, a( 1, 2 ), 1, csr, snr )
247  CALL drot( 2, b( 1, 1 ), 1, b( 1, 2 ), 1, csr, snr )
248  csl = one
249  snl = zero
250  a( 2, 1 ) = zero
251  b( 2, 1 ) = zero
252  b( 2, 2 ) = zero
253  wi = zero
254 *
255  ELSE
256 *
257 * B is nonsingular, first compute the eigenvalues of (A,B)
258 *
259  CALL dlag2( a, lda, b, ldb, safmin, scale1, scale2, wr1, wr2,
260  $ wi )
261 *
262  IF( wi.EQ.zero ) THEN
263 *
264 * two real eigenvalues, compute s*A-w*B
265 *
266  h1 = scale1*a( 1, 1 ) - wr1*b( 1, 1 )
267  h2 = scale1*a( 1, 2 ) - wr1*b( 1, 2 )
268  h3 = scale1*a( 2, 2 ) - wr1*b( 2, 2 )
269 *
270  rr = dlapy2( h1, h2 )
271  qq = dlapy2( scale1*a( 2, 1 ), h3 )
272 *
273  IF( rr.GT.qq ) THEN
274 *
275 * find right rotation matrix to zero 1,1 element of
276 * (sA - wB)
277 *
278  CALL dlartg( h2, h1, csr, snr, t )
279 *
280  ELSE
281 *
282 * find right rotation matrix to zero 2,1 element of
283 * (sA - wB)
284 *
285  CALL dlartg( h3, scale1*a( 2, 1 ), csr, snr, t )
286 *
287  END IF
288 *
289  snr = -snr
290  CALL drot( 2, a( 1, 1 ), 1, a( 1, 2 ), 1, csr, snr )
291  CALL drot( 2, b( 1, 1 ), 1, b( 1, 2 ), 1, csr, snr )
292 *
293 * compute inf norms of A and B
294 *
295  h1 = max( abs( a( 1, 1 ) )+abs( a( 1, 2 ) ),
296  $ abs( a( 2, 1 ) )+abs( a( 2, 2 ) ) )
297  h2 = max( abs( b( 1, 1 ) )+abs( b( 1, 2 ) ),
298  $ abs( b( 2, 1 ) )+abs( b( 2, 2 ) ) )
299 *
300  IF( ( scale1*h1 ).GE.abs( wr1 )*h2 ) THEN
301 *
302 * find left rotation matrix Q to zero out B(2,1)
303 *
304  CALL dlartg( b( 1, 1 ), b( 2, 1 ), csl, snl, r )
305 *
306  ELSE
307 *
308 * find left rotation matrix Q to zero out A(2,1)
309 *
310  CALL dlartg( a( 1, 1 ), a( 2, 1 ), csl, snl, r )
311 *
312  END IF
313 *
314  CALL drot( 2, a( 1, 1 ), lda, a( 2, 1 ), lda, csl, snl )
315  CALL drot( 2, b( 1, 1 ), ldb, b( 2, 1 ), ldb, csl, snl )
316 *
317  a( 2, 1 ) = zero
318  b( 2, 1 ) = zero
319 *
320  ELSE
321 *
322 * a pair of complex conjugate eigenvalues
323 * first compute the SVD of the matrix B
324 *
325  CALL dlasv2( b( 1, 1 ), b( 1, 2 ), b( 2, 2 ), r, t, snr,
326  $ csr, snl, csl )
327 *
328 * Form (A,B) := Q(A,B)Z**T where Q is left rotation matrix and
329 * Z is right rotation matrix computed from DLASV2
330 *
331  CALL drot( 2, a( 1, 1 ), lda, a( 2, 1 ), lda, csl, snl )
332  CALL drot( 2, b( 1, 1 ), ldb, b( 2, 1 ), ldb, csl, snl )
333  CALL drot( 2, a( 1, 1 ), 1, a( 1, 2 ), 1, csr, snr )
334  CALL drot( 2, b( 1, 1 ), 1, b( 1, 2 ), 1, csr, snr )
335 *
336  b( 2, 1 ) = zero
337  b( 1, 2 ) = zero
338 *
339  END IF
340 *
341  END IF
342 *
343 * Unscaling
344 *
345  a( 1, 1 ) = anorm*a( 1, 1 )
346  a( 2, 1 ) = anorm*a( 2, 1 )
347  a( 1, 2 ) = anorm*a( 1, 2 )
348  a( 2, 2 ) = anorm*a( 2, 2 )
349  b( 1, 1 ) = bnorm*b( 1, 1 )
350  b( 2, 1 ) = bnorm*b( 2, 1 )
351  b( 1, 2 ) = bnorm*b( 1, 2 )
352  b( 2, 2 ) = bnorm*b( 2, 2 )
353 *
354  IF( wi.EQ.zero ) THEN
355  alphar( 1 ) = a( 1, 1 )
356  alphar( 2 ) = a( 2, 2 )
357  alphai( 1 ) = zero
358  alphai( 2 ) = zero
359  beta( 1 ) = b( 1, 1 )
360  beta( 2 ) = b( 2, 2 )
361  ELSE
362  alphar( 1 ) = anorm*wr1 / scale1 / bnorm
363  alphai( 1 ) = anorm*wi / scale1 / bnorm
364  alphar( 2 ) = alphar( 1 )
365  alphai( 2 ) = -alphai( 1 )
366  beta( 1 ) = one
367  beta( 2 ) = one
368  END IF
369 *
370  return
371 *
372 * End of DLAGV2
373 *
374  END