1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
|
<HTML>
<HEAD><TITLE>SB03PD - SLICOT Library Routine Documentation</TITLE>
</HEAD>
<BODY>
<H2><A Name="SB03PD">SB03PD</A></H2>
<H3>
Solution of discrete-time Lyapunov equations and separation estimation
</H3>
<A HREF ="#Specification"><B>[Specification]</B></A>
<A HREF ="#Arguments"><B>[Arguments]</B></A>
<A HREF ="#Method"><B>[Method]</B></A>
<A HREF ="#References"><B>[References]</B></A>
<A HREF ="#Comments"><B>[Comments]</B></A>
<A HREF ="#Example"><B>[Example]</B></A>
<P>
<B><FONT SIZE="+1">Purpose</FONT></B>
<PRE>
To solve the real discrete Lyapunov matrix equation
op(A)'*X*op(A) - X = scale*C
and/or estimate the quantity, called separation,
sepd(op(A),op(A)') = min norm(op(A)'*X*op(A) - X)/norm(X)
where op(A) = A or A' (A**T) and C is symmetric (C = C').
(A' denotes the transpose of the matrix A.) A is N-by-N, the right
hand side C and the solution X are N-by-N, and scale is an output
scale factor, set less than or equal to 1 to avoid overflow in X.
</PRE>
<A name="Specification"><B><FONT SIZE="+1">Specification</FONT></B></A>
<PRE>
SUBROUTINE SB03PD( JOB, FACT, TRANA, N, A, LDA, U, LDU, C, LDC,
$ SCALE, SEPD, FERR, WR, WI, IWORK, DWORK,
$ LDWORK, INFO )
C .. Scalar Arguments ..
CHARACTER FACT, JOB, TRANA
INTEGER INFO, LDA, LDC, LDU, LDWORK, N
DOUBLE PRECISION FERR, SCALE, SEPD
C .. Array Arguments ..
INTEGER IWORK( * )
DOUBLE PRECISION A( LDA, * ), C( LDC, * ), DWORK( * ),
$ U( LDU, * ), WI( * ), WR( * )
</PRE>
<A name="Arguments"><B><FONT SIZE="+1">Arguments</FONT></B></A>
<P>
<B>Mode Parameters</B>
<PRE>
JOB CHARACTER*1
Specifies the computation to be performed, as follows:
= 'X': Compute the solution only;
= 'S': Compute the separation only;
= 'B': Compute both the solution and the separation.
FACT CHARACTER*1
Specifies whether or not the real Schur factorization
of the matrix A is supplied on entry, as follows:
= 'F': On entry, A and U contain the factors from the
real Schur factorization of the matrix A;
= 'N': The Schur factorization of A will be computed
and the factors will be stored in A and U.
TRANA CHARACTER*1
Specifies the form of op(A) to be used, as follows:
= 'N': op(A) = A (No transpose);
= 'T': op(A) = A**T (Transpose);
= 'C': op(A) = A**T (Conjugate transpose = Transpose).
</PRE>
<B>Input/Output Parameters</B>
<PRE>
N (input) INTEGER
The order of the matrices A, X, and C. N >= 0.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading N-by-N part of this array must
contain the matrix A. If FACT = 'F', then A contains
an upper quasi-triangular matrix in Schur canonical form.
On exit, if INFO = 0 or INFO = N+1, the leading N-by-N
part of this array contains the upper quasi-triangular
matrix in Schur canonical form from the Shur factorization
of A. The contents of array A is not modified if
FACT = 'F'.
LDA INTEGER
The leading dimension of array A. LDA >= MAX(1,N).
U (input or output) DOUBLE PRECISION array, dimension
(LDU,N)
If FACT = 'F', then U is an input argument and on entry
it must contain the orthogonal matrix U from the real
Schur factorization of A.
If FACT = 'N', then U is an output argument and on exit,
if INFO = 0 or INFO = N+1, it contains the orthogonal
N-by-N matrix from the real Schur factorization of A.
LDU INTEGER
The leading dimension of array U. LDU >= MAX(1,N).
C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
On entry with JOB = 'X' or 'B', the leading N-by-N part of
this array must contain the symmetric matrix C.
On exit with JOB = 'X' or 'B', if INFO = 0 or INFO = N+1,
the leading N-by-N part of C has been overwritten by the
symmetric solution matrix X.
If JOB = 'S', C is not referenced.
LDC INTEGER
The leading dimension of array C.
LDC >= 1, if JOB = 'S';
LDC >= MAX(1,N), otherwise.
SCALE (output) DOUBLE PRECISION
The scale factor, scale, set less than or equal to 1 to
prevent the solution overflowing.
SEPD (output) DOUBLE PRECISION
If JOB = 'S' or JOB = 'B', and INFO = 0 or INFO = N+1,
SEPD contains the estimate in the 1-norm of
sepd(op(A),op(A)').
If JOB = 'X' or N = 0, SEPD is not referenced.
FERR (output) DOUBLE PRECISION
If JOB = 'B', and INFO = 0 or INFO = N+1, FERR contains
an estimated forward error bound for the solution X.
If XTRUE is the true solution, FERR bounds the relative
error in the computed solution, measured in the Frobenius
norm: norm(X - XTRUE)/norm(XTRUE).
If JOB = 'X' or JOB = 'S', FERR is not referenced.
WR (output) DOUBLE PRECISION array, dimension (N)
WI (output) DOUBLE PRECISION array, dimension (N)
If FACT = 'N', and INFO = 0 or INFO = N+1, WR and WI
contain the real and imaginary parts, respectively, of the
eigenvalues of A.
If FACT = 'F', WR and WI are not referenced.
</PRE>
<B>Workspace</B>
<PRE>
IWORK INTEGER array, dimension (N*N)
This array is not referenced if JOB = 'X'.
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0 or INFO = N+1, DWORK(1) returns the
optimal value of LDWORK.
LDWORK INTEGER
The length of the array DWORK. LDWORK >= 1 and
If JOB = 'X' then
If FACT = 'F', LDWORK >= MAX(N*N,2*N);
If FACT = 'N', LDWORK >= MAX(N*N,3*N).
If JOB = 'S' or JOB = 'B' then
LDWORK >= 2*N*N + 2*N.
For optimum performance LDWORK should be larger.
</PRE>
<B>Error Indicator</B>
<PRE>
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value;
> 0: if INFO = i, the QR algorithm failed to compute all
the eigenvalues (see LAPACK Library routine DGEES);
elements i+1:n of WR and WI contain eigenvalues
which have converged, and A contains the partially
converged Schur form;
= N+1: if matrix A has almost reciprocal eigenvalues;
perturbed values were used to solve the equation
(but the matrix A is unchanged).
</PRE>
<A name="Method"><B><FONT SIZE="+1">Method</FONT></B></A>
<PRE>
After reducing matrix A to real Schur canonical form (if needed),
a discrete-time version of the Bartels-Stewart algorithm is used.
A set of equivalent linear algebraic systems of equations of order
at most four are formed and solved using Gaussian elimination with
complete pivoting.
</PRE>
<A name="References"><B><FONT SIZE="+1">References</FONT></B></A>
<PRE>
[1] Barraud, A.Y. T
A numerical algorithm to solve A XA - X = Q.
IEEE Trans. Auto. Contr., AC-22, pp. 883-885, 1977.
[2] Bartels, R.H. and Stewart, G.W. T
Solution of the matrix equation A X + XB = C.
Comm. A.C.M., 15, pp. 820-826, 1972.
</PRE>
<A name="Numerical Aspects"><B><FONT SIZE="+1">Numerical Aspects</FONT></B></A>
<PRE> 3
The algorithm requires 0(N ) operations.
</PRE>
<A name="Comments"><B><FONT SIZE="+1">Further Comments</FONT></B></A>
<PRE>
SEPD is defined as
sepd( op(A), op(A)' ) = sigma_min( T )
where sigma_min(T) is the smallest singular value of the
N*N-by-N*N matrix
T = kprod( op(A)', op(A)' ) - I(N**2).
I(N**2) is an N*N-by-N*N identity matrix, and kprod denotes the
Kronecker product. The program estimates sigma_min(T) by the
reciprocal of an estimate of the 1-norm of inverse(T). The true
reciprocal 1-norm of inverse(T) cannot differ from sigma_min(T) by
more than a factor of N.
When SEPD is small, small changes in A, C can cause large changes
in the solution of the equation. An approximate bound on the
maximum relative error in the computed solution is
EPS * norm(A)**2 / SEPD
where EPS is the machine precision.
</PRE>
<A name="Example"><B><FONT SIZE="+1">Example</FONT></B></A>
<P>
<B>Program Text</B>
<PRE>
None
</PRE>
<B>Program Data</B>
<PRE>
None
</PRE>
<B>Program Results</B>
<PRE>
None
</PRE>
<HR>
<p>
<A HREF=..\libindex.html><B>Return to index</B></A></BODY>
</HTML>
|
