|
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 |
|
|
28 ## SYS_INTERNAL accesses members of system data structure |
|
|
29 |
|
|
30 is_digital(Asys); # make sure it's pure digital/continuous |
|
|
31 |
|
|
32 exit_1 = 0; # exit_1 = 1 or 2 on exit of loop |
|
|
33 |
|
|
34 if(Asys.n + Asys.nz == 0) |
|
|
35 exit_1 = 2; # there are no finite zeros |
|
|
36 endif |
|
|
37 |
|
|
38 while (! exit_1) |
|
|
39 [Q,R,Pi] = qr(Asys.d); # compress rows of D |
|
|
40 Asys.d = Q'*Asys.d; |
|
|
41 Asys.c = Q'*Asys.c; |
|
|
42 |
|
|
43 ## check row norms of Asys.d |
|
|
44 [sig,tau] = zgrownorm(Asys.d,meps); |
|
|
45 |
|
|
46 ## disp("=======================================") |
|
|
47 ## disp(["zgreduce: meps=",num2str(meps), ", sig=",num2str(sig), ... |
|
|
48 ## ", tau=",num2str(tau)]) |
|
|
49 ## sysout(Asys) |
|
|
50 |
|
|
51 if(tau == 0) |
|
|
52 exit_1 = 1; # exit_1 - reduction complete and correct |
|
|
53 else |
|
|
54 Cb = Db = []; |
|
|
55 if(sig) |
|
|
56 Cb = Asys.c(1:sig,:); |
|
|
57 Db = Asys.d(1:sig,:); |
|
|
58 endif |
|
|
59 Ct =Asys.c(sig+(1:tau),:); |
|
|
60 |
|
|
61 ## compress columns of Ct |
|
|
62 [pp,nn] = size(Ct); |
|
|
63 rvec = nn:-1:1; |
|
|
64 [V,Sj,Pi] = qr(Ct'); |
|
|
65 V = V(:,rvec); |
|
|
66 [rho,gnu] = zgrownorm(Sj,meps); |
|
|
67 |
|
|
68 ## disp(["zgreduce: rho=",num2str(rho),", gnu=",num2str(gnu)]) |
|
|
69 ## Cb |
|
|
70 ## Db |
|
|
71 ## Ct |
|
|
72 ## Sj' |
|
|
73 |
|
|
74 if(rho == 0) |
|
|
75 exit_1 = 1; # exit_1 - reduction complete and correct |
|
|
76 elseif(gnu == 0) |
|
|
77 exit_1 = 2; # there are no zeros at all |
|
|
78 else |
|
|
79 mu = rho + sig; |
|
|
80 |
|
|
81 ## update system with Q |
|
|
82 M = [Asys.a , Asys.b ]; |
|
|
83 [nn,mm] = size(Asys.b); |
|
|
84 |
|
|
85 pp = rows(Asys.d); |
|
|
86 Vm =[V,zeros(nn,mm) ; zeros(mm,nn), eye(mm)]; |
|
|
87 if(sig) |
|
|
88 M = [M; Cb, Db]; |
|
|
89 Vs =[V',zeros(nn,sig) ; zeros(sig,nn), eye(sig)]; |
|
|
90 else |
|
|
91 Vs = V'; |
|
|
92 endif |
|
|
93 ## disp("zgreduce: before transform: M="); |
|
|
94 ## M |
|
|
95 ## Vs |
|
|
96 ## Vm |
|
|
97 |
|
|
98 M = Vs*M*Vm; |
|
|
99 |
|
|
100 ## disp("zgreduce: after transform: M="); |
|
|
101 ## M |
|
|
102 |
|
|
103 ## disp("debugging code:") |
|
|
104 ## Mtmp = [Asys.a Asys.b; Asys.c Asys.d] |
|
|
105 ## Vl = [V', zeros(nn,mm); zeros(mm,nn),Q] |
|
|
106 ## Vr =[V,zeros(nn,mm) ; zeros(mm,nn), eye(mm)]; |
|
|
107 ## Mtmpf = Vl*Mtmp*Vr |
|
|
108 |
|
|
109 idx = 1:gnu; |
|
|
110 jdx = nn + (1:mm); |
|
|
111 sdx = gnu + (1:mu); |
|
|
112 |
|
|
113 Asys.a = M(idx,idx); |
|
|
114 Asys.b = M(idx,jdx); |
|
|
115 Asys.c = M(sdx,idx); |
|
|
116 Asys.d = M(sdx,jdx); |
|
|
117 |
|
|
118 ## disp(["zgreduce: resulting system: nn =",num2str(nn)," mu=",num2str(mu)]) |
|
|
119 ## sysout(Asys) |
|
|
120 ## idx |
|
|
121 ## jdx |
|
|
122 ## sdx |
|
|
123 endif |
|
|
124 endif |
|
|
125 endwhile |
|
|
126 |
|
|
127 ## disp(["zgreduce: while loop done: exit_1=",num2str(exit_1)]); |
|
|
128 |
|
|
129 if(exit_1 == 2) |
|
|
130 ## there are no zeros at all! |
|
|
131 Asys.a = Asys.b = Asys.c = []; |
|
|
132 endif |
|
|
133 |
|
|
134 ## update dimensions |
|
|
135 if(is_digital(Asys)) |
|
|
136 Asys.nz = rows(Asys.a); |
|
|
137 else |
|
|
138 Asys.n = rows(Asys.a); |
|
|
139 endif |
|
|
140 |
|
|
141 retsys = Asys; |
|
|
142 |
|
|
143 endfunction |
