--- a/extra/control-devel/devel/test_ab09jd.m	Thu Oct 20 16:56:06 2011 +0000
+++ b/extra/control-devel/devel/test_ab09jd.m	Thu Oct 20 19:50:05 2011 +0000
@@ -29,10 +29,34 @@
 
 dv = [       1 ];
 
+[ar, br, cr, dr] = slab09jd (a, b, c, d, 0, 0, 0, 1, 0.0, \
+                             1, av, bv, cv, dv, \
+                             0, [], [], [], [], \
+                             2, 1e-1, 1e-14)
+                             
+                             
+%{                             
+                             0, 0.0, \
+                             1, 0, 2, 0, 0, 1, \
+                             1e-1, 1e-14)
 
+  [ar, br, cr, dr, nr] = slab09jd (a, b, c, d, dt, scaled, nr, ordsel, alpha, \
+                                   jobv, av, bv, cv, dv, \
+                                   jobw, aw, bw, cw, dw, \
+                                   jobinv, tol1, tol2);
+%}
+%{
+sys = ss (a, b, c, d);
+sysv = ss (av, bv, cv, dv);
+
+sysr = hnamodred (sys, 0, sysv, [])
+%}
+
+%{
 [ar, br, cr, dr] = slab09jd (a, b, c, d, av, bv, cv, dv, [], [], [], [], 0, 0.0, \
                              1, 0, 2, 0, 0, 1, \
                              1e-1, 1e-14)
+%}
 
 %{
  The reduced state dynamics matrix Ar is 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/extra/control-devel/inst/hnamodred.m	Thu Oct 20 19:50:05 2011 +0000
@@ -0,0 +1,94 @@
+function [sysr, nr] = hnamodred (sys, nr = 0, sysv = [], sysw = [])
+
+  [a, b, c, d, tsam, scaled] = ssdata (sys);
+  dt = isdt (sys);
+  
+  if (isempty (sysv))
+    av = bv = cv = dv = [];
+    jobv = 0;
+  else
+    sysv = ss (sysv);
+    [av, bv, cv, dv] = ssdata (sysv);
+    jobv = 1;
+  endif
+
+  if (isempty (sysw))
+    aw = bw = cw = dw = [];
+    jobw = 0;
+  else
+    sysw = ss (sysw);
+    [aw, bw, cw, dw] = ssdata (sysw);
+    jobw = 1;
+  endif
+
+  jobinv = 2;
+  tol1 = 1e-1; 
+  tol2 = 1e-14;
+  alpha = 0.0;
+  ordsel = 1;
+
+  [ar, br, cr, dr, nr] = slab09jd (a, b, c, d, dt, scaled, nr, ordsel, alpha, \
+                                   jobv, av, bv, cv, dv, \
+                                   jobw, aw, bw, cw, dw, \
+                                   jobinv, tol1, tol2);
+
+  sysr = ss (ar, br, cr, dr, tsam);
+
+endfunction
+
+
+%!shared Mo, Me
+%! a =  [ -3.8637   -7.4641   -9.1416   -7.4641   -3.8637   -1.0000
+%!         1.0000,         0         0         0         0         0
+%!              0    1.0000         0         0         0         0
+%!              0         0    1.0000         0         0         0
+%!              0         0         0    1.0000         0         0
+%!              0         0         0         0    1.0000         0 ];
+%!
+%! b =  [       1
+%!              0
+%!              0
+%!              0
+%!              0
+%!              0 ];
+%!
+%! c =  [       0         0         0         0         0         1 ];
+%!
+%! d =  [       0 ];
+%!
+%! sys = ss (a, b, c, d);
+%!
+%! av = [  0.2000   -1.0000
+%!         1.0000         0 ];
+%!
+%! bv = [       1
+%!              0 ];
+%!
+%! cv = [ -1.8000         0 ];
+%!
+%! dv = [       1 ];
+%!
+%! sysv = ss (av, bv, cv, dv);
+%!
+%! sysr = hnamodred (sys, 0, sysv, []);
+%! [ao, bo, co, do] = ssdata (sysr);
+%!
+%! ae = [ -0.2391   0.3072   1.1630   1.1967
+%!        -2.9709  -0.2391   2.6270   3.1027
+%!         0.0000   0.0000  -0.5137  -1.2842
+%!         0.0000   0.0000   0.1519  -0.5137 ];
+%!
+%! be = [ -1.0497
+%!        -3.7052
+%!         0.8223
+%!         0.7435 ];
+%!
+%! ce = [ -0.4466   0.0143  -0.4780  -0.2013 ];
+%!
+%! de = [  0.0219 ];
+%!
+%! Mo = [ao, bo; co, do];
+%! Me = [ae, be; ce, de];
+%!
+%!assert (Mo, Me, 1e-4);
+

--- a/extra/control-devel/src/slab09jd.cc	Thu Oct 20 16:56:06 2011 +0000
+++ b/extra/control-devel/src/slab09jd.cc	Thu Oct 20 19:50:05 2011 +0000
@@ -87,26 +87,25 @@
         Matrix c = args(2).matrix_value ();
         Matrix d = args(3).matrix_value ();
         
-        Matrix av = args(4).matrix_value ();
-        Matrix bv = args(5).matrix_value ();
-        Matrix cv = args(6).matrix_value ();
-        Matrix dv = args(7).matrix_value ();
-        
-        Matrix aw = args(8).matrix_value ();
-        Matrix bw = args(9).matrix_value ();
-        Matrix cw = args(10).matrix_value ();
-        Matrix dw = args(11).matrix_value ();
+        const int idico = args(4).int_value ();
+        const int iequil = args(5).int_value ();
+        int nr = args(6).int_value ();
+        const int iordsel = args(7).int_value ();
+        double alpha = args(8).double_value ();
+                 
+        const int ijobv = args(9).int_value ();       
+        Matrix av = args(10).matrix_value ();
+        Matrix bv = args(11).matrix_value ();
+        Matrix cv = args(12).matrix_value ();
+        Matrix dv = args(13).matrix_value ();
 
-        int nr = args(12).int_value ();
-        double alpha = args(13).double_value ();
-        
-        const int ijobv = args(14).int_value ();
-        const int ijobw = args(15).int_value ();
-        const int ijobinv = args(16).int_value ();
-        const int idico = args(17).int_value ();
-        const int iequil = args(18).int_value ();
-        const int iordsel = args(19).int_value ();
-        
+        const int ijobw = args(14).int_value ();        
+        Matrix aw = args(15).matrix_value ();
+        Matrix bw = args(16).matrix_value ();
+        Matrix cw = args(17).matrix_value ();
+        Matrix dw = args(18).matrix_value ();
+
+        const int ijobinv = args(19).int_value ();
         double tol1 = args(20).double_value ();
         double tol2 = args(21).double_value ();
 
@@ -308,6 +307,7 @@
         retval(1) = b;
         retval(2) = c;
         retval(3) = d;
+        retval(4) = octave_value (nr);
         // retval(0) = hsv;
         // retval(1) = octave_value (ns);
     }