-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathfstdemo2.m
223 lines (151 loc) · 6.04 KB
/
fstdemo2.m
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
echo on
% This demonstration shows the use of the SLICOT-based Matlab function
% for computing the generator and the Cholesky factor for the inverse
% of a positive definite block Toeplitz matrix.
echo off
% RELEASE 2.0 of SLICOT Basic Systems and Control Toolbox.
% Based on SLICOT RELEASE 5.7, Copyright (c) 2002-2020 NICONET e.V.
%
% Contributor:
% V. Sima, Research Institute for Informatics, Bucharest, Dec. 2000.
%
% Revisions: V. Sima, March 2009.
%
echo on
global pause_wait % This could be used in pause(n) command.
% If pause_wait < 0, standard command pause is used (default).
% Any key should then be pressed to continue.
if ~exist('pause_wait', 'var') || isempty(pause_wait), pause_wait = -1; end
if pause_wait < 0, pause, else pause(pause_wait), end
% Consider an n-by-n block Toeplitz matrix, BT, with k-by-k blocks.
% Let k = 2, n = 3. The first block-row of BT has n = 3 blocks:
% the block i has elements containing the figure i, i = 1, 2, 3.
k = 2; n = 3; Tr = [ [ 100 1; 1 100] [26 27; 28 29] [36 37; 38 39]];
Tc = Tr'; BT = btoeplitz(Tc,Tr);
echo off
disp(' '), disp('The first block-row, Tr, is')
Tr
echo on
if pause_wait < 0, pause, else pause(pause_wait), end
echo off
disp(' '), disp('The first block-column, Tc, is')
Tc
echo on
if pause_wait < 0, pause, else pause(pause_wait), end
echo off
disp(' '), disp('The block Toeplitz matrix, BT, is')
BT
echo on
if pause_wait < 0, pause, else pause(pause_wait), end
% Available SLICOT-based Matlab functions: fstchol, fstgen, fstupd, fstsol.
% These functions call the mexfile fstoep. Use function fstgen here.
if pause_wait < 0, pause, else pause(pause_wait), end
more on
help fstgen
more off
if pause_wait < 0, pause, else pause(pause_wait), end
% Call fstgen, given the first block-row of BT.
G = fstgen(Tr);
if pause_wait < 0, pause, else pause(pause_wait), end
echo off
disp(' '), disp('The generator of the inverse of BT, G, is')
G
echo on
if pause_wait < 0, pause, else pause(pause_wait), end
% The generator can be used to compute an approximate inverse of BT:
approxinv = btoeplitz(G(1:k,:)', G(1:k,1:k) * eye(k,size(G,2))) ...
* btoeplitz(G(1:k,:)', G(1:k,1:k) * eye(k,size(G,2)))' ...
- btoeplitz(G(k+1:2*k,:)', G(k+1:2*k,1:k) * eye(k,size(G,2))) ...
* btoeplitz(G(k+1:2*k,:)', G(k+1:2*k,1:k) * eye(k,size(G,2)))';
if pause_wait < 0, pause, else pause(pause_wait), end
echo off
disp(' '),
disp(['The error norm: norm( approxinv * BT - eye(size(BT)) ) = ',...
num2str(norm( approxinv * BT - eye(size(BT)) ))])
echo on
if pause_wait < 0, pause, else pause(pause_wait), end
% Compute also the Cholesky factor, L, for the inverse of BT.
[G,L] = fstgen(Tr);
if pause_wait < 0, pause, else pause(pause_wait), end
echo off
disp(' '), disp('The lower Cholesky factor of inv(BT), L, is')
L
echo on
if pause_wait < 0, pause, else pause(pause_wait), end
echo off
disp(' '), disp('Check with Matlab function chol, Lm = inv(chol(BT))''')
Lm = inv(chol(BT))';
Lm
disp(' '),
disp(['The relative error norm: norm(L - Lm)/norm(Lm) = ',...
num2str(norm(L - Lm)/norm(Lm))])
echo on
if pause_wait < 0, pause, else pause(pause_wait), end
% Use fstgen to solve the systems X*BT = B.
% First, generate X, compute B, and then solve the equations.
% The columns of X contain the first n*k*m natural numbers.
m = n; X0 = 1 : k*n*m; X0 = reshape(X0,m,k*n); B = X0*BT;
if pause_wait < 0, pause, else pause(pause_wait), end
[G,L,X] = fstgen(Tr,B);
if pause_wait < 0, pause, else pause(pause_wait), end
echo off
disp(' '), disp('The solution of X*BT = B is')
X
echo on
if pause_wait < 0, pause, else pause(pause_wait), end
echo off
disp(' '),
disp(['The relative error norm: norm(X - X0)/norm(X0) = ',...
num2str(norm(X - X0)/norm(X0))])
echo on
if pause_wait < 0, pause, else pause(pause_wait), end
% Call fstgen, given the first block-column of BT.
G = fstgen(Tc);
if pause_wait < 0, pause, else pause(pause_wait), end
% The generator can be used to compute an approximate inverse of BT:
approxinv = btoeplitz(G(:,1:k), G(1:k,1:k) * eye(k,size(G,1))) ...
* btoeplitz(G(:,1:k), G(1:k,1:k) * eye(k,size(G,1)))' ...
- btoeplitz(G(:,k+1:2*k), G(1:k,k+1:2*k) * eye(k,size(G,1))) ...
* btoeplitz(G(:,k+1:2*k), G(1:k,k+1:2*k) * eye(k,size(G,1)))';
if pause_wait < 0, pause, else pause(pause_wait), end
echo off
disp(' '),
disp(['The error norm: norm( approxinv * BT - eye(size(BT)) ) = ',...
num2str(norm( approxinv * BT - eye(size(BT)) ))])
echo on
if pause_wait < 0, pause, else pause(pause_wait), end
% Compute also the upper Cholesky factor, L, for the inverse of BT.
[G,L] = fstgen(Tc);
if pause_wait < 0, pause, else pause(pause_wait), end
echo off
disp(' '), disp('The uper Cholesky factor of inv(BT), L, is')
L
echo on
if pause_wait < 0, pause, else pause(pause_wait), end
echo off
disp(' '), disp('Check with Matlab function chol, Lm = inv(chol(BT))')
Lm = inv(chol(BT));
disp(' '),
disp(['The relative error norm: norm(L - Lm)/norm(Lm) = ',...
num2str(norm(L - Lm)/norm(Lm))])
echo on
if pause_wait < 0, pause, else pause(pause_wait), end
% Use fstgen to solve the systems BT*X = B.
% First, generate X, compute B, and then solve the equations.
% The columns of X contain the first n*k*m natural numbers.
X0 = 1 : k*n*m; X0 = reshape(X0,k*n,m); B = BT*X0;
if pause_wait < 0, pause, else pause(pause_wait), end
[G,L,X] = fstgen(Tc,B);
if pause_wait < 0, pause, else pause(pause_wait), end
echo off
disp(' '), disp('The solution of BT*X = B is')
X
echo on
if pause_wait < 0, pause, else pause(pause_wait), end
echo off
disp(' '),
disp(['The relative error norm: norm(X - X0)/norm(X0) = ',...
num2str(norm(X - X0)/norm(X0))])
echo on
if pause_wait < 0, pause, else pause(pause_wait), end
echo off