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