LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
slagv2.f
1 *> \brief \b SLAGV2 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 SLAGV2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slagv2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slagv2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slagv2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLAGV2( A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL,
22 * CSR, SNR )
23 *
24 * .. Scalar Arguments ..
25 * INTEGER LDA, LDB
26 * REAL CSL, CSR, SNL, SNR
27 * ..
28 * .. Array Arguments ..
29 * REAL A( LDA, * ), ALPHAI( 2 ), ALPHAR( 2 ),
30 * $ B( LDB, * ), BETA( 2 )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> SLAGV2 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 REAL 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 REAL 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 REAL array, dimension (2)
100 *> \endverbatim
101 *>
102 *> \param[out] ALPHAI
103 *> \verbatim
104 *> ALPHAI is REAL array, dimension (2)
105 *> \endverbatim
106 *>
107 *> \param[out] BETA
108 *> \verbatim
109 *> BETA is REAL 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 REAL
118 *> The cosine of the left rotation matrix.
119 *> \endverbatim
120 *>
121 *> \param[out] SNL
122 *> \verbatim
123 *> SNL is REAL
124 *> The sine of the left rotation matrix.
125 *> \endverbatim
126 *>
127 *> \param[out] CSR
128 *> \verbatim
129 *> CSR is REAL
130 *> The cosine of the right rotation matrix.
131 *> \endverbatim
132 *>
133 *> \param[out] SNR
134 *> \verbatim
135 *> SNR is REAL
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 slagv2( 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  REAL CSL, CSR, SNL, SNR
165 * ..
166 * .. Array Arguments ..
167  REAL A( lda, * ), ALPHAI( 2 ), ALPHAR( 2 ),
168  $ b( ldb, * ), beta( 2 )
169 * ..
170 *
171 * =====================================================================
172 *
173 * .. Parameters ..
174  REAL ZERO, ONE
175  parameter( zero = 0.0e+0, one = 1.0e+0 )
176 * ..
177 * .. Local Scalars ..
178  REAL 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 slag2, slartg, slasv2, srot
184 * ..
185 * .. External Functions ..
186  REAL SLAMCH, SLAPY2
187  EXTERNAL slamch, slapy2
188 * ..
189 * .. Intrinsic Functions ..
190  INTRINSIC abs, max
191 * ..
192 * .. Executable Statements ..
193 *
194  safmin = slamch( 'S' )
195  ulp = slamch( '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 slartg( a( 1, 1 ), a( 2, 1 ), csl, snl, r )
231  csr = one
232  snr = zero
233  CALL srot( 2, a( 1, 1 ), lda, a( 2, 1 ), lda, csl, snl )
234  CALL srot( 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 slartg( a( 2, 2 ), a( 2, 1 ), csr, snr, t )
242  snr = -snr
243  CALL srot( 2, a( 1, 1 ), 1, a( 1, 2 ), 1, csr, snr )
244  CALL srot( 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 slag2( 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 = slapy2( h1, h2 )
268  qq = slapy2( 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 slartg( 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 slartg( h3, scale1*a( 2, 1 ), csr, snr, t )
283 *
284  END IF
285 *
286  snr = -snr
287  CALL srot( 2, a( 1, 1 ), 1, a( 1, 2 ), 1, csr, snr )
288  CALL srot( 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 slartg( 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 slartg( a( 1, 1 ), a( 2, 1 ), csl, snl, r )
308 *
309  END IF
310 *
311  CALL srot( 2, a( 1, 1 ), lda, a( 2, 1 ), lda, csl, snl )
312  CALL srot( 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 slasv2( 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 SLASV2
327 *
328  CALL srot( 2, a( 1, 1 ), lda, a( 2, 1 ), lda, csl, snl )
329  CALL srot( 2, b( 1, 1 ), ldb, b( 2, 1 ), ldb, csl, snl )
330  CALL srot( 2, a( 1, 1 ), 1, a( 1, 2 ), 1, csr, snr )
331  CALL srot( 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 SLAGV2
370 *
371  END
slasv2
subroutine slasv2(F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL)
SLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
Definition: slasv2.f:138
slag2
subroutine slag2(A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2, WI)
SLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
Definition: slag2.f:156
srot
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:92
slagv2
subroutine slagv2(A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL, CSR, SNR)
SLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A...
Definition: slagv2.f:157