# Copyright (C) 1996 A. Scottedward Hodel 
#
# This file is part of Octave. 
#
# Octave is free software; you can redistribute it and/or modify it 
# under the terms of the GNU General Public License as published by the 
# Free Software Foundation; either version 2, or (at your option) any 
# later version. 
# 
# Octave is distributed in the hope that it will be useful, but WITHOUT 
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
# for more details.
# 
# You should have received a copy of the GNU General Public License 
# along with Octave; see the file COPYING.  If not, write to the Free 
# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. 
 
function retsys = zgreduce(Asys,meps)
# function retsys = zgreduce(Asys,meps)
# implementation of procedure REDUCE in (Emami-Naeini and Van Dooren, 
# Automatica, # 1982).
#
# used internally in tzero; minimal argument checking performed

#$Revision: 1.1.1.1 $
# SYS_INTERNAL accesses members of system data structure

is_digital(Asys);		# make sure it's pure digital/continuous

exit_1 = 0;			# exit_1 = 1 or 2 on exit of loop

if(Asys.n + Asys.nz == 0)
  exit_1 = 2;			# there are no finite zeros
endif

while (! exit_1)
  [Q,R,Pi] = qr(Asys.d);		# compress rows of D
  Asys.d = Q'*Asys.d;
  Asys.c = Q'*Asys.c;

  # check row norms of Asys.d
  [sig,tau] = zgrownorm(Asys.d,meps);

  #disp("=======================================")
  #disp(["zgreduce: meps=",num2str(meps), ", sig=",num2str(sig), ...
  #	", tau=",num2str(tau)])
  #sysout(Asys)

  if(tau == 0)
    exit_1 = 1;		# exit_1 - reduction complete and correct
  else
    Cb = Db = [];
    if(sig)
      Cb = Asys.c(1:sig,:);
      Db = Asys.d(1:sig,:);
    endif
    Ct =Asys.c(sig+(1:tau),:);

    # compress columns of Ct
    [pp,nn] = size(Ct);
    rvec = nn:-1:1;
    [V,Sj,Pi] = qr(Ct');
    V = V(:,rvec);
    [rho,gnu] = zgrownorm(Sj,meps);

    #disp(["zgreduce: rho=",num2str(rho),", gnu=",num2str(gnu)])
    #Cb
    #Db
    #Ct
    #Sj'

    if(rho == 0)
      exit_1 = 1;	# exit_1 - reduction complete and correct
    elseif(gnu == 0)
      exit_1 = 2;	# there are no zeros at all
    else
      mu = rho + sig;

      # update system with Q
      M = [Asys.a , Asys.b ];
      [nn,mm] = size(Asys.b);

      pp = rows(Asys.d);
      Vm =[V,zeros(nn,mm) ; zeros(mm,nn), eye(mm)];
      if(sig)
        M = [M; Cb, Db];
        Vs =[V',zeros(nn,sig) ; zeros(sig,nn), eye(sig)];
      else
        Vs = V';
      endif
      #disp("zgreduce: before transform: M=");
      #M
      #Vs   
      #Vm

      M = Vs*M*Vm;
      
      #disp("zgreduce: after transform: M=");
      #M

      #disp("debugging code:")
      #Mtmp = [Asys.a Asys.b; Asys.c Asys.d]
      #Vl = [V', zeros(nn,mm); zeros(mm,nn),Q]
      #Vr =[V,zeros(nn,mm) ; zeros(mm,nn), eye(mm)];
      #Mtmpf = Vl*Mtmp*Vr

      idx = 1:gnu;
      jdx = nn + (1:mm);
      sdx = gnu + (1:mu);

      Asys.a = M(idx,idx);
      Asys.b = M(idx,jdx);
      Asys.c = M(sdx,idx);
      Asys.d = M(sdx,jdx);

      #disp(["zgreduce: resulting system: nn =",num2str(nn)," mu=",num2str(mu)])
      #sysout(Asys)
      #idx
      #jdx
      #sdx
    endif
  endif
endwhile

#disp(["zgreduce: while loop done: exit_1=",num2str(exit_1)]);

if(exit_1 == 2)
  # there are no zeros at all!
  Asys.a = Asys.b = Asys.c = [];
endif

# update dimensions
if(is_digital(Asys))
  Asys.nz = rows(Asys.a);
else
  Asys.n = rows(Asys.a);
endif

retsys = Asys;

endfunction