# HG changeset patch
# User jwe
# Date 1194538654 0
# Node ID b01db194c5267faed650f25125f62f28aa74c7b4
# Parent  a184bc985c3730539ed95271e205c8ec946ba08c
[project @ 2007-11-08 16:17:34 by jwe]

diff -r a184bc985c37 -r b01db194c526 scripts/ChangeLog
--- a/scripts/ChangeLog	Thu Nov 08 15:55:02 2007 +0000
+++ b/scripts/ChangeLog	Thu Nov 08 16:17:34 2007 +0000
@@ -3,7 +3,12 @@
 	* control/hinf/h2syn.m, control/hinf/hinf_ctr.m,
 	control/hinf/hinfnorm.m, control/hinf/hinfsyn.m,
 	control/hinf/hinfsyn_chk.m, control/hinf/is_dgkf.m,
-	control/hinf/wgt1o.m: Style fixes.
+	control/hinf/wgt1o.m, control/util/__outlist__.m,
+	control/util/__zgpbal__.m, control/util/axis2dlim.m,
+	control/util/prompt.m, control/util/sortcom.m,
+	control/util/zgfmul.m, control/util/zgfslv.m,
+	control/util/zginit.m, control/util/zgreduce.m,
+	control/util/zgrownorm.m, control/util/zgscal.m: Style fixes.
 
 2007-11-08  David Bateman  <dbateman@free.fr>
 
diff -r a184bc985c37 -r b01db194c526 scripts/control/util/__outlist__.m
--- a/scripts/control/util/__outlist__.m	Thu Nov 08 15:55:02 2007 +0000
+++ b/scripts/control/util/__outlist__.m	Thu Nov 08 16:17:34 2007 +0000
@@ -49,33 +49,33 @@
 
 function str_val = __outlist__ (name_list, tabchar, yd, ilist)
 
-  if( nargin < 1 | nargin > 4 )
+  if (nargin < 1 || nargin > 4)
     print_usage ();
   endif
 
-  m = length(name_list);
-  if(nargin < 4)
+  m = length (name_list);
+  if (nargin < 4)
     ilist = 1:m;
   endif
-  if(nargin ==1)
+  if (nargin == 1)
     tabchar = "";
   endif
 
-  if(nargin < 3)
-    yd = zeros(1,m);
-  elseif(isempty(yd))
-    yd = zeros(1,m);
+  if (nargin < 3)
+    yd = zeros (1, m);
+  elseif (isempty (yd))
+    yd = zeros (1, m);
   endif
 
   str_val = "";
-  dstr = {""," (discrete)"};
-  if((m >= 1) && (iscell(name_list)))
-    for ii=1:m
-	str_val = sprintf("%s%s%d: %s%s\n",str_val,tabchar, ilist(ii), ...
-			  name_list{ii},dstr{yd(ii)+1});
+  dstr = {"", " (discrete)"};
+  if (m >= 1 && iscell (name_list))
+    for ii = 1:m
+      str_val = sprintf ("%s%s%d: %s%s\n", str_val, tabchar, ilist(ii),
+			 name_list{ii}, dstr{yd(ii)+1});
     endfor
   else
-    str_val = sprintf("%sNone",tabchar);
+    str_val = sprintf ("%sNone", tabchar);
   endif
 
 endfunction
diff -r a184bc985c37 -r b01db194c526 scripts/control/util/__zgpbal__.m
--- a/scripts/control/util/__zgpbal__.m	Thu Nov 08 15:55:02 2007 +0000
+++ b/scripts/control/util/__zgpbal__.m	Thu Nov 08 16:17:34 2007 +0000
@@ -47,67 +47,67 @@
 
 function retsys = __zgpbal__ (Asys)
 
-  if( (nargin != 1) | (!isstruct(Asys)))
+  if (nargin != 1 || ! isstruct (Asys))
     print_usage ();
   endif
 
-  Asys = sysupdate(Asys,"ss");
-  [a,b,c,d] = sys2ss(Asys);
+  Asys = sysupdate (Asys, "ss");
+  [a, b, c, d] = sys2ss (Asys);
 
-  [nn,mm,pp] = abcddim(a,b,c,d);
+  [nn, mm, pp] = abcddim (a, b, c, d);
 
   np1 = nn+1;
   nmp = nn+mm+pp;
 
   ## set up log vector zz, incidence matrix ff
-  zz = zginit(a,b,c,d);
+  zz = zginit (a, b, c, d);
 
   ## disp("__zgpbal__: zginit returns")
   ## zz
   ## disp("/__zgpbal__")
 
-  if (norm(zz))
+  if (norm (zz))
     ## generalized conjugate gradient approach
-    xx = zgscal(a,b,c,d,zz,nn,mm,pp);
+    xx = zgscal (a, b, c, d, zz, nn, mm, pp);
 
-    for i=1:nmp
-      xx(i) = floor(xx(i)+0.5);
+    for i = 1:nmp
+      xx(i) = floor (xx(i)+0.5);
       xx(i) = 2.0^xx(i);
     endfor
 
     ## now scale a
     ## block 1: a = sigma a inv(sigma)
-    for i=1:nn
+    for i = 1:nn
       a(i,1:nn) = a(i,1:nn)*xx(i);
       a(1:nn,i) = a(1:nn,i)/xx(i);
     endfor
     ## block 2: b= sigma a phi
-    for j=1:mm
+    for j = 1:mm
       j1 = j+nn;
       b(1:nn,j) = b(1:nn,j)*xx(j1);
     endfor
-    for i=1:nn
+    for i = 1:nn
       b(i,1:mm) = b(i,1:mm)*xx(i);
     endfor
-    for i=1:pp
+    for i = 1:pp
       i1 = i+nn+mm;
       ## block 3: c = psi C inv(sigma)
       c(i,1:nn) = c(i,1:nn)*xx(i1);
     endfor
-    for j=1:nn
+    for j = 1:nn
       c(1:pp,j) = c(1:pp,j)/xx(j);
     endfor
     ## block 4: d = psi D phi
-    for j=1:mm
+    for j = 1:mm
       j1 = j+nn;
       d(1:pp,j) = d(1:pp,j)*xx(j1);
     endfor
-    for i=1:pp
+    for i = 1:pp
       i1 = i + nn + mm;
       d(i,1:mm) = d(i,1:mm)*xx(i1);
     endfor
   endif
 
-  retsys = ss(a,b,c,d);
+  retsys = ss (a, b, c, d);
+
 endfunction
-
diff -r a184bc985c37 -r b01db194c526 scripts/control/util/axis2dlim.m
--- a/scripts/control/util/axis2dlim.m	Thu Nov 08 15:55:02 2007 +0000
+++ b/scripts/control/util/axis2dlim.m	Thu Nov 08 16:17:34 2007 +0000
@@ -49,24 +49,24 @@
   delv = (maxv-minv)/2;      # breadth of the plot
   midv = (minv + maxv)/2;    # midpoint of the plot
   axmid = [midv(1), midv(1), midv(2), midv(2)];
-  axdel = [-0.1, 0.1,-0.1,0.1];   # default plot width (if less than 2-d data)
-  if(max(delv) == 0)
-    if(midv(1) != 0)
-      axdel(1:2) = [-0.1*midv(1),0.1*midv(1)];
+  axdel = [-0.1, 0.1, -0.1, 0.1];   # default plot width (if less than 2-d data)
+  if (max (delv) == 0)
+    if (midv(1) != 0)
+      axdel(1:2) = [-0.1*midv(1), 0.1*midv(1)];
     endif
-    if(midv(2) != 0)
-      axdel(3:4) = [-0.1*midv(2),0.1*midv(2)];
+    if (midv(2) != 0)
+      axdel(3:4) = [-0.1*midv(2), 0.1*midv(2)];
     endif
   else
     ## they're at least one-dimensional
     tolv = max(1e-8, 1e-8*abs(midv));
-    if(abs(delv(1)) >= tolv(1))
+    if (abs (delv(1)) >= tolv(1))
       axdel(1:2) = 1.1*[-delv(1),delv(1)];
     endif
-    if(abs(delv(2)) >= tolv(2))
+    if (abs (delv(2)) >= tolv(2))
       axdel(3:4) = 1.1*[-delv(2),delv(2)];
     endif
   endif
   axvec = axmid + axdel;
+
 endfunction
-
diff -r a184bc985c37 -r b01db194c526 scripts/control/util/prompt.m
--- a/scripts/control/util/prompt.m	Thu Nov 08 15:55:02 2007 +0000
+++ b/scripts/control/util/prompt.m	Thu Nov 08 16:17:34 2007 +0000
@@ -41,7 +41,7 @@
     print_usage ();
   elseif (nargin == 0)
     str = "\n ---- Press a key to continue ---";
-  elseif (! ischar (str) )
+  elseif (! ischar (str))
     error ("prompt: input must be a string");
   endif
 
diff -r a184bc985c37 -r b01db194c526 scripts/control/util/sortcom.m
--- a/scripts/control/util/sortcom.m	Thu Nov 08 15:55:02 2007 +0000
+++ b/scripts/control/util/sortcom.m	Thu Nov 08 16:17:34 2007 +0000
@@ -53,42 +53,47 @@
 
 function [yy, idx] = sortcom (xx, opt)
 
-  if( nargin < 1 | nargin > 2 )
+  if (nargin < 1 || nargin > 2)
      print_usage ();
-  elseif( !(isvector(xx) | isempty(xx) ))
-    error("sortcom: first argument must be a vector");
+  elseif (! (isvector (xx) || isempty (xx)))
+    error ("sortcom: first argument must be a vector");
   endif
 
-  if(nargin == 1)         opt = "re";
+  if (nargin == 1)
+    opt = "re";
   else
-    if (!ischar(opt))
-      error("sortcom: second argument must be a string");
+    if (! ischar (opt))
+      error ("sortcom: second argument must be a string");
     endif
   endif
 
-  if(isempty(xx))
+  if (isempty (xx))
     yy = idx = [];
   else
-    if(strcmp(opt,"re"))        datavec = real(xx);
-    elseif(strcmp(opt,"im"))    datavec = imag(xx);
-    elseif(strcmp(opt,"mag"))   datavec = abs(xx);
-    else                        error(["sortcom: invalid option = ", opt])
+    if (strcmp (opt, "re"))
+      datavec = real (xx);
+    elseif (strcmp (opt, "im"))
+      datavec = imag (xx);
+    elseif (strcmp (opt, "mag"))
+      datavec = abs (xx);
+    else
+      error ("sortcom: invalid option = %s", opt);
     endif
 
-    [datavec,idx] = sort(datavec);
+    [datavec, idx] = sort (datavec);
     yy= xx(idx);
 
-    if(strcmp(opt,"re") | strcmp(opt,"mag"))
+    if (strcmp (opt, "re") || strcmp (opt, "mag"))
       ## sort so that complex conjugate pairs appear together
 
-      ddiff = diff(datavec);
-      zidx = find(ddiff == 0);
+      ddiff = diff (datavec);
+      zidx = find (ddiff == 0);
 
       ## sort common datavec values
-      if(!isempty(zidx))
-        for iv=create_set(datavec(zidx))
-          vidx = find(datavec == iv);
-          [vals,imidx] = sort(imag(yy(vidx)));
+      if (! isempty (zidx))
+        for iv = create_set (datavec(zidx))
+          vidx = find (datavec == iv);
+          [vals, imidx] = sort (imag (yy(vidx)));
           yy(vidx)  = yy(vidx(imidx));
           idx(vidx) = idx(vidx(imidx));
         endfor
diff -r a184bc985c37 -r b01db194c526 scripts/control/util/zgfmul.m
--- a/scripts/control/util/zgfmul.m	Thu Nov 08 15:55:02 2007 +0000
+++ b/scripts/control/util/zgfmul.m	Thu Nov 08 16:17:34 2007 +0000
@@ -37,48 +37,67 @@
     print_usage ();
   endif 
 
-  [n,m] = size(b);
-  [p,m1] = size(c);
+  [n, m] = size (b);
+  [p, m1] = size (c);
   nm = n+m;
-  y = zeros(nm+p,1);
+  y = zeros (nm+p, 1);
 
   ## construct F column by column
-  for jj=1:n
-    Fj = zeros(nm+p,1);
+  for jj = 1:n
+    Fj = zeros (nm+p, 1);
 
     ## rows 1:n: F1
-    aridx = complement(jj,find(a(jj,:) != 0));
-    acidx = complement(jj,find(a(:,jj) != 0));
-    bidx = find(b(jj,:) != 0);
-    cidx = find(c(:,jj) != 0);
+    aridx = complement (jj, find (a(jj,:) != 0));
+    acidx = complement (jj, find (a(:,jj) != 0));
+    bidx = find (b(jj,:) != 0);
+    cidx = find (c(:,jj) != 0);
 
     Fj(aridx) = Fj(aridx) - 1;      # off diagonal entries of F1
     Fj(acidx) = Fj(acidx) - 1;
     ## diagonal entry of F1
-    Fj(jj) = length(aridx)+length(acidx) + length(bidx) + length(cidx);
+    Fj(jj) = length (aridx) + length (acidx) + length (bidx) + length (cidx);
 
-    if(!isempty(bidx)) Fj(n+bidx) = 1;     endif # B' incidence
-    if(!isempty(cidx)) Fj(n+m+cidx) = -1;  endif # -C incidence
+    ## B' incidence
+    if (! isempty (bidx))
+      Fj(n+bidx) = 1;
+    endif
+
+    ## -C incidence
+    if (! isempty (cidx))
+      Fj(n+m+cidx) = -1;
+    endif
     y = y + x(jj)*Fj;   # multiply by corresponding entry of x
   endfor
 
-  for jj=1:m
-    Fj = zeros(nm+p,1);
-    bidx = find(b(:,jj) != 0);
-    if(!isempty(bidx)) Fj(bidx) = 1; endif     # B incidence
-    didx = find(d(:,jj) != 0);
-    if(!isempty(didx)) Fj(n+m+didx) = 1; endif # D incidence
+  for jj = 1:m
+    Fj = zeros (nm+p, 1);
+    bidx = find (b(:,jj) != 0);
+    ## B incidence
+    if (! isempty (bidx))
+      Fj(bidx) = 1;
+    endif
+    didx = find (d(:,jj) != 0);
+    ## D incidence
+    if (! isempty (didx))
+      Fj(n+m+didx) = 1;
+    endif
     Fj(n+jj) = length(bidx) + length(didx);         # F2 is diagonal
     y = y + x(n+jj)*Fj;   # multiply by corresponding entry of x
   endfor
 
-  for jj=1:p
-    Fj = zeros(nm+p,1);
-    cidx = find(c(jj,:) != 0);
-    if(!isempty(cidx)) Fj(cidx) = -1; endif  # -C' incidence
+  for jj = 1:p
+    Fj = zeros (nm+p, 1);
+    cidx = find (c(jj,:) != 0);
+    ## -C' incidence
+    if (! isempty (cidx))
+      Fj(cidx) = -1;
+    endif
     didx = find(d(jj,:) != 0);
-    if(!isempty(didx)) Fj(n+didx) = 1;  endif # D' incidence
-    Fj(n+m+jj) = length(cidx) + length(didx);     # F2 is diagonal
+    ## D' incidence
+    if (! isempty (didx))
+      Fj(n+didx) = 1;
+    endif
+    Fj(n+m+jj) = length (cidx) + length (didx);     # F2 is diagonal
     y = y + x(n+m+jj)*Fj;   # multiply by corresponding entry of x
   endfor
 
diff -r a184bc985c37 -r b01db194c526 scripts/control/util/zgfslv.m
--- a/scripts/control/util/zgfslv.m	Thu Nov 08 15:55:02 2007 +0000
+++ b/scripts/control/util/zgfslv.m	Thu Nov 08 16:17:34 2007 +0000
@@ -32,41 +32,51 @@
   endif
 
   nmp = n+m+p;
-  gam1 = (2*n)+m+p;    gam2 = n+p;     gam3 = n+m;
+  gam1 = (2*n)+m+p;
+  gam2 = n+p;
+  gam3 = n+m;
 
-  G1 = givens(sqrt(m),-sqrt(p))';
-  G2 = givens(m+p,sqrt(n*(m+p)))';
+  G1 = givens (sqrt (m), -sqrt (p))';
+  G2 = givens (m+p, sqrt (n*(m+p)))';
 
   x = b;
 
   ## 1) U1 e^n = sqrt(n)e_1^n
   ## 2) U2 e^m = sqrt(m)e_1^m
   ## 3) U3 e^p = sqrt(p)e_1^p
-  xdx1 = 1:n; xdx2 = n+(1:m); xdx3 = n+m+(1:p);
-  x(xdx1,1) = zgshsr(x(xdx1,1));
-  x(xdx2,1) = zgshsr(x(xdx2,1));
-  x(xdx3,1) = zgshsr(x(xdx3,1));
+  xdx1 = 1:n;
+  xdx2 = n+(1:m);
+  xdx3 = n+m+(1:p);
+
+  x(xdx1,1) = zgshsr (x(xdx1,1));
+  x(xdx2,1) = zgshsr (x(xdx2,1));
+  x(xdx3,1) = zgshsr (x(xdx3,1));
 
   ## 4) Givens rotations to reduce stray non-zero elements
-  idx1 = [n+1,n+m+1];     idx2 = [1,n+1];
+  idx1 = [n+1, n+m+1];
+  idx2 = [1, n+1];
+
   x(idx1) = G1'*x(idx1);
   x(idx2) = G2'*x(idx2);
 
   ## 6) Scale x, then back-transform to get x
-  en = ones(n,1);  em = ones(m,1);   ep = ones(p,1);
-  lam = [gam1*en;gam2*em;gam3*ep];
+  en = ones (n, 1);
+  em = ones (m, 1);
+  ep = ones (p, 1);
+  lam = [gam1*en; gam2*em; gam3*ep];
   lam(1) = n+m+p;
   lam(n+1) = 1;       # dummy value to avoid divide by zero
-  lam(n+m+1)=n+m+p;
+  lam(n+m+1) = n+m+p;
 
-  x = x ./ lam;       x(n+1) = 0;  # minimum norm solution
+  x = x ./ lam;
+  x(n+1) = 0;  # minimum norm solution
 
   ## back transform now.
   x(idx2) = G2*x(idx2);
   x(idx1) = G1*x(idx1);
-  x(xdx3,1) = zgshsr(x(xdx3,1));
-  x(xdx2,1) = zgshsr(x(xdx2,1));
-  x(xdx1,1) = zgshsr(x(xdx1,1));
+  x(xdx3,1) = zgshsr (x(xdx3,1));
+  x(xdx2,1) = zgshsr (x(xdx2,1));
+  x(xdx1,1) = zgshsr (x(xdx1,1));
 
 endfunction
 
diff -r a184bc985c37 -r b01db194c526 scripts/control/util/zginit.m
--- a/scripts/control/util/zginit.m	Thu Nov 08 15:55:02 2007 +0000
+++ b/scripts/control/util/zginit.m	Thu Nov 08 16:17:34 2007 +0000
@@ -38,58 +38,61 @@
     print_usage ();
   endif
 
-  [nn,mm] = size(b);
-  [pp,mm] = size(d);
+  [nn, mm] = size (b);
+  [pp, mm] = size (d);
 
   nmp = nn+mm+pp;
 
   ## set up log vector zz
-  zz = zeros(nmp,1);
+  zz = zeros (nmp, 1);
 
   ## zz part 1:
-  for i=1:nn
+  for i = 1:nn
     ## nonzero off diagonal entries of a
-    if(nn > 1)
-      nidx = complement(i,1:nn);
-      a_row_i = a(i,nidx);                 a_col_i = a(nidx,i);
-      arnz = a_row_i(find(a_row_i != 0));  acnz = a_col_i(find(a_col_i != 0));
+    if (nn > 1)
+      nidx = complement (i, 1:nn);
+      a_row_i = a(i,nidx);
+      a_col_i = a(nidx,i);
+      arnz = a_row_i(find (a_row_i != 0));
+      acnz = a_col_i(find (a_col_i != 0));
     else
       arnz = acnz = [];
     endif
 
     ## row of b
-    bidx = find(b(i,:) != 0);
+    bidx = find (b(i,:) != 0);
     b_row_i = b(i,bidx);
 
     ## column of c
-    cidx = find(c(:,i) != 0);
+    cidx = find (c(:,i) != 0);
     c_col_i = c(cidx,i);
 
     ## sum the entries
-    zz(i) = sum(log(abs(acnz))) - sum(log(abs(arnz))) ...
-            - sum(log(abs(b_row_i))) + sum(log(abs(c_col_i)));
+    zz(i) = sum (log (abs (acnz))) - sum (log (abs (arnz))) ...
+        - sum (log (abs (b_row_i))) + sum (log (abs (c_col_i)));
   endfor
 
   ## zz part 2:
-  bd = [b;d];
-  for i=1:mm
+  bd = [b; d];
+  for i = 1:mm
     i1 = i+nn;
 
     ## column of [b;d]
-    bdidx = find(bd(:,i) != 0);
+    bdidx = find (bd(:,i) != 0);
     bd_col_i = bd(bdidx,i);
-    zz(i1) = sum(log(abs(bd_col_i)));
+    zz(i1) = sum (log (abs(bd_col_i)));
   endfor
 
   ## zz part 3:
   cd = [c, d];
-  for i=1:pp
+  for i = 1:pp
     i1 = i+nn+mm;
-    cdidx = find(cd(i,:) != 0);
+    cdidx = find (cd(i,:) != 0);
     cd_row_i = cd(i,cdidx);
-    zz(i1) = -sum(log(abs(cd_row_i)));
+    zz(i1) = -sum (log (abs (cd_row_i)));
   endfor
 
   ## now set zz as log base 2
-  zz = zz*(1/log(2));
+  zz *= (1 / log (2));
+
 endfunction
diff -r a184bc985c37 -r b01db194c526 scripts/control/util/zgreduce.m
--- a/scripts/control/util/zgreduce.m	Thu Nov 08 15:55:02 2007 +0000
+++ b/scripts/control/util/zgreduce.m	Thu Nov 08 16:17:34 2007 +0000
@@ -31,43 +31,43 @@
 
   ## SYS_INTERNAL accesses members of system data structure
 
-  is_digital(Asys);             # make sure it's pure digital/continuous
+  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)
+  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
+    [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);
+    [sig, tau] = zgrownorm (Asys.d, meps);
 
     ## disp("=======================================")
     ## disp(["zgreduce: meps=",num2str(meps), ", sig=",num2str(sig), ...
     ##   ", tau=",num2str(tau)])
     ## sysout(Asys)
 
-    if(tau == 0)
+    if (tau == 0)
       exit_1 = 1;               # exit_1 - reduction complete and correct
     else
       Cb = Db = [];
-      if(sig)
+      if (sig)
         Cb = Asys.c(1:sig,:);
         Db = Asys.d(1:sig,:);
       endif
-      Ct =Asys.c(sig+(1:tau),:);
+      Ct = Asys.c(sig+(1:tau),:);
 
       ## compress columns of Ct
-      [pp,nn] = size(Ct);
+      [pp, nn] = size (Ct);
       rvec = nn:-1:1;
-      [V,Sj,Pi] = qr(Ct');
+      [V, Sj, Pi] = qr (Ct');
       V = V(:,rvec);
-      [rho,gnu] = zgrownorm(Sj,meps);
+      [rho, gnu] = zgrownorm (Sj, meps);
 
       ## disp(["zgreduce: rho=",num2str(rho),", gnu=",num2str(gnu)])
       ## Cb
@@ -75,22 +75,22 @@
       ## Ct
       ## Sj'
 
-      if(rho == 0)
+      if (rho == 0)
         exit_1 = 1;     # exit_1 - reduction complete and correct
-      elseif(gnu == 0)
+      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);
+        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)
+        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)];
+          Vs =[V', zeros(nn,sig); zeros(sig,nn), eye(sig)];
         else
           Vs = V';
         endif
@@ -130,16 +130,16 @@
 
   ## disp(["zgreduce: while loop done: exit_1=",num2str(exit_1)]);
 
-  if(exit_1 == 2)
+  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);
+  if (is_digital (Asys))
+    Asys.nz = rows (Asys.a);
   else
-    Asys.n = rows(Asys.a);
+    Asys.n = rows (Asys.a);
   endif
 
   retsys = Asys;
diff -r a184bc985c37 -r b01db194c526 scripts/control/util/zgrownorm.m
--- a/scripts/control/util/zgrownorm.m	Thu Nov 08 15:55:02 2007 +0000
+++ b/scripts/control/util/zgrownorm.m	Thu Nov 08 16:17:34 2007 +0000
@@ -38,4 +38,3 @@
   tau = sum (rownorm <= meps);
 
 endfunction
-
diff -r a184bc985c37 -r b01db194c526 scripts/control/util/zgscal.m
--- a/scripts/control/util/zgscal.m	Thu Nov 08 15:55:02 2007 +0000
+++ b/scripts/control/util/zgscal.m	Thu Nov 08 16:17:34 2007 +0000
@@ -44,27 +44,27 @@
   nmp = n+m+p;
 
   ## x_0 = x_{-1} = 0, r_0 = z
-  x = zeros(nmp,1);
+  x = zeros (nmp, 1);
   xk1 = x;
   xk2 = x;
   rk1 = z;
   k = 0;
 
   ## construct balancing least squares problem
-  F = eye(nmp);
-  for kk=1:nmp
-    F(1:nmp,kk) = zgfmul(a,b,c,d,F(:,kk));
+  F = eye (nmp);
+  for kk = 1:nmp
+    F(1:nmp,kk) = zgfmul (a, b, c, d, F(:,kk));
   endfor
 
-  [U,H,k1] = krylov(F,z,nmp,1e-12,1);
-  if(!issquare(H))
-    if(columns(H) != k1)
-      error("zgscal(tzero): k1=%d, columns(H)=%d",k1,columns(H));
-    elseif(rows(H) != k1+1)
-      error("zgscal: k1=%d, rows(H) = %d",k1,rows(H));
-    elseif ( norm(H(k1+1,:)) > 1e-12*norm(H,"inf") )
+  [U, H, k1] = krylov (F, z, nmp, 1e-12, 1);
+  if (! issquare (H))
+    if (columns (H) != k1)
+      error ("zgscal(tzero): k1=%d, columns(H)=%d", k1, columns (H));
+    elseif (rows (H) != k1+1)
+      error ("zgscal: k1=%d, rows(H) = %d", k1, rows (H));
+    elseif (norm (H(k1+1,:)) > 1e-12*norm (H, "inf"))
       zgscal_last_row_of_H = H(k1+1,:)
-      error("zgscal: last row of H nonzero (norm(H)=%e)",norm(H,"inf"))
+      error ("zgscal: last row of H nonzero (norm(H)=%e)", norm (H, "inf"))
     endif
     H = H(1:k1,1:k1);
     U = U(:,1:k1);
@@ -72,8 +72,8 @@
 
   ## tridiagonal H can still be rank deficient, so do permuted qr
   ## factorization
-  [qq,rr,pp] = qr(H);   # H = qq*rr*pp'
-  nn = rank(rr);
+  [qq, rr, pp] = qr (H);   # H = qq*rr*pp'
+  nn = rank (rr);
   qq = qq(:,1:nn);
   rr = rr(1:nn,:);            # rr may not be square, but "\" does least
   xx = U*pp*(rr\qq'*(U'*z));  # squares solution, so this works
@@ -86,38 +86,41 @@
   ## so for now I'm solving it with the krylov routine.
 
   ## initialize residual error norm
-  rnorm = norm(rk1,1);
+  rnorm = norm (rk1, 1);
 
   xnorm = 0;
-  fnorm = 1e-12 * norm([a,b;c,d],1);
+  fnorm = 1e-12 * norm ([a, b; c, d], 1);
 
-  ## dummy defines for MATHTOOLS compiler
-  gamk2 = 0;      omega1 = 0;      ztmz2 = 0;
+  gamk2 = 0;
+  omega1 = 0;
+  ztmz2 = 0;
 
   ## do until small changes to x
   len_x = length(x);
-  while ((k < 2*len_x) & (xnorm> 0.5) & (rnorm>fnorm))|(k == 0)
-    k = k+1;
+  while ((k < 2*len_x && xnorm > 0.5 && rnorm > fnorm) || k == 0)
+    k++;
 
     ## solve F_d z_{k-1} = r_{k-1}
-    zk1= zgfslv(n,m,p,rk1);
+    zk1= zgfslv (n, m, p, rk1);
 
     ## Generalized CG iteration
     ## gamk1 = (zk1'*F_d*zk1)/(zk1'*F*zk1);
     ztMz1 = zk1'*rk1;
-    gamk1 = ztMz1/(zk1'*zgfmul(a,b,c,d,zk1));
+    gamk1 = ztMz1/(zk1'*zgfmul (a, b, c, d, zk1));
 
-    if(rem(k,len_x) == 1) omega = 1;
-    else                  omega = 1/(1-gamk1*ztMz1/(gamk2*omega1*ztmz2));
+    if (rem (k, len_x) == 1)
+      omega = 1;
+    else
+      omega = 1/(1-gamk1*ztMz1/(gamk2*omega1*ztmz2));
     endif
 
     ## store x in xk2 to save space
     xk2 = xk2 + omega*(gamk1*zk1 + xk1 - xk2);
 
     ## compute new residual error: rk = z - F xk, check end conditions
-    rk1 = z - zgfmul(a,b,c,d,xk2);
-    rnorm = norm(rk1);
-    xnorm = max(abs(xk1 - xk2));
+    rk1 = z - zgfmul (a, b, c, d, xk2);
+    rnorm = norm (rk1);
+    xnorm = max (abs (xk1 - xk2));
 
     ## printf("zgscal: k=%d, gamk1=%e, gamk2=%e, \nztMz1=%e ztmz2=%e\n", ...
     ##   k,gamk1, gamk2, ztMz1, ztmz2);
@@ -129,21 +132,22 @@
     gamk2 = gamk1;
     omega1 = omega;
     ztmz2 = ztMz1;
-    [xk1,xk2] = swap(xk1,xk2);
+    [xk1, xk2] = swap (xk1, xk2);
   endwhile
   x = xk2;
 
   ## check convergence
-  if (xnorm> 0.5 & rnorm>fnorm)
-    warning("zgscal(tzero): GCG iteration failed; solving with pinv");
+  if (xnorm> 0.5 && rnorm > fnorm)
+    warning ("zgscal(tzero): GCG iteration failed; solving with pinv");
 
     ## perform brute force least squares; construct F
-    Am = eye(nmp);
-    for ii=1:nmp
-      Am(:,ii) = zgfmul(a,b,c,d,Am(:,ii));
+    Am = eye (nmp);
+    for ii = 1:nmp
+      Am(:,ii) = zgfmul (a, b, c, d, Am(:,ii));
     endfor
 
     ## now solve with qr factorization
-    x = pinv(Am)*z;
+    x = pinv (Am) * z;
   endif
+
 endfunction