NAG Library Routine Document

F08QHF (DTRSYL)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

+− 9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

1 Purpose

F08QHF (DTRSYL) solves the real quasi-triangular Sylvester matrix equation.

2 Specification

SUBROUTINE F08QHF (

TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, LDC, SCALE, INFO)

INTEGER	ISGN, M, N, LDA, LDB, LDC, INFO
REAL (KIND=nag_wp)	A(LDA,), B(LDB,), C(LDC,*), SCALE
CHARACTER(1)	TRANA, TRANB

The routine may be called by its LAPACK name dtrsyl.

3 Description

F08QHF (DTRSYL) solves the real Sylvester matrix equation

op (A) X \pm X op (B) = α C,

where

op (A) = A

A^{T}

, and the matrices

A

and

B

are upper quasi-triangular matrices in canonical Schur form (as returned by F08PEF (DHSEQR));

α

is a scale factor (

\leq 1

) determined by the routine to avoid overflow in

X

;

A

m

m

and

B

n

n

while the right-hand side matrix

C

and the solution matrix

X

are both

m

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

α_{i} \pm β_{j} \neq 0

, where

\{α_{i}\}

and

\{β_{j}\}

are the eigenvalues of

A

and

B

respectively and the sign (

+

-

) is the same as that used in the equation to be solved.

4 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

A X - X B = C

Numerical Analysis Report University of Manchester

5 Parameters

1: TRANA – CHARACTER(1)Input

On entry: specifies the option

op (A)

$TRANA ='N'$: $op (A) = A$ .
$TRANA ='T'$ or $'C'$: $op (A) = A^{T}$ .

Constraint:

TRANA ='N'

'T'

'C'

2: TRANB – CHARACTER(1)Input

On entry: specifies the option

op (B)

$TRANB ='N'$: $op (B) = B$ .
$TRANB ='T'$ or $'C'$: $op (B) = B^{T}$ .

Constraint:

TRANB ='N'

'T'

'C'

3: ISGN – INTEGERInput

On entry: indicates the form of the Sylvester equation.

$ISGN = +1$: The equation is of the form $op (A) X + X op (B) = α C$ .
$ISGN = -1$: The equation is of the form $op (A) X - X op (B) = α C$ .

Constraint:

ISGN = +1

-1

4: M – INTEGERInput

On entry:

m

, the order of the matrix

A

, and the number of rows in the matrices

X

and

C

Constraint:

M \geq 0

5: N – INTEGERInput

On entry:

n

, the order of the matrix

B

, and the number of columns in the matrices

X

and

C

Constraint:

N \geq 0

6: A(LDA, $*$ ) – REAL (KIND=nag_wp) arrayInput

Note: the second dimension of the array A must be at least

\max (1, M)

On entry: the

m

m

upper quasi-triangular matrix

A

in canonical Schur form, as returned by F08PEF (DHSEQR).

7: LDA – INTEGERInput

On entry: the first dimension of the array A as declared in the (sub)program from which F08QHF (DTRSYL) is called.

Constraint:

LDA \geq \max (1, M)

8: B(LDB, $*$ ) – REAL (KIND=nag_wp) arrayInput

Note: the second dimension of the array B must be at least

\max (1, N)

On entry: the

n

n

upper quasi-triangular matrix

B

in canonical Schur form, as returned by F08PEF (DHSEQR).

9: LDB – INTEGERInput

On entry: the first dimension of the array B as declared in the (sub)program from which F08QHF (DTRSYL) is called.

Constraint:

LDB \geq \max (1, N)

10: C(LDC, $*$ ) – REAL (KIND=nag_wp) arrayInput/Output

Note: the second dimension of the array C must be at least

\max (1, N)

On entry: the

m

n

right-hand side matrix

C

On exit: C is overwritten by the solution matrix

X

11: LDC – INTEGERInput

On entry: the first dimension of the array C as declared in the (sub)program from which F08QHF (DTRSYL) is called.

Constraint:

LDC \geq \max (1, M)

12: SCALE – REAL (KIND=nag_wp)Output

On exit: the value of the scale factor

α

13: INFO – INTEGEROutput

On exit:

INFO = 0

unless the routine detects an error (see Section 6).

6 Error Indicators and Warnings

$INFO < 0$: If $INFO = - i$ , argument $i$ had an illegal value. An explanatory message is output, and execution of the program is terminated.

$INFO = 1$: $A$ and $B$ have common or close eigenvalues, perturbed values of which were used to solve the equation.

7 Accuracy

Consider the equation

A X - X B = C

. (To apply the remarks to the equation

A X + X B = C

, simply replace

B

- B

Let

\tilde{X}

be the computed solution and

R

the residual matrix:

R = C - (A \tilde{X} - \tilde{X} B) .

Then the residual is always small:

{‖R‖}_{F} = O (ε) ({‖A‖}_{F} + {‖B‖}_{F}) {‖\tilde{X}‖}_{F} .

However,

\tilde{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:

{‖\tilde{X} - X‖}_{F} \leq \frac{{‖R‖}_{F}}{sep (A, B)}

but this may be a considerable over estimate. See Golub and Van Loan (1996) for a definition of

sep (A, B)

, and Higham (1992) for further details.

These remarks also apply to the solution of a general Sylvester equation, as described in Section 8.

8 Further Comments

The total number of floating point operations is approximately

m n (m + n)

To solve the general real Sylvester equation

A X \pm X B = C

where

A

and

B

are general nonsymmetric matrices,

A

and

B

must first be reduced to Schur form (by calling F08PAF (DGEES), for example):

A = Q_{1} \tilde{A} Q_{1}^{T} and B = Q_{2} \tilde{B} Q_{2}^{T}

where

\tilde{A}

and

\tilde{B}

are upper quasi-triangular and

Q_{1}

and

Q_{2}

are orthogonal. The original equation may then be transformed to:

\tilde{A} \tilde{X} \pm \tilde{X} \tilde{B} = \tilde{C}

where

\tilde{X} = Q_{1}^{T} X Q_{2}

and

\tilde{C} = Q_{1}^{T} C Q_{2}

\tilde{C}

may be computed by matrix multiplication; F08QHF (DTRSYL) may be used to solve the transformed equation; and the solution to the original equation can be obtained as

X = Q_{1} \tilde{X} Q_{2}^{T}

The complex analogue of this routine is F08QVF (ZTRSYL).

9 Example

This example solves the Sylvester equation

A X + X B = C

, where

A = (\begin{array}{r} 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 \end{array}),

B = (\begin{array}{r} - 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 \end{array})

and

C = (\begin{array}{r} 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 \end{array}) .

NAG Library Routine DocumentF08QHF (DTRSYL)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

NAG Library Routine Document

F08QHF (DTRSYL)