LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dlagv2()

subroutine dlagv2 ( double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( 2 ) alphar,
double precision, dimension( 2 ) alphai,
double precision, dimension( 2 ) beta,
double precision csl,
double precision snl,
double precision csr,
double precision snr )

DLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A,B) where B is upper triangular.

Download DLAGV2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DLAGV2 computes the Generalized Schur factorization of a real 2-by-2
!> matrix pencil (A,B) where B is upper triangular. This routine
!> computes orthogonal (rotation) matrices given by CSL, SNL and CSR,
!> SNR such that
!>
!> 1) if the pencil (A,B) has two real eigenvalues (include 0/0 or 1/0
!>    types), then
!>
!>    [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]
!>    [  0  a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]
!>
!>    [ b11 b12 ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]
!>    [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ],
!>
!> 2) if the pencil (A,B) has a pair of complex conjugate eigenvalues,
!>    then
!>
!>    [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]
!>    [ a21 a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]
!>
!>    [ b11  0  ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]
!>    [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ]
!>
!>    where b11 >= b22 > 0.
!>
!>
Parameters
[in,out]A
!>          A is DOUBLE PRECISION array, dimension (LDA, 2)
!>          On entry, the 2 x 2 matrix A.
!>          On exit, A is overwritten by the ``A-part'' of the
!>          generalized Schur form.
!>
[in]LDA
!>          LDA is INTEGER
!>          THe leading dimension of the array A.  LDA >= 2.
!>
[in,out]B
!>          B is DOUBLE PRECISION array, dimension (LDB, 2)
!>          On entry, the upper triangular 2 x 2 matrix B.
!>          On exit, B is overwritten by the ``B-part'' of the
!>          generalized Schur form.
!>
[in]LDB
!>          LDB is INTEGER
!>          THe leading dimension of the array B.  LDB >= 2.
!>
[out]ALPHAR
!>          ALPHAR is DOUBLE PRECISION array, dimension (2)
!>
[out]ALPHAI
!>          ALPHAI is DOUBLE PRECISION array, dimension (2)
!>
[out]BETA
!>          BETA is DOUBLE PRECISION array, dimension (2)
!>          (ALPHAR(k)+i*ALPHAI(k))/BETA(k) are the eigenvalues of the
!>          pencil (A,B), k=1,2, i = sqrt(-1).  Note that BETA(k) may
!>          be zero.
!>
[out]CSL
!>          CSL is DOUBLE PRECISION
!>          The cosine of the left rotation matrix.
!>
[out]SNL
!>          SNL is DOUBLE PRECISION
!>          The sine of the left rotation matrix.
!>
[out]CSR
!>          CSR is DOUBLE PRECISION
!>          The cosine of the right rotation matrix.
!>
[out]SNR
!>          SNR is DOUBLE PRECISION
!>          The sine of the right rotation matrix.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA

Definition at line 153 of file dlagv2.f.

156*
157* -- LAPACK auxiliary routine --
158* -- LAPACK is a software package provided by Univ. of Tennessee, --
159* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
160*
161* .. Scalar Arguments ..
162 INTEGER LDA, LDB
163 DOUBLE PRECISION CSL, CSR, SNL, SNR
164* ..
165* .. Array Arguments ..
166 DOUBLE PRECISION A( LDA, * ), ALPHAI( 2 ), ALPHAR( 2 ),
167 $ B( LDB, * ), BETA( 2 )
168* ..
169*
170* =====================================================================
171*
172* .. Parameters ..
173 DOUBLE PRECISION ZERO, ONE
174 parameter( zero = 0.0d+0, one = 1.0d+0 )
175* ..
176* .. Local Scalars ..
177 DOUBLE PRECISION ANORM, ASCALE, BNORM, BSCALE, H1, H2, H3, QQ,
178 $ R, RR, SAFMIN, SCALE1, SCALE2, T, ULP, WI, WR1,
179 $ WR2
180* ..
181* .. External Subroutines ..
182 EXTERNAL dlag2, dlartg, dlasv2, drot
183* ..
184* .. External Functions ..
185 DOUBLE PRECISION DLAMCH, DLAPY2
186 EXTERNAL dlamch, dlapy2
187* ..
188* .. Intrinsic Functions ..
189 INTRINSIC abs, max
190* ..
191* .. Executable Statements ..
192*
193 safmin = dlamch( 'S' )
194 ulp = dlamch( 'P' )
195*
196* Scale A
197*
198 anorm = max( abs( a( 1, 1 ) )+abs( a( 2, 1 ) ),
199 $ abs( a( 1, 2 ) )+abs( a( 2, 2 ) ), safmin )
200 ascale = one / anorm
201 a( 1, 1 ) = ascale*a( 1, 1 )
202 a( 1, 2 ) = ascale*a( 1, 2 )
203 a( 2, 1 ) = ascale*a( 2, 1 )
204 a( 2, 2 ) = ascale*a( 2, 2 )
205*
206* Scale B
207*
208 bnorm = max( abs( b( 1, 1 ) ), abs( b( 1, 2 ) )+abs( b( 2, 2 ) ),
209 $ safmin )
210 bscale = one / bnorm
211 b( 1, 1 ) = bscale*b( 1, 1 )
212 b( 1, 2 ) = bscale*b( 1, 2 )
213 b( 2, 2 ) = bscale*b( 2, 2 )
214*
215* Check if A can be deflated
216*
217 IF( abs( a( 2, 1 ) ).LE.ulp ) THEN
218 csl = one
219 snl = zero
220 csr = one
221 snr = zero
222 a( 2, 1 ) = zero
223 b( 2, 1 ) = zero
224 wi = zero
225*
226* Check if B is singular
227*
228 ELSE IF( abs( b( 1, 1 ) ).LE.ulp ) THEN
229 CALL dlartg( a( 1, 1 ), a( 2, 1 ), csl, snl, r )
230 csr = one
231 snr = zero
232 CALL drot( 2, a( 1, 1 ), lda, a( 2, 1 ), lda, csl, snl )
233 CALL drot( 2, b( 1, 1 ), ldb, b( 2, 1 ), ldb, csl, snl )
234 a( 2, 1 ) = zero
235 b( 1, 1 ) = zero
236 b( 2, 1 ) = zero
237 wi = zero
238*
239 ELSE IF( abs( b( 2, 2 ) ).LE.ulp ) THEN
240 CALL dlartg( a( 2, 2 ), a( 2, 1 ), csr, snr, t )
241 snr = -snr
242 CALL drot( 2, a( 1, 1 ), 1, a( 1, 2 ), 1, csr, snr )
243 CALL drot( 2, b( 1, 1 ), 1, b( 1, 2 ), 1, csr, snr )
244 csl = one
245 snl = zero
246 a( 2, 1 ) = zero
247 b( 2, 1 ) = zero
248 b( 2, 2 ) = zero
249 wi = zero
250*
251 ELSE
252*
253* B is nonsingular, first compute the eigenvalues of (A,B)
254*
255 CALL dlag2( a, lda, b, ldb, safmin, scale1, scale2, wr1,
256 $ 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*
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:154
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlapy2
double precision function dlapy2(x, y)
DLAPY2 returns sqrt(x2+y2).
Definition dlapy2.f:61
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:134
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
Here is the call graph for this function:
Here is the caller graph for this function: