dtgsen(3)

NAME

DTGSEN - reorder the generalized real Schur decomposition

of a real matrix pair (A, B) (in terms of an orthonormal equiva

lence trans- formation Q' * (A, B) * Z), so that a selected clus

ter of eigenvalues appears in the leading diagonal blocks of the

upper quasi-triangular matrix A and the upper triangular B

SYNOPSIS

SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A,  LDA,
B,  LDB,  ALPHAR,  ALPHAI,  BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF,
WORK, LWORK, IWORK, LIWORK, INFO )
    LOGICAL        WANTQ, WANTZ
    INTEGER        IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK,
LWORK, M, N
    DOUBLE         PRECISION PL, PR
    LOGICAL        SELECT( * )
    INTEGER        IWORK( * )
    DOUBLE         PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ), B( LDB, * ), BETA( * ), DIF( * ), Q( LDQ, * ), WORK( *
), Z( LDZ, * )

PURPOSE

DTGSEN reorders the generalized real Schur decomposition

of a real matrix pair (A, B) (in terms of an orthonormal equiva

lence trans- formation Q' * (A, B) * Z), so that a selected clus

ter of eigenvalues appears in the leading diagonal blocks of the

upper quasi-triangular matrix A and the upper triangular B. The

leading columns of Q and Z form orthonormal bases of the corre

sponding left and right eigen- spaces (deflating subspaces). (A,

B) must be in generalized real Schur canonical form (as returned

by DGGES), i.e. A is block upper triangular with 1-by-1 and

2-by-2 diagonal blocks. B is upper triangular.

DTGSEN also computes the generalized eigenvalues

w(j) = (ALPHAR(j) + i*ALPHAI(j))/BETA(j)

of the reordered matrix pair (A, B).

Optionally, DTGSEN computes the estimates of reciprocal

condition numbers for eigenvalues and eigenspaces. These are Di

fu[(A11,B11), (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the

separation(s) between the matrix pairs (A11, B11) and (A22,B22)

that correspond to the selected cluster and the eigenvalues out

side the cluster, resp., and norms of "projections" onto left and

right eigenspaces w.r.t. the selected cluster in the

(1,1)-block.

ARGUMENTS

IJOB (input) INTEGER

Specifies whether condition numbers are required

for the cluster of eigenvalues (PL and PR) or the deflating sub

spaces (Difu and Difl):
=0: Only reorder w.r.t. SELECT. No extras.
=1: Reciprocal of norms of "projections" onto left

and right eigenspaces w.r.t. the selected cluster (PL and PR).

=2: Upper bounds on Difu and Difl. F-norm-based estimate
(DIF(1:2)).
=3: Estimate of Difu and Difl. 1-norm-based esti

mate
(DIF(1:2)). About 5 times as expensive as IJOB =

2. =4: Compute PL, PR and DIF (i.e. 0, 1 and 2 above): Economic

version to get it all. =5: Compute PL, PR and DIF (i.e. 0, 1 and

3 above)

WANTQ (input) LOGICAL

WANTZ (input) LOGICAL

SELECT (input) LOGICAL array, dimension (N)

SELECT specifies the eigenvalues in the selected

cluster. To select a real eigenvalue w(j), SELECT(j) must be set

to w(j) and w(j+1), corresponding to a 2-by-2 diagonal block, ei

ther SELECT(j) or SELECT(j+1) or both must be set to either both

included in the cluster or both excluded.

N (input) INTEGER

The order of the matrices A and B. N >= 0.

A (input/output) DOUBLE PRECISION array, dimen

sion(LDA,N)

On entry, the upper quasi-triangular matrix A,

with (A, B) in generalized real Schur canonical form. On exit, A

is overwritten by the reordered matrix A.

LDA (input) INTEGER

The leading dimension of the array A. LDA >=

max(1,N).

B (input/output) DOUBLE PRECISION array, dimen

sion(LDB,N)

On entry, the upper triangular matrix B, with (A,

B) in generalized real Schur canonical form. On exit, B is over

written by the reordered matrix B.

LDB (input) INTEGER

The leading dimension of the array B. LDB >=

max(1,N).

ALPHAR (output) DOUBLE PRECISION array, dimension (N)

ALPHAI (output) DOUBLE PRECISION array, dimension

(N) BETA (output) DOUBLE PRECISION array, dimension (N) On ex

it, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will be the

generalized eigenvalues. ALPHAR(j) + ALPHAI(j)*i and BE

TA(j),j=1,...,N are the diagonals of the complex Schur form

(S,T) that would result if the 2-by-2 diagonal blocks of the real

generalized Schur form of (A,B) were further reduced to triangu

lar form using complex unitary transformations. If ALPHAI(j) is

zero, then the j-th eigenvalue is real; if positive, then the j

th and (j+1)-st eigenvalues are a complex conjugate pair, with

ALPHAI(j+1) negative.

Q (input/output) DOUBLE PRECISION array, dimension

(LDQ,N)

On entry, if WANTQ = .TRUE., Q is an N-by-N ma

trix. On exit, Q has been postmultiplied by the left orthogonal

transformation matrix which reorder (A, B); The leading M columns

of Q form orthonormal bases for the specified pair of left

eigenspaces (deflating subspaces). If WANTQ = .FALSE., Q is not

referenced.

LDQ (input) INTEGER

The leading dimension of the array Q. LDQ >= 1;

and if WANTQ = .TRUE., LDQ >= N.

Z (input/output) DOUBLE PRECISION array, dimension

(LDZ,N)

On entry, if WANTZ = .TRUE., Z is an N-by-N ma

trix. On exit, Z has been postmultiplied by the left orthogonal

transformation matrix which reorder (A, B); The leading M columns

of Z form orthonormal bases for the specified pair of left

eigenspaces (deflating subspaces). If WANTZ = .FALSE., Z is not

referenced.

LDZ (input) INTEGER

The leading dimension of the array Z. LDZ >= 1; If

WANTZ = .TRUE., LDZ >= N.

M (output) INTEGER

The dimension of the specified pair of left and

right eigen- spaces (deflating subspaces). 0 <= M <= N.

PL, PR (output) DOUBLE PRECISION If IJOB = 1, 4

or 5, PL, PR are lower bounds on the reciprocal of the norm of

"projections" onto left and right eigenspaces with respect to the

selected cluster. 0 < PL, PR <= 1. If M = 0 or M = N, PL = PR

= 1. If IJOB = 0, 2 or 3, PL and PR are not referenced.

DIF (output) DOUBLE PRECISION array, dimension (2).

If IJOB >= 2, DIF(1:2) store the estimates of Difu

and Difl.
If IJOB = 2 or 4, DIF(1:2) are F-norm-based upper

bounds on
Difu and Difl. If IJOB = 3 or 5, DIF(1:2) are

1-norm-based estimates of Difu and Difl. If M = 0 or N, DIF(1:2)

= F-norm([A, B]). If IJOB = 0 or 1, DIF is not referenced.

WORK (workspace/output) DOUBLE PRECISION array, dimen

sion (LWORK)

IF IJOB = 0, WORK is not referenced. Otherwise,

on exit, if INFO = 0, WORK(1) returns the optimal LWORK.

LWORK (input) INTEGER

The dimension of the array WORK. LWORK >= 4*N+16.

If IJOB = 1, 2 or 4, LWORK >= MAX(4*N+16, 2*M*(N-M)). If IJOB =

3 or 5, LWORK >= MAX(4*N+16, 4*M*(N-M)).

If LWORK = -1, then a workspace query is assumed;

the routine only calculates the optimal size of the WORK array,

returns this value as the first entry of the WORK array, and no

error message related to LWORK is issued by XERBLA.

IWORK (workspace/output) INTEGER array, dimension (LI

WORK)

IF IJOB = 0, IWORK is not referenced. Otherwise,

on exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.

LIWORK (input) INTEGER

The dimension of the array IWORK. LIWORK >= 1. If

IJOB = 1, 2 or 4, LIWORK >= N+6. If IJOB = 3 or 5, LIWORK >=

MAX(2*M*(N-M), N+6).

If LIWORK = -1, then a workspace query is assumed;

the routine only calculates the optimal size of the IWORK array,

returns this value as the first entry of the IWORK array, and no

error message related to LIWORK is issued by XERBLA.

INFO (output) INTEGER

=0: Successful exit.
<0: If INFO = -i, the i-th argument had an illegal

value.
=1: Reordering of (A, B) failed because the trans

formed matrix pair (A, B) would be too far from generalized Schur

form; the problem is very ill-conditioned. (A, B) may have been

partially reordered. If requested, 0 is returned in DIF(*), PL

and PR.

FURTHER DETAILS

DTGSEN first collects the selected eigenvalues by comput

ing orthogonal U and W that move them to the top left corner of

(A, B). In other words, the selected eigenvalues are the eigen

values of (A11, B11) in:

U'*(A, B)*W = (A11 A12) (B11 B12) n1

( 0 A22),( 0 B22) n2

n1 n2 n1 n2

where N = n1+n2 and U' means the transpose of U. The first

n1 columns of U and W span the specified pair of left and right

eigenspaces (deflating subspaces) of (A, B).

If (A, B) has been obtained from the generalized real

Schur decomposition of a matrix pair (C, D) = Q*(A, B)*Z', then

the reordered generalized real Schur form of (C, D) is given by

(C, D) = (Q*U)*(U'*(A, B)*W)*(Z*W)',

and the first n1 columns of Q*U and Z*W span the corre

sponding deflating subspaces of (C, D) (Q and Z store Q*U and

Z*W, resp.).

Note that if the selected eigenvalue is sufficiently ill

conditioned, then its value may differ significantly from its

value before reordering.

The reciprocal condition numbers of the left and right

eigenspaces spanned by the first n1 columns of U and W (or Q*U

and Z*W) may be returned in DIF(1:2), corresponding to Difu and

Difl, resp.

The Difu and Difl are defined as:

Difu[(A11, B11), (A22, B22)] = sigma-min( Zu )

and

Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11,

B11)],

where sigma-min(Zu) is the smallest singular value of the

(2*n1*n2)-by-(2*n1*n2) matrix

Zu = [ kron(In2, A11) -kron(A22', In1) ]

[ kron(In2, B11) -kron(B22', In1) ].

Here, Inx is the identity matrix of size nx and A22' is

the transpose of A22. kron(X, Y) is the Kronecker product between

the matrices X and Y.

When DIF(2) is small, small changes in (A, B) can cause

large changes in the deflating subspace. An approximate (asymp

totic) bound on the maximum angular error in the computed deflat

ing subspaces is

EPS * norm((A, B)) / DIF(2),

where EPS is the machine precision.

The reciprocal norm of the projectors on the left and

right eigenspaces associated with (A11, B11) may be returned in

PL and PR. They are computed as follows. First we compute L and

R so that P*(A, B)*Q is block diagonal, where

P = ( I -L ) n1 Q = ( I R ) n1

( 0 I ) n2 and ( 0 I ) n2

n1 n2 n1 n2

and (L, R) is the solution to the generalized Sylvester

equation

A11*R - L*A22 = -A12
B11*R - L*B22 = -B12

Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F

norm(R)**2+1)**(-1/2). An approximate (asymptotic) bound on the

average absolute error of the selected eigenvalues is

EPS * norm((A, B)) / PL.

There are also global error bounds which valid for pertur

bations up to a certain restriction: A lower bound (x) on the

smallest F-norm(E,F) for which an eigenvalue of (A11, B11) may

move and coalesce with an eigenvalue of (A22, B22) under pertur

bation (E,F), (i.e. (A + E, B + F), is

x = min(Difu,Di

fl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)).

An approximate bound on x can be computed from DIF(1:2),

PL and PR.

If y = ( F-norm(E,F) / x) <= 1, the angles between the

perturbed (L', R') and unperturbed (L, R) left and right deflat

ing subspaces associated with the selected cluster in the

(1,1)-blocks can be bounded as

max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL *

PL)**(1/2))
max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR *

PR)**(1/2))

See LAPACK User's Guide section 4.11 or the following ref

erences for more information.

Note that if the default method for computing the Frobe

nius-norm- based estimate DIF is not wanted (see DLATDF), then

the parameter IDIFJB (see below) should be changed from 3 to 4

(routine DLATDF (IJOB = 2 will be used)). See DTGSYL for more de

tails.

Based on contributions by

Bo Kagstrom and Peter Poromaa, Department of Computing

Science,
Umea University, S-901 87 Umea, Sweden.

References
==========

[1] B. Kagstrom; A Direct Method for Reordering Eigenval

ues in the

Generalized Real Schur Form of a Regular Matrix Pair

(A, B), in
M.S. Moonen et al (eds), Linear Algebra for Large

Scale and
Real-Time Applications, Kluwer Academic Publ. 1993, pp

195-218.

[2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with

Specified

Eigenvalues of a Regular Matrix Pair (A, B) and Condi

tion
Estimation: Theory, Algorithms and Software,
Report UMINF - 94.04, Department of Computing Science,

Umea
University, S-901 87 Umea, Sweden, 1994. Also as LA

PACK Working
Note 87. To appear in Numerical Algorithms, 1996.

[3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms

and Software

for Solving the Generalized Sylvester Equation and Es

timating the
Separation between Regular Matrix Pairs, Report UMINF

- 93.23,
Department of Computing Science, Umea University,

S-901 87 Umea,
Sweden, December 1993, Revised April 1994, Also as LA

PACK Working
Note 75. To appear in ACM Trans. on Math. Software,

Vol 22, No 1,
1996.

LAPACK version 3.0 15 June 2000

DTGSEN(3)