F08QHF (DTRSYL) (PDF version)
F08 Chapter Contents
F08 Chapter Introduction
NAG Library Manual

NAG Library Routine Document

F08QHF (DTRSYL)

Note:  before using this routine, please read the Users' Note for your implementation to check the interpretation of bold italicised terms and other implementation-dependent details.

+ Contents

    1  Purpose
    7  Accuracy

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
opAX ± XopB = αC ,
where opA = A  or AT, and the matrices A and B are upper quasi-triangular matrices in canonical Schur form (as returned by F08PEF (DHSEQR)); α is a scale factor (1) determined by the routine 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 αi±βj0, where αi and βj are the eigenvalues of A and B respectively and the sign (+ or -) 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 AX-XB=C Numerical Analysis Report University of Manchester

5  Parameters

1:     TRANA – CHARACTER(1)Input
On entry: specifies the option opA.
TRANA='N'
opA=A.
TRANA='T' or 'C'
opA=AT.
Constraint: TRANA='N', 'T' or 'C'.
2:     TRANB – CHARACTER(1)Input
On entry: specifies the option opB.
TRANB='N'
opB=B.
TRANB='T' or 'C'
opB=BT.
Constraint: TRANB='N', 'T' or 'C'.
3:     ISGN – INTEGERInput
On entry: indicates the form of the Sylvester equation.
ISGN=+1
The equation is of the form opAX+XopB=αC.
ISGN=-1
The equation is of the form opAX-XopB=αC.
Constraint: ISGN=+1 or -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: M0.
5:     N – INTEGERInput
On entry: n, the order of the matrix B, and the number of columns in the matrices X and C.
Constraint: N0.
6:     A(LDA,*) – REAL (KIND=nag_wp) arrayInput
Note: the second dimension of the array A must be at least max1,M.
On entry: the m by 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: LDAmax1,M.
8:     B(LDB,*) – REAL (KIND=nag_wp) arrayInput
Note: the second dimension of the array B must be at least max1,N.
On entry: the n by 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: LDBmax1,N.
10:   C(LDC,*) – REAL (KIND=nag_wp) arrayInput/Output
Note: the second dimension of the array C must be at least max1,N.
On entry: the m by 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: LDCmax1,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 AX-XB=C. (To apply the remarks to the equation AX+XB=C, simply replace B by -B.)
Let 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, 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 sepA,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 mnm+n.
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 F08PAF (DGEES), for example):
A = Q1 A~ Q1T   and   B = Q2 B~ Q2T
where A~ and B~ are upper quasi-triangular and Q1 and Q2 are orthogonal. The original equation may then be transformed to:
A~ X~ ± X~ B~ = C~
where X~ = Q1T X Q2  and C~ = Q1T C Q2 . 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 = Q1 X~ Q2T .
The complex analogue of this routine is F08QVF (ZTRSYL).

9  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 .

9.1  Program Text

Program Text (f08qhfe.f90)

9.2  Program Data

Program Data (f08qhfe.d)

9.3  Program Results

Program Results (f08qhfe.r)


F08QHF (DTRSYL) (PDF version)
F08 Chapter Contents
F08 Chapter Introduction
NAG Library Manual

© The Numerical Algorithms Group Ltd, Oxford, UK. 2012