LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
September 2012
Contributors:
Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA

Definition at line 159 of file dlagv2.f.

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 *
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
drot
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:53
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:158
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:140
dlapy2
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
Definition: dlapy2.f:65
dlartg
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f:99

Here is the call graph for this function:

Here is the caller graph for this function: