home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts / control / zgreduce.m < prev    next >
Encoding:
Text File  |  1999-12-15  |  3.5 KB  |  143 lines
  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. 
  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
  24.  
  25. function retsys = zgreduce (Asys, meps)
  26.  
  27.   ## SYS_INTERNAL accesses members of system data structure
  28.  
  29.   is_digital(Asys);        # make sure it's pure digital/continuous
  30.  
  31.   exit_1 = 0;            # exit_1 = 1 or 2 on exit of loop
  32.  
  33.   if(Asys.n + Asys.nz == 0)
  34.     exit_1 = 2;            # there are no finite zeros
  35.   endif
  36.  
  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;
  41.  
  42.     ## check row norms of Asys.d
  43.     [sig,tau] = zgrownorm(Asys.d,meps);
  44.  
  45.     ## disp("=======================================")
  46.     ## disp(["zgreduce: meps=",num2str(meps), ", sig=",num2str(sig), ...
  47.     ##     ", tau=",num2str(tau)])
  48.     ## sysout(Asys)
  49.  
  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),:);
  59.  
  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);
  66.  
  67.       ## disp(["zgreduce: rho=",num2str(rho),", gnu=",num2str(gnu)])
  68.       ## Cb
  69.       ## Db
  70.       ## Ct
  71.       ## Sj'
  72.  
  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;
  79.  
  80.     ## update system with Q
  81.     M = [Asys.a , Asys.b ];
  82.     [nn,mm] = size(Asys.b);
  83.  
  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
  96.  
  97.     M = Vs*M*Vm;
  98.  
  99.     ## disp("zgreduce: after transform: M=");
  100.     ## M
  101.  
  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
  107.  
  108.     idx = 1:gnu;
  109.     jdx = nn + (1:mm);
  110.     sdx = gnu + (1:mu);
  111.  
  112.     Asys.a = M(idx,idx);
  113.     Asys.b = M(idx,jdx);
  114.     Asys.c = M(sdx,idx);
  115.     Asys.d = M(sdx,jdx);
  116.  
  117.     ## disp(["zgreduce: resulting system: nn =",num2str(nn)," mu=",num2str(mu)])
  118.     ## sysout(Asys)
  119.     ## idx
  120.     ## jdx
  121.     ## sdx
  122.       endif
  123.     endif
  124.   endwhile
  125.  
  126.   ## disp(["zgreduce: while loop done: exit_1=",num2str(exit_1)]);
  127.  
  128.   if(exit_1 == 2)
  129.     ## there are no zeros at all!
  130.     Asys.a = Asys.b = Asys.c = [];
  131.   endif
  132.  
  133.   ## update dimensions
  134.   if(is_digital(Asys))
  135.     Asys.nz = rows(Asys.a);
  136.   else
  137.     Asys.n = rows(Asys.a);
  138.   endif
  139.  
  140.   retsys = Asys;
  141.  
  142. endfunction
  143.