home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / zgreduce.m < prev    next >
Text File  |  1999-04-29  |  4KB  |  142 lines
  1. # Copyright (C) 1996 A. Scottedward Hodel 
  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, 675 Mass Ave, Cambridge, MA 02139, USA. 
  18.  
  19. function retsys = zgreduce(Asys,meps)
  20. # function retsys = zgreduce(Asys,meps)
  21. # implementation of procedure REDUCE in (Emami-Naeini and Van Dooren, 
  22. # Automatica, # 1982).
  23. #
  24. # used internally in tzero; minimal argument checking performed
  25.  
  26. # SYS_INTERNAL accesses members of system data structure
  27.  
  28. is_digit(Asys);        # make sure it's pure digital/continuous
  29.  
  30. exit_1 = 0;            # exit_1 = 1 or 2 on exit of loop
  31.  
  32. if(Asys.n + Asys.nz == 0)
  33.   exit_1 = 2;            # there are no finite zeros
  34. endif
  35.  
  36. while (! exit_1)
  37.   [Q,R,Pi] = qr(Asys.d);        # compress rows of D
  38.   Asys.d = Q'*Asys.d;
  39.   Asys.c = Q'*Asys.c;
  40.  
  41.   # check row norms of Asys.d
  42.   [sig,tau] = zgrownor(Asys.d,meps);
  43.  
  44.   #disp("=======================================")
  45.   #disp(["zgreduce: meps=",num2str(meps), ", sig=",num2str(sig), ...
  46.   #    ", tau=",num2str(tau)])
  47.   #sysout(Asys)
  48.  
  49.   if(tau == 0)
  50.     exit_1 = 1;        # exit_1 - reduction complete and correct
  51.   else
  52.     Cb = Db = [];
  53.     if(sig)
  54.       Cb = Asys.c(1:sig,:);
  55.       Db = Asys.d(1:sig,:);
  56.     endif
  57.     Ct =Asys.c(sig+(1:tau),:);
  58.  
  59.     # compress columns of Ct
  60.     [pp,nn] = size(Ct);
  61.     rvec = nn:-1:1;
  62.     [V,Sj,Pi] = qr(Ct');
  63.     V = V(:,rvec);
  64.     [rho,gnu] = zgrownor(Sj,meps);
  65.  
  66.     #disp(["zgreduce: rho=",num2str(rho),", gnu=",num2str(gnu)])
  67.     #Cb
  68.     #Db
  69.     #Ct
  70.     #Sj'
  71.  
  72.     if(rho == 0)
  73.       exit_1 = 1;    # exit_1 - reduction complete and correct
  74.     elseif(gnu == 0)
  75.       exit_1 = 2;    # there are no zeros at all
  76.     else
  77.       mu = rho + sig;
  78.  
  79.       # update system with Q
  80.       M = [Asys.a , Asys.b ];
  81.       [nn,mm] = size(Asys.b);
  82.  
  83.       pp = rows(Asys.d);
  84.       Vm =[V,zeros(nn,mm) ; zeros(mm,nn), eye(mm)];
  85.       if(sig)
  86.         M = [M; Cb, Db];
  87.         Vs =[V',zeros(nn,sig) ; zeros(sig,nn), eye(sig)];
  88.       else
  89.         Vs = V';
  90.       endif
  91.       #disp("zgreduce: before transform: M=");
  92.       #M
  93.       #Vs   
  94.       #Vm
  95.  
  96.       M = Vs*M*Vm;
  97.       
  98.       #disp("zgreduce: after transform: M=");
  99.       #M
  100.  
  101.       #disp("debugging code:")
  102.       #Mtmp = [Asys.a Asys.b; Asys.c Asys.d]
  103.       #Vl = [V', zeros(nn,mm); zeros(mm,nn),Q]
  104.       #Vr =[V,zeros(nn,mm) ; zeros(mm,nn), eye(mm)];
  105.       #Mtmpf = Vl*Mtmp*Vr
  106.  
  107.       idx = 1:gnu;
  108.       jdx = nn + (1:mm);
  109.       sdx = gnu + (1:mu);
  110.  
  111.       Asys.a = M(idx,idx);
  112.       Asys.b = M(idx,jdx);
  113.       Asys.c = M(sdx,idx);
  114.       Asys.d = M(sdx,jdx);
  115.  
  116.       #disp(["zgreduce: resulting system: nn =",num2str(nn)," mu=",num2str(mu)])
  117.       #sysout(Asys)
  118.       #idx
  119.       #jdx
  120.       #sdx
  121.     endif
  122.   endif
  123. endwhile
  124.  
  125. #disp(["zgreduce: while loop done: exit_1=",num2str(exit_1)]);
  126.  
  127. if(exit_1 == 2)
  128.   # there are no zeros at all!
  129.   Asys.a = Asys.b = Asys.c = [];
  130. endif
  131.  
  132. # update dimensions
  133. if(is_digit(Asys))
  134.   Asys.nz = rows(Asys.a);
  135. else
  136.   Asys.n = rows(Asys.a);
  137. endif
  138.  
  139. retsys = Asys;
  140.  
  141. endfunction
  142.