## Copyright (C) 1996 Auburn University.  All rights reserved.
##
## 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. 

## -*- texinfo -*-
## @deftypefn {Function File } { retsys = } zgreduce(@var{Asys},@var{meps})
## Implementation of procedure REDUCE in (Emami-Naeini and Van Dooren, 
## Automatica, # 1982).
## @end deftypefn
 
function retsys = zgreduce (Asys, meps)

  ## 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