|
7017
|
1 ## Copyright (C) 1996, 2000, 2005, 2007 |
|
|
2 ## Auburn University. All rights reserved. |
|
3431
|
3 ## |
|
|
4 ## This file is part of Octave. |
|
|
5 ## |
|
|
6 ## Octave is free software; you can redistribute it and/or modify it |
|
7016
|
7 ## under the terms of the GNU General Public License as published by |
|
|
8 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
9 ## your option) any later version. |
|
3431
|
10 ## |
|
7016
|
11 ## Octave is distributed in the hope that it will be useful, but |
|
|
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
14 ## General Public License for more details. |
|
3431
|
15 ## |
|
|
16 ## You should have received a copy of the GNU General Public License |
|
7016
|
17 ## along with Octave; see the file COPYING. If not, see |
|
|
18 ## <http://www.gnu.org/licenses/>. |
|
3431
|
19 |
|
|
20 ## -*- texinfo -*- |
|
3502
|
21 ## @deftypefn {Function File} {} zgreduce (@var{sys}, @var{meps}) |
|
3431
|
22 ## Implementation of procedure REDUCE in (Emami-Naeini and Van Dooren, |
|
|
23 ## Automatica, # 1982). |
|
|
24 ## @end deftypefn |
|
|
25 |
|
|
26 function retsys = zgreduce (Asys, meps) |
|
|
27 |
|
7125
|
28 if (nargin != 2) |
|
|
29 print_usage (); |
|
|
30 endif |
|
|
31 |
|
3431
|
32 ## SYS_INTERNAL accesses members of system data structure |
|
|
33 |
|
7132
|
34 is_digital (Asys); # make sure it's pure digital/continuous |
|
3431
|
35 |
|
|
36 exit_1 = 0; # exit_1 = 1 or 2 on exit of loop |
|
|
37 |
|
7132
|
38 if (Asys.n + Asys.nz == 0) |
|
3431
|
39 exit_1 = 2; # there are no finite zeros |
|
|
40 endif |
|
|
41 |
|
|
42 while (! exit_1) |
|
7132
|
43 [Q, R, Pi] = qr (Asys.d); # compress rows of D |
|
3431
|
44 Asys.d = Q'*Asys.d; |
|
|
45 Asys.c = Q'*Asys.c; |
|
|
46 |
|
|
47 ## check row norms of Asys.d |
|
7132
|
48 [sig, tau] = zgrownorm (Asys.d, meps); |
|
3431
|
49 |
|
|
50 ## disp("=======================================") |
|
|
51 ## disp(["zgreduce: meps=",num2str(meps), ", sig=",num2str(sig), ... |
|
|
52 ## ", tau=",num2str(tau)]) |
|
|
53 ## sysout(Asys) |
|
|
54 |
|
7132
|
55 if (tau == 0) |
|
3431
|
56 exit_1 = 1; # exit_1 - reduction complete and correct |
|
|
57 else |
|
|
58 Cb = Db = []; |
|
7132
|
59 if (sig) |
|
3431
|
60 Cb = Asys.c(1:sig,:); |
|
|
61 Db = Asys.d(1:sig,:); |
|
|
62 endif |
|
7132
|
63 Ct = Asys.c(sig+(1:tau),:); |
|
3431
|
64 |
|
|
65 ## compress columns of Ct |
|
7132
|
66 [pp, nn] = size (Ct); |
|
3431
|
67 rvec = nn:-1:1; |
|
7132
|
68 [V, Sj, Pi] = qr (Ct'); |
|
3431
|
69 V = V(:,rvec); |
|
7132
|
70 [rho, gnu] = zgrownorm (Sj, meps); |
|
3431
|
71 |
|
|
72 ## disp(["zgreduce: rho=",num2str(rho),", gnu=",num2str(gnu)]) |
|
|
73 ## Cb |
|
|
74 ## Db |
|
|
75 ## Ct |
|
|
76 ## Sj' |
|
|
77 |
|
7132
|
78 if (rho == 0) |
|
3431
|
79 exit_1 = 1; # exit_1 - reduction complete and correct |
|
7132
|
80 elseif (gnu == 0) |
|
3431
|
81 exit_1 = 2; # there are no zeros at all |
|
|
82 else |
|
|
83 mu = rho + sig; |
|
|
84 |
|
|
85 ## update system with Q |
|
7132
|
86 M = [Asys.a, Asys.b ]; |
|
|
87 [nn, mm] = size (Asys.b); |
|
3431
|
88 |
|
|
89 pp = rows(Asys.d); |
|
7132
|
90 Vm =[V, zeros(nn,mm); zeros(mm,nn), eye(mm)]; |
|
|
91 if (sig) |
|
3431
|
92 M = [M; Cb, Db]; |
|
7132
|
93 Vs =[V', zeros(nn,sig); zeros(sig,nn), eye(sig)]; |
|
3431
|
94 else |
|
|
95 Vs = V'; |
|
|
96 endif |
|
|
97 ## disp("zgreduce: before transform: M="); |
|
|
98 ## M |
|
|
99 ## Vs |
|
|
100 ## Vm |
|
|
101 |
|
|
102 M = Vs*M*Vm; |
|
|
103 |
|
|
104 ## disp("zgreduce: after transform: M="); |
|
|
105 ## M |
|
|
106 |
|
|
107 ## disp("debugging code:") |
|
|
108 ## Mtmp = [Asys.a Asys.b; Asys.c Asys.d] |
|
|
109 ## Vl = [V', zeros(nn,mm); zeros(mm,nn),Q] |
|
|
110 ## Vr =[V,zeros(nn,mm) ; zeros(mm,nn), eye(mm)]; |
|
|
111 ## Mtmpf = Vl*Mtmp*Vr |
|
|
112 |
|
|
113 idx = 1:gnu; |
|
|
114 jdx = nn + (1:mm); |
|
|
115 sdx = gnu + (1:mu); |
|
|
116 |
|
|
117 Asys.a = M(idx,idx); |
|
|
118 Asys.b = M(idx,jdx); |
|
|
119 Asys.c = M(sdx,idx); |
|
|
120 Asys.d = M(sdx,jdx); |
|
|
121 |
|
|
122 ## disp(["zgreduce: resulting system: nn =",num2str(nn)," mu=",num2str(mu)]) |
|
|
123 ## sysout(Asys) |
|
|
124 ## idx |
|
|
125 ## jdx |
|
|
126 ## sdx |
|
|
127 endif |
|
|
128 endif |
|
|
129 endwhile |
|
|
130 |
|
|
131 ## disp(["zgreduce: while loop done: exit_1=",num2str(exit_1)]); |
|
|
132 |
|
7132
|
133 if (exit_1 == 2) |
|
3431
|
134 ## there are no zeros at all! |
|
|
135 Asys.a = Asys.b = Asys.c = []; |
|
|
136 endif |
|
|
137 |
|
|
138 ## update dimensions |
|
7132
|
139 if (is_digital (Asys)) |
|
|
140 Asys.nz = rows (Asys.a); |
|
3431
|
141 else |
|
7132
|
142 Asys.n = rows (Asys.a); |
|
3431
|
143 endif |
|
|
144 |
|
|
145 retsys = Asys; |
|
|
146 |
|
|
147 endfunction |
