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_dtrsyl (f08qh)

## Purpose

nag_lapack_dtrsyl (f08qh) solves the real quasi-triangular Sylvester matrix equation.

## Syntax

[c, scale, info] = f08qh(trana, tranb, isgn, a, b, c, 'm', m, 'n', n)
[c, scale, info] = nag_lapack_dtrsyl(trana, tranb, isgn, a, b, c, 'm', m, 'n', n)

## Description

nag_lapack_dtrsyl (f08qh) solves the real Sylvester matrix equation
 $opAX ± XopB = αC ,$
where $\mathrm{op}\left(A\right)=A$ or ${A}^{\mathrm{T}}$, and the matrices $A$ and $B$ are upper quasi-triangular matrices in canonical Schur form (as returned by nag_lapack_dhseqr (f08pe)); $\alpha$ is a scale factor ($\text{}\le 1$) determined by the function to avoid overflow in $X$; $A$ is $m$ by $m$ and $B$ is $n$ by $n$ while the right-hand side matrix $C$ and the solution matrix $X$ are both $m$ by $n$. The matrix $X$ is obtained by a straightforward process of back-substitution (see Golub and Van Loan (1996)).
Note that the equation has a unique solution if and only if ${\alpha }_{i}±{\beta }_{j}\ne 0$, where $\left\{{\alpha }_{i}\right\}$ and $\left\{{\beta }_{j}\right\}$ are the eigenvalues of $A$ and $B$ respectively and the sign ($+$ or $-$) is the same as that used in the equation to be solved.

## References

Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore
Higham N J (1992) Perturbation theory and backward error for $AX-XB=C$ Numerical Analysis Report University of Manchester

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{trana}$ – string (length ≥ 1)
Specifies the option $\mathrm{op}\left(A\right)$.
${\mathbf{trana}}=\text{'N'}$
$\mathrm{op}\left(A\right)=A$.
${\mathbf{trana}}=\text{'T'}$ or $\text{'C'}$
$\mathrm{op}\left(A\right)={A}^{\mathrm{T}}$.
Constraint: ${\mathbf{trana}}=\text{'N'}$, $\text{'T'}$ or $\text{'C'}$.
2:     $\mathrm{tranb}$ – string (length ≥ 1)
Specifies the option $\mathrm{op}\left(B\right)$.
${\mathbf{tranb}}=\text{'N'}$
$\mathrm{op}\left(B\right)=B$.
${\mathbf{tranb}}=\text{'T'}$ or $\text{'C'}$
$\mathrm{op}\left(B\right)={B}^{\mathrm{T}}$.
Constraint: ${\mathbf{tranb}}=\text{'N'}$, $\text{'T'}$ or $\text{'C'}$.
3:     $\mathrm{isgn}$int64int32nag_int scalar
Indicates the form of the Sylvester equation.
${\mathbf{isgn}}=+1$
The equation is of the form $\mathrm{op}\left(A\right)X+X\mathrm{op}\left(B\right)=\alpha C$.
${\mathbf{isgn}}=-1$
The equation is of the form $\mathrm{op}\left(A\right)X-X\mathrm{op}\left(B\right)=\alpha C$.
Constraint: ${\mathbf{isgn}}=+1$ or $-1$.
4:     $\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 $m$ by $m$ upper quasi-triangular matrix $A$ in canonical Schur form, as returned by nag_lapack_dhseqr (f08pe).
5:     $\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 $n$ by $n$ upper quasi-triangular matrix $B$ in canonical Schur form, as returned by nag_lapack_dhseqr (f08pe).
6:     $\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)$.
The $m$ by $n$ right-hand side matrix $C$.

### Optional Input Parameters

1:     $\mathrm{m}$int64int32nag_int scalar
Default: the first dimension of the array a and the second dimension of the array a. (An error is raised if these dimensions are not equal.)
$m$, the order of the matrix $A$, and the number of rows in the matrices $X$ and $C$.
Constraint: ${\mathbf{m}}\ge 0$.
2:     $\mathrm{n}$int64int32nag_int scalar
Default: the first dimension of the array b and the second dimension of the array b. (An error is raised if these dimensions are not equal.)
$n$, the order of the matrix $B$, and the number of columns in the matrices $X$ and $C$.
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)$.
c stores the solution matrix $X$.
2:     $\mathrm{scale}$ – double scalar
The value of the scale factor $\alpha$.
3:     $\mathrm{info}$int64int32nag_int scalar
${\mathbf{info}}=0$ unless the function detects an error (see Error Indicators and Warnings).

## Error Indicators and Warnings

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

${\mathbf{info}}=-i$
If ${\mathbf{info}}=-i$, parameter $i$ had an illegal value on entry. The parameters are numbered as follows:
1: trana, 2: tranb, 3: isgn, 4: m, 5: n, 6: a, 7: lda, 8: b, 9: ldb, 10: c, 11: ldc, 12: scale, 13: 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.
W  ${\mathbf{info}}=1$
$A$ and $B$ have common or close eigenvalues, perturbed values of which were used to solve the equation.

## Accuracy

Consider the equation $AX-XB=C$. (To apply the remarks to the equation $AX+XB=C$, simply replace $B$ by $-B$.)
Let $\stackrel{~}{X}$ be the computed solution and $R$ the residual matrix:
 $R = C - AX~ - X~B .$
Then the residual is always small:
 $RF = Oε AF + BF X~F .$
However, $\stackrel{~}{X}$ is not necessarily the exact solution of a slightly perturbed equation; in other words, the solution is not backwards stable.
For the forward error, the following bound holds:
 $X~ - X F ≤ RF sep A,B$
but this may be a considerable over estimate. See Golub and Van Loan (1996) for a definition of $\mathit{sep}\left(A,B\right)$, and Higham (1992) for further details.
These remarks also apply to the solution of a general Sylvester equation, as described in Further Comments.

The total number of floating-point operations is approximately $mn\left(m+n\right)$.
To solve the general real Sylvester equation
 $AX ± XB = C$
where $A$ and $B$ are general nonsymmetric matrices, $A$ and $B$ must first be reduced to Schur form (by calling nag_lapack_dgees (f08pa), for example):
 $A = Q1 A~ Q1T and B = Q2 B~ Q2T$
where $\stackrel{~}{A}$ and $\stackrel{~}{B}$ are upper quasi-triangular and ${Q}_{1}$ and ${Q}_{2}$ are orthogonal. The original equation may then be transformed to:
 $A~ X~ ± X~ B~ = C~$
where $\stackrel{~}{X}={Q}_{1}^{\mathrm{T}}X{Q}_{2}$ and $\stackrel{~}{C}={Q}_{1}^{\mathrm{T}}C{Q}_{2}$. $\stackrel{~}{C}$ may be computed by matrix multiplication; nag_lapack_dtrsyl (f08qh) may be used to solve the transformed equation; and the solution to the original equation can be obtained as $X={Q}_{1}\stackrel{~}{X}{Q}_{2}^{\mathrm{T}}$.
The complex analogue of this function is nag_lapack_ztrsyl (f08qv).

## Example

This example solves the Sylvester equation $AX+XB=C$, where
 $A = 0.10 0.50 0.68 -0.21 -0.50 0.10 -0.24 0.67 0.00 0.00 0.19 -0.35 0.00 0.00 0.00 -0.72 ,$
 $B = -0.99 -0.17 0.39 0.58 0.00 0.48 -0.84 -0.15 0.00 0.00 0.75 0.25 0.00 0.00 -0.25 0.75$
and
 $C = 0.63 -0.56 0.08 -0.23 -0.45 -0.31 0.27 1.21 0.20 -0.35 0.41 0.84 0.49 -0.05 -0.52 -0.08 .$
```function f08qh_example

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

% Quasi-triangular (Schur form) matrices A and B, and general matrix C
a = [ 0.10,  0.50,  0.68, -0.21;
-0.50,  0.10, -0.24,  0.67;
0,     0,     0.19, -0.35;
0,     0,     0,    -0.72];
b = [-0.99, -0.17,  0.39,  0.58;
0,     0.48, -0.84, -0.15;
0,     0,     0.75,  0.25;
0,     0,    -0.25,  0.75];
c = [ 0.63, -0.56,  0.08, -0.23;
-0.45, -0.31,  0.27,  1.21;
0.20, -0.35,  0.41,  0.84;
0.49, -0.05, -0.52, -0.08];

% Solve Sylvester equation AX + XB = C for X
trana = 'No transpose';
tranb = 'No transpose';
isgn = int64(1);
[X, scale, info] = f08qh( ...
trana, tranb, isgn, a, b, c);

disp('Solution matrix');
disp(X);

```
```f08qh example results

Solution matrix
-0.4209    0.1764    0.2438   -0.9577
0.5600   -0.8337   -0.7221    0.5386
-0.1246   -0.3392    0.6221    0.8691
-0.2865    0.4113    0.5535    0.3174

```