Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_lapack_dtgsyl (f08yh)

## Purpose

nag_lapack_dtgsyl (f08yh) solves the generalized real quasi-triangular Sylvester equations.

## Syntax

[c, f, scale, dif, info] = f08yh(trans, ijob, a, b, c, d, e, f, 'm', m, 'n', n)
[c, f, scale, dif, info] = nag_lapack_dtgsyl(trans, ijob, a, b, c, d, e, f, 'm', m, 'n', n)

## Description

nag_lapack_dtgsyl (f08yh) solves either the generalized real Sylvester equations
 $AR-LB =αC DR-LE =αF,$ (1)
or the equations
 $ATR+DTL =αC RBT+LET =-αF,$ (2)
where the pair $\left(A,D\right)$ are given $m$ by $m$ matrices in real generalized Schur form, $\left(B,E\right)$ are given $n$ by $n$ matrices in real generalized Schur form and $\left(C,F\right)$ are given $m$ by $n$ matrices. The pair $\left(R,L\right)$ are the $m$ by $n$ solution matrices, and $\alpha$ is an output scaling factor determined by the function to avoid overflow in computing $\left(R,L\right)$.
Equations (1) are equivalent to equations of the form
 $Zx=αb ,$
where
 $Z = I⊗A-BT⊗I I⊗D-ET⊗I$
and $\otimes$ is the Kronecker product. Equations (2) are then equivalent to
 $ZTy = αb .$
The pair $\left(S,T\right)$ are in real generalized Schur form if $S$ is block upper triangular with $1$ by $1$ and $2$ by $2$ diagonal blocks on the diagonal and $T$ is upper triangular as returned, for example, by nag_lapack_dgges (f08xa), or nag_lapack_dhgeqz (f08xe) with ${\mathbf{job}}=\text{'S'}$.
Optionally, the function estimates $\mathrm{Dif}\left[\left(A,D\right),\left(B,E\right)\right]$, the separation between the matrix pairs $\left(A,D\right)$ and $\left(B,E\right)$, which is the smallest singular value of $Z$. The estimate can be based on either the Frobenius norm, or the $1$-norm. The $1$-norm estimate can be three to ten times more expensive than the Frobenius norm estimate, but makes the condition estimation uniform with the nonsymmetric eigenproblem. The Frobenius norm estimate provides a low cost, but equally reliable estimate. For more information see Sections 2.4.8.3 and 4.11.1.3 of Anderson et al. (1999) and Kågström and Poromaa (1996).

## References

Anderson E, Bai Z, Bischof C, Blackford S, Demmel J, Dongarra J J, Du Croz J J, Greenbaum A, Hammarling S, McKenney A and Sorensen D (1999) LAPACK Users' Guide (3rd Edition) SIAM, Philadelphia http://www.netlib.org/lapack/lug
Kågström B (1994) A perturbation analysis of the generalized Sylvester equation $\left(AR-LB,DR-LE\right)=\left(c,F\right)$ SIAM J. Matrix Anal. Appl. 15 1045–1060
Kågström B and Poromaa P (1996) LAPACK-style algorithms and software for solving the generalized Sylvester equation and estimating the separation between regular matrix pairs ACM Trans. Math. Software 22 78–103

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{trans}$ – string (length ≥ 1)
If ${\mathbf{trans}}=\text{'N'}$, solve the generalized Sylvester equation (1).
If ${\mathbf{trans}}=\text{'T'}$, solve the ‘transposed’ system (2).
Constraint: ${\mathbf{trans}}=\text{'N'}$ or $\text{'T'}$.
2:     $\mathrm{ijob}$int64int32nag_int scalar
Specifies what kind of functionality is to be performed when ${\mathbf{trans}}=\text{'N'}$.
${\mathbf{ijob}}=0$
Solve (1) only.
${\mathbf{ijob}}=1$
The functionality of ${\mathbf{ijob}}=0$ and $3$.
${\mathbf{ijob}}=2$
The functionality of ${\mathbf{ijob}}=0$ and $4$.
${\mathbf{ijob}}=3$
Only an estimate of $\mathrm{Dif}\left[\left(A,D\right),\left(B,E\right)\right]$ is computed based on the Frobenius norm.
${\mathbf{ijob}}=4$
Only an estimate of $\mathrm{Dif}\left[\left(A,D\right),\left(B,E\right)\right]$ is computed based on the $1$-norm.
If ${\mathbf{trans}}=\text{'T'}$, ijob is not referenced.
Constraint: if ${\mathbf{trans}}=\text{'N'}$, $0\le {\mathbf{ijob}}\le 4$.
3:     $\mathrm{a}\left(\mathit{lda},:\right)$ – double array
The first dimension of the array a must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$.
The second dimension of the array a must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$.
The upper quasi-triangular matrix $A$.
4:     $\mathrm{b}\left(\mathit{ldb},:\right)$ – double array
The first dimension of the array b must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The second dimension of the array b must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The upper quasi-triangular matrix $B$.
5:     $\mathrm{c}\left(\mathit{ldc},:\right)$ – double array
The first dimension of the array c must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$.
The second dimension of the array c must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
Contains the right-hand-side matrix $C$.
6:     $\mathrm{d}\left(\mathit{ldd},:\right)$ – double array
The first dimension of the array d must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$.
The second dimension of the array d must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$.
The upper triangular matrix $D$.
7:     $\mathrm{e}\left(\mathit{lde},:\right)$ – double array
The first dimension of the array e must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The second dimension of the array e must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The upper triangular matrix $E$.
8:     $\mathrm{f}\left(\mathit{ldf},:\right)$ – double array
The first dimension of the array f must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$.
The second dimension of the array f must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
Contains the right-hand side matrix $F$.

### Optional Input Parameters

1:     $\mathrm{m}$int64int32nag_int scalar
Default: the first dimension of the arrays a, c, d, f and the second dimension of the arrays a, d.
$m$, the order of the matrices $A$ and $D$, and the row dimension of the matrices $C$, $F$, $R$ and $L$.
Constraint: ${\mathbf{m}}\ge 0$.
2:     $\mathrm{n}$int64int32nag_int scalar
Default: the first dimension of the arrays b, e and the second dimension of the arrays b, c, e, f.
$n$, the order of the matrices $B$ and $E$, and the column dimension of the matrices $C$, $F$, $R$ and $L$.
Constraint: ${\mathbf{n}}\ge 0$.

### Output Parameters

1:     $\mathrm{c}\left(\mathit{ldc},:\right)$ – double array
The first dimension of the array c will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$.
The second dimension of the array c will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
If ${\mathbf{ijob}}=0$, $1$ or $2$, c stores the solution matrix $R$.
If ${\mathbf{trans}}=\text{'N'}$ and ${\mathbf{ijob}}=3$ or $4$, c holds $R$, the solution achieved during the computation of the Dif estimate.
2:     $\mathrm{f}\left(\mathit{ldf},:\right)$ – double array
The first dimension of the array f will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$.
The second dimension of the array f will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
If ${\mathbf{ijob}}=0$, $1$ or $2$, f stores the solution matrix $L$.
If ${\mathbf{trans}}=\text{'N'}$ and ${\mathbf{ijob}}=3$ or $4$, f holds $L$, the solution achieved during the computation of the Dif estimate.
3:     $\mathrm{scale}$ – double scalar
$\alpha$, the scaling factor in (1) or (2).
If $0<{\mathbf{scale}}<1$, c and f hold the solutions $R$ and $L$, respectively, to a slightly perturbed system but the input arrays a, b, d and e have not been changed.
If ${\mathbf{scale}}=0$, c and f hold the solutions $R$ and $L$, respectively, to the homogeneous system with $C=F=0$. In this case dif is not referenced.
Normally, ${\mathbf{scale}}=1$.
4:     $\mathrm{dif}$ – double scalar
The estimate of $\mathrm{Dif}$. If ${\mathbf{ijob}}=0$, dif is not referenced.
5:     $\mathrm{info}$int64int32nag_int scalar
${\mathbf{info}}=0$ unless the function detects an error (see Error Indicators and Warnings).

## Error Indicators and Warnings

${\mathbf{info}}=-i$
If ${\mathbf{info}}=-i$, parameter $i$ had an illegal value on entry. The parameters are numbered as follows:
1: trans, 2: ijob, 3: m, 4: n, 5: a, 6: lda, 7: b, 8: ldb, 9: c, 10: ldc, 11: d, 12: ldd, 13: e, 14: lde, 15: f, 16: ldf, 17: scale, 18: dif, 19: work, 20: lwork, 21: iwork, 22: info.
It is possible that info refers to a parameter that is omitted from the MATLAB interface. This usually indicates that an error in one of the other input parameters has caused an incorrect value to be inferred.
${\mathbf{info}}>0$
$\left(A,D\right)$ and $\left(B,E\right)$ have common or close eigenvalues and so no solution could be computed.

## Accuracy

See Kågström (1994) for a perturbation analysis of the generalized Sylvester equation.

The total number of floating-point operations needed to solve the generalized Sylvester equations is approximately $2mn\left(n+m\right)$. The Frobenius norm estimate of $\mathrm{Dif}$ does not require additional significant computation, but the $1$-norm estimate is typically five times more expensive.
The complex analogue of this function is nag_lapack_ztgsyl (f08yv).

## Example

This example solves the generalized Sylvester equations
 $AR-LB =αC DR-LE =αF,$
where
 $A = 4.0 1.0 1.0 2.0 0.0 3.0 4.0 1.0 0.0 1.0 3.0 1.0 0.0 0.0 0.0 6.0 , B= 1.0 1.0 1.0 1.0 0.0 3.0 4.0 1.0 0.0 1.0 3.0 1.0 0.0 0.0 0.0 4.0 ,$
 $D = 2.0 1.0 1.0 3.0 0.0 1.0 2.0 1.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 2.0 , E= 1.0 1.0 1.0 2.0 0.0 1.0 4.0 1.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 1.0 ,$
 $C = -4.0 7.0 1.0 12.0 -9.0 2.0 -2.0 -2.0 -4.0 2.0 -2.0 8.0 -7.0 7.0 -6.0 19.0 and F= -7.0 5.0 0.0 7.0 -5.0 1.0 -8.0 0.0 -1.0 2.0 -3.0 5.0 -3.0 2.0 0.0 5.0 .$
```function f08yh_example

fprintf('f08yh example results\n\n');

% Given generalized Schur form pairs (A,D) and (B,E), and matrices C and F:
% solve AR - LB = alpha*C
%       DR - LE = alpha*F
% for R, L and alpha

% Schur form pairs
A = [4, 1, 1, 2;
0, 3, 4, 1;
0, 1, 3, 1;
0, 0, 0, 6];
D = [2, 1, 1, 3;
0, 1, 2, 1;
0, 0, 1, 1;
0, 0, 0, 2];
B = [1, 1, 1, 1;
0, 3, 4, 1;
0, 1, 3, 1;
0, 0, 0, 4];
E = [1, 1, 1, 2;
0, 1, 4, 1;
0, 0, 1, 1;
0, 0, 0, 1];
% Right hand matrices
C = [-4, 7,  1, 12;
-9, 2, -2, -2;
-4, 2, -2,  8;
-7, 7, -6, 19];
F = [-7, 5,  0,  7;
-5, 1, -8,  0;
-1, 2, -3,  5;
-3, 2,  0,  5];

% Solve
trans = 'No transpose';
ijob = int64(0);
[R, L, alpha, ~, info] = f08yh( ...
trans, ijob, A, B, C, D, E, F);

[ifail] = x04ca( ...
'Gen', ' ', R, 'Solution matrix R');

fprintf('\n');
[ifail] = x04ca( ...
'Gen', ' ', L, 'Solution matrix L');

fprintf('\nalpha = %10.2e\n', alpha);

```
```f08yh example results

Solution matrix R
1          2          3          4
1      1.0000     1.0000     1.0000     1.0000
2     -1.0000     2.0000    -1.0000    -1.0000
3     -1.0000     1.0000     3.0000     1.0000
4     -1.0000     1.0000    -1.0000     4.0000

Solution matrix L
1          2          3          4
1      4.0000    -1.0000     1.0000    -1.0000
2      1.0000     3.0000    -1.0000     1.0000
3     -1.0000     1.0000     2.0000    -1.0000
4      1.0000    -1.0000     1.0000     1.0000

alpha =   1.00e+00
```