REDUCE CHARACTER*1
Indicates whether the matrix pairs (A,D) and/or (B,E) are
to be reduced to generalized Schur form as follows:
= 'R': The matrix pairs (A,D) and (B,E) are to be reduced
to generalized (real) Schur canonical form;
= 'A': The matrix pair (A,D) only is to be reduced
to generalized (real) Schur canonical form,
and the matrix pair (B,E) already is in this form;
= 'B': The matrix pair (B,E) only is to be reduced
to generalized (real) Schur canonical form,
and the matrix pair (A,D) already is in this form;
= 'N': The matrix pairs (A,D) and (B,E) are already in
generalized (real) Schur canonical form, as
produced by LAPACK routine DGGES.
TRANS CHARACTER*1
Indicates which of the equations, (1) or (2), is to be
solved as follows:
= 'N': The generalized Sylvester equation (1) is to be
solved;
= 'T': The "transposed" generalized Sylvester equation
(2) is to be solved.
JOBD CHARACTER*1
Indicates whether the Dif estimator is to be computed as
follows:
= '1': Only the one-norm-based Dif estimate is computed
and stored in DIF;
= '2': Only the Frobenius norm-based Dif estimate is
computed and stored in DIF;
= 'D': The equation (1) is solved and the one-norm-based
Dif estimate is computed and stored in DIF;
= 'F': The equation (1) is solved and the Frobenius norm-
based Dif estimate is computed and stored in DIF;
= 'N': The Dif estimator is not required and hence DIF is
not referenced. (Solve either (1) or (2) only.)
JOBD is not referenced if TRANS = 'T'.

Input/Output Parameters

M (input) INTEGER
The order of the matrices A and D and the number of rows
of the matrices C, F, R and L. M >= 0.
N (input) INTEGER
The order of the matrices B and E and the number of
columns of the matrices C, F, R and L. N >= 0.
No computations are performed if N = 0 or M = 0, but SCALE
and DIF (if JOB <> 'N') are set to 1.
A (input/output) DOUBLE PRECISION array, dimension (LDA,M)
On entry, the leading M-by-M part of this array must
contain the coefficient matrix A of the equation; A must
be in upper quasi-triangular form if REDUCE = 'B' or 'N'.
On exit, the leading M-by-M part of this array contains
the upper quasi-triangular form of A.
LDA INTEGER
The leading dimension of array A. LDA >= MAX(1,M).
B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
On entry, the leading N-by-N part of this array must
contain the coefficient matrix B of the equation; B must
be in upper quasi-triangular form if REDUCE = 'A' or 'N'.
On exit, the leading N-by-N part of this array contains
the upper quasi-triangular form of B.
LDB INTEGER
The leading dimension of array B. LDB >= MAX(1,N).
C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
On entry, the leading M-by-N part of this array must
contain the right-hand side matrix C of the first equation
in (1) or (2).
On exit, if JOBD = 'N', 'D' or 'F', the leading M-by-N
part of this array contains the solution matrix R of the
problem; if JOBD = '1' or '2' and TRANS = 'N', the leading
M-by-N part of this array contains the solution matrix R
achieved during the computation of the Dif estimate.
LDC INTEGER
The leading dimension of array C. LDC >= MAX(1,M).
D (input/output) DOUBLE PRECISION array, dimension (LDD,M)
On entry, the leading M-by-M part of this array must
contain the coefficient matrix D of the equation; D must
be in upper triangular form if REDUCE = 'B' or 'N'.
On exit, the leading M-by-M part of this array contains
the upper triangular form of D.
LDD INTEGER
The leading dimension of array D. LDD >= MAX(1,M).
E (input/output) DOUBLE PRECISION array, dimension (LDE,N)
On entry, the leading N-by-N part of this array must
contain the coefficient matrix E of the equation; E must
be in upper triangular form if REDUCE = 'A' or 'N'.
On exit, the leading N-by-N part of this array contains
the upper triangular form of E.
LDE INTEGER
The leading dimension of array E. LDE >= MAX(1,N).
F (input/output) DOUBLE PRECISION array, dimension (LDF,N)
On entry, the leading M-by-N part of this array must
contain the right-hand side matrix F of the second
equation in (1) or (2).
On exit, if JOBD = 'N', 'D' or 'F', the leading M-by-N
part of this array contains the solution matrix L of the
problem; if JOBD = '1' or '2' and TRANS = 'N', the leading
M-by-N part of this array contains the solution matrix L
achieved during the computation of the Dif estimate.
LDF INTEGER
The leading dimension of array F. LDF >= MAX(1,M).
SCALE (output) DOUBLE PRECISION
The scaling factor in (1) or (2). If 0 < SCALE < 1, C and
F hold the solutions R and L, respectively, to a slightly
perturbed system, but the computed generalized (real)
Schur canonical form matrices P'*A*Q, U'*B*V, P'*D*Q and
U'*E*V (see METHOD), or input matrices A, B, D, and E (if
already reduced to these forms), have not been changed.
If SCALE = 0, C and F hold the solutions R and L,
respectively, to the homogeneous system with C = F = 0.
Normally, SCALE = 1.
DIF (output) DOUBLE PRECISION
If TRANS = 'N' and JOBD <> 'N', then DIF contains the
value of the Dif estimator, which is an upper bound of
-1
Dif[(A,D),(B,E)] = sigma_min(Z) = 1/||Z ||, in either the
one-norm, or Frobenius norm, respectively (see METHOD).
Otherwise, DIF is not referenced.
P (output) DOUBLE PRECISION array, dimension (LDP,*)
If REDUCE = 'R' or 'A', then the leading M-by-M part of
this array contains the (left) transformation matrix used
to reduce (A,D) to generalized Schur form.
Otherwise, P is not referenced and can be supplied as a
dummy array (i.e. set parameter LDP = 1 and declare this
array to be P(1,1) in the calling program).
LDP INTEGER
The leading dimension of array P.
LDP >= MAX(1,M) if REDUCE = 'R' or 'A',
LDP >= 1 if REDUCE = 'B' or 'N'.
Q (output) DOUBLE PRECISION array, dimension (LDQ,*)
If REDUCE = 'R' or 'A', then the leading M-by-M part of
this array contains the (right) transformation matrix used
to reduce (A,D) to generalized Schur form.
Otherwise, Q is not referenced and can be supplied as a
dummy array (i.e. set parameter LDQ = 1 and declare this
array to be Q(1,1) in the calling program).
LDQ INTEGER
The leading dimension of array Q.
LDQ >= MAX(1,M) if REDUCE = 'R' or 'A',
LDQ >= 1 if REDUCE = 'B' or 'N'.
U (output) DOUBLE PRECISION array, dimension (LDU,*)
If REDUCE = 'R' or 'B', then the leading N-by-N part of
this array contains the (left) transformation matrix used
to reduce (B,E) to generalized Schur form.
Otherwise, U is not referenced and can be supplied as a
dummy array (i.e. set parameter LDU = 1 and declare this
array to be U(1,1) in the calling program).
LDU INTEGER
The leading dimension of array U.
LDU >= MAX(1,N) if REDUCE = 'R' or 'B',
LDU >= 1 if REDUCE = 'A' or 'N'.
V (output) DOUBLE PRECISION array, dimension (LDV,*)
If REDUCE = 'R' or 'B', then the leading N-by-N part of
this array contains the (right) transformation matrix used
to reduce (B,E) to generalized Schur form.
Otherwise, V is not referenced and can be supplied as a
dummy array (i.e. set parameter LDV = 1 and declare this
array to be V(1,1) in the calling program).
LDV INTEGER
The leading dimension of array V.
LDV >= MAX(1,N) if REDUCE = 'R' or 'B',
LDV >= 1 if REDUCE = 'A' or 'N'.

Workspace

IWORK INTEGER array, dimension (M+N+6)
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, DWORK(1) returns the optimal value
of LDWORK.
LDWORK INTEGER
The length of the array DWORK.
If TRANS = 'N' and JOBD = 'D' or 'F', then
LDWORK = MAX(1,11*MN,10*MN+23,2*M*N) if REDUCE = 'R';
LDWORK = MAX(1,11*M, 10*M+23, 2*M*N) if REDUCE = 'A';
LDWORK = MAX(1,11*N, 10*N+23, 2*M*N) if REDUCE = 'B';
LDWORK = MAX(1,2*M*N) if REDUCE = 'N',
where MN = max(M,N).
Otherwise, the term 2*M*N above should be omitted.
For optimum performance LDWORK should be larger.
If LDWORK = -1 a workspace query is assumed; the
routine only calculates the optimal size of the DWORK
array, returns this value as the first entry of the DWORK
array, and no error message is issued by XERBLA.

Error Indicator

INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value;
= 1: if REDUCE <> 'N' and either (A,D) and/or (B,E)
cannot be reduced to generalized Schur form;
= 2: if REDUCE = 'N' and either A or B is not in
upper quasi-triangular form;
= 3: if a singular matrix was encountered during the
computation of the solution matrices R and L, that
is (A,D) and (B,E) have common or close eigenvalues.

The algorithm is backward stable. A reliable estimate for the
condition number of Z in the Frobenius norm, is (see [1])
K(Z) = SQRT( ||A||**2 + ||B||**2 + ||C||**2 + ||D||**2 )/DIF.
If mu is an upper bound on the relative error of the elements of
the matrices A, B, C, D, E and F, then the relative error in the
actual solution is approximately mu * K(Z).
The relative error in the computed solution (due to rounding
errors) is approximately EPS * K(Z), where EPS is the machine
precision (see LAPACK Library routine DLAMCH).