new file mode 100644
--- /dev/null
+++ b/scripts/control/zgreduce.m
@@ -0,0 +1,142 @@
+# 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