--- a/extra/control-oo/inst/care.m	Fri Feb 12 22:52:54 2010 +0000
+++ b/extra/control-oo/inst/care.m	Sat Feb 13 13:06:39 2010 +0000
@@ -1,4 +1,4 @@
-## Copyright (C) 2009   Lukas F. Reichlin
+## Copyright (C) 2009 - 2010   Lukas F. Reichlin
 ##
 ## This file is part of LTI Syncope.
 ##
@@ -21,56 +21,55 @@
 ## Return unique stabilizing solution x of the continuous-time
 ## riccati equation as well as the closed-loop poles l and the
 ## corresponding gain matrix g.
+## Uses SLICOT SB02OD by courtesy of NICONET e.V.
+## <http://www.slicot.org>
 ## @end deftypefn
 
 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
 ## Created: November 2009
-## Version: 0.1
+## Version: 0.2
 
-function [x, l, g] = care (a, b, q, r, s = [], opt = "B")
+function [x, l, g] = care (a, b, q, r, s = [])
 
-  if (nargin < 4 || nargin > 6)
+  if (nargin < 4 || nargin > 5)
     print_usage ();
   endif
 
-  if (nargin == 6)
-    if (ischar (opt))
-      opt = upper (opt(1));
-      if (opt != "B" && opt != N && opt != "P" && opt != "S")
-        warning ("dare: opt has invalid value ""%s""; setting to ""B""", opt);
-        opt = "B";
-      endif
-    else
-      warning ("dare: invalid argument opt, setting to ""B""");
-      opt = "B";
-    endif
+  asize = issquare (a);
+  [brows, bcols] = size (b);
+  qsize = issquare (q);
+  rsize = issquare (r);
+
+  if (! asize)
+    error ("care: a is not square");
   endif
 
-  [n, m] = size (b);
-  p = issquare (q);
-  m1 = issquare (r);
-
-  if (! m1)
-    error ("care: r is not square");
-  elseif (m1 != m)
-    error ("care: b, r are not conformable");
+  if (! qsize)
+    error ("care: q is not square");
   endif
 
-  if (! p)
-    q = q' * q;
+  if (! rsize)
+    error ("care: r is not square");
+  endif
+  
+  if (asize != brows)
+    error ("care: a, b are not conformable");
+  endif
+  
+  if (rsize != bcols)
+    error ("care: b, r are not conformable");
   endif
 
   ## incorporate cross term into a and q
   if (isempty (s))
-    s = zeros (n, m);
     ao = a;
     qo = q;
   else
-    [n2, m2] = size (s);
+    [srows, scols] = size (s);  % [n2, m2]
 
-    if (n2 != n || m2 != m)
+    if (srows != brows || scols != bcols)
       error ("care: s (%dx%d) must be identically dimensioned with b (%dx%d)",
-              n2, m2, n, m);
+              srows, scols, brows, bcols);
     endif
 
     ao = a - (b/r)*s';
@@ -86,26 +85,38 @@
   dflag = isdetectable (ao, qo, [], 0);
 
   if (dflag == 0)
-    warning ("care: a and q not detectable");
+    warning ("care: (a,q) not detectable");
   elseif (dflag == -1)
-    error ("care: a and q have poles on imaginary axis");
+    error ("care: (a,q) has poles on imaginary axis");
   endif
 
   ## to allow lqe design, don't force
   ## qo to be positive semi-definite
 
+  ## checking positive definiteness
   if (isdefinite (r) <= 0)
     error ("care: r must be positive definite");
   endif
 
-  ## unique stabilizing solution
-  x = are (ao, (b/r)*b', qo, opt);
-
-  ## corresponding gain matrix
-  g = r \ (b'*x + s');
+  ## solve the riccati equation
+  if (isempty (s))
+    ## unique stabilizing solution
+    x = slsb02od (a, b, q, r, b, false, false);
+    
+    ## corresponding gain matrix
+    g = r \ (b'*x);
+  else
+    ## unique stabilizing solution
+    x = slsb02od (a, b, q, r, s, false, true);
+    
+    ## corresponding gain matrix
+    g = r \ (b'*x + s');
+  endif
 
   ## closed-loop poles
   l = eig (a - b*g);
+  
+  ## TODO: use alphar, alphai and beta from SB02OD
 
 endfunction
 
@@ -133,4 +144,33 @@
 %!
 %!assert (x, xe, 1e-4);
 %!assert (l, le, 1e-4);
+%!assert (g, ge, 1e-4);
+
+%!shared x, l, g, xe, le, ge
+%! a = [ 0.0  1.0
+%!       0.0  0.0];
+%!
+%! b = [ 0.0
+%!       1.0];
+%!
+%! c = [ 1.0  0.0
+%!       0.0  1.0
+%!       0.0  0.0];
+%!
+%! d = [ 0.0
+%!       0.0
+%!       1.0];
+%!
+%! [x, l, g] = care (a, b, c'*c, d'*d);
+%!
+%! xe = [ 1.7321   1.0000
+%!        1.0000   1.7321];
+%!
+%! le = [-0.8660 + 0.5000i
+%!       -0.8660 - 0.5000i];
+%!
+%! ge = [ 1.0000   1.7321];
+%!
+%!assert (x, xe, 1e-4);
+%!assert (l, le, 1e-4);
 %!assert (g, ge, 1e-4);
\ No newline at end of file

--- a/extra/control-oo/inst/dare.m	Fri Feb 12 22:52:54 2010 +0000
+++ b/extra/control-oo/inst/dare.m	Sat Feb 13 13:06:39 2010 +0000
@@ -1,135 +1,81 @@
-## Copyright (C) 1996, 1997, 2000, 2002, 2003, 2004, 2005, 2006, 2007
-##               Auburn University.  All rights reserved.
+## Copyright (C) 2009 - 2010   Lukas F. Reichlin
 ##
+## This file is part of LTI Syncope.
 ##
-## This program 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 3 of the License, or (at
-## your option) any later version.
+## LTI Syncope 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 3 of the License, or
+## (at your option) any later version.
 ##
-## This program 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.
+## LTI Syncope 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 this program; see the file COPYING.  If not, see
-## <http://www.gnu.org/licenses/>.
+## along with this program. If not, see <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {[@var{x}, @var{l}, @var{g}] =} dare (@var{a}, @var{b}, @var{q}, @var{r})
 ## @deftypefnx {Function File} {[@var{x}, @var{l}, @var{g}] =} dare (@var{a}, @var{b}, @var{q}, @var{r}, @var{s})
-## @deftypefnx {Function File} {[@var{x}, @var{l}, @var{g}] =} dare (@var{a}, @var{b}, @var{q}, @var{r}, @var{s}, @var{opt})
-## Return the solution, @var{x} of the discrete-time algebraic Riccati
-## equation
-## @iftex
-## @tex
-## $$
-## A^TXA - X + A^TXB  (R + B^TXB)^{-1} B^TXA + Q = 0
-## $$
-## @end tex
-## @end iftex
-## @ifinfo
-## @example
-## a' x a - x + a' x b (r + b' x b)^(-1) b' x a + q = 0
-## @end example
-## @end ifinfo
-## @noindent
-##
-## @strong{Inputs}
-## @table @var
-## @item a
-## @var{n} by @var{n} matrix;
-##
-## @item b
-## @var{n} by @var{m} matrix;
-##
-## @item q
-## @var{n} by @var{n} matrix, symmetric positive semidefinite, or a @var{p} by @var{n} matrix,
-## In the latter case @math{q:=q'*q} is used;
-##
-## @item r
-## @var{m} by @var{m}, symmetric positive definite (invertible);
-##
-## @item opt
-## (optional argument; default = @code{"B"}):
-## String option passed to @code{balance} prior to ordered @var{QZ} decomposition.
-## @end table
-##
-## @strong{Output}
-## @table @var
-## @item x
-## solution of @acronym{DARE}.
-## @end table
-##
-## @strong{Method}
-## Generalized eigenvalue approach (Van Dooren; @acronym{SIAM} J.
-##  Sci. Stat. Comput., Vol 2) applied  to the appropriate symplectic pencil.
-##
-##  See also: Ran and Rodman, @cite{Stable Hermitian Solutions of Discrete
-##  Algebraic Riccati Equations}, Mathematics of Control, Signals and
-##  Systems, Vol 5, no 2 (1992), pp 165--194.
-## @seealso{balance, are}
+## Return unique stabilizing solution x of the discrete-time
+## riccati equation as well as the closed-loop poles l and the
+## corresponding gain matrix g.
+## Uses SLICOT SB02OD by courtesy of NICONET e.V.
+## <http://www.slicot.org>
 ## @end deftypefn
 
-## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu>
-## Created: August 1993
-## Adapted-By: jwe
+## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
+## Created: October 2009
+## Version: 0.2
 
-## Adapted-By: Lukas Reichlin <lukas.reichlin@gmail.com>
-## Date: October 2009
-## Version: 0.1
+function [x, l, g] = dare (a, b, q, r, s = [])
 
-function [x, l, g] = dare (a, b, q, r, s = [], opt = "B")
-
-  if (nargin < 4 || nargin > 6)
+  if (nargin < 4 || nargin > 5)
     print_usage ();
   endif
 
-  if (nargin == 6)
-    if (ischar (opt))
-      opt = upper (opt(1));
-      if (opt != "B" && opt != N && opt != "P" && opt != "S")
-        warning ("dare: opt has invalid value ""%s""; setting to ""B""", opt);
-        opt = "B";
-      endif
-    else
-      warning ("dare: invalid argument opt, setting to ""B""");
-      opt = "B";
-    endif
+  asize = issquare (a);
+  [brows, bcols] = size (b);
+  qsize = issquare (q);
+  rsize = issquare (r);
+
+  if (! asize)
+    error ("dare: a is not square");
   endif
 
-  [n, m] = size (b);
-  p = issquare (q);
-  m1 = issquare (r);
+  if (! qsize)
+    error ("dare: q is not square");
+  endif
 
-  if (! m1)
+  if (! rsize)
     error ("dare: r is not square");
-  elseif (m1 != m)
+  endif
+  
+  if (asize != brows)
+    error ("dare: a, b are not conformable");
+  endif
+  
+  if (rsize != bcols)
     error ("dare: b, r are not conformable");
   endif
 
-  if (! p)
-    q = q' * q;
-  endif
-
-  ## incorporate cross term into a and q.
+  ## incorporate cross term into a and q
   if (isempty (s))
-    s = zeros (n, m);
     ao = a;
     qo = q;
   else
-    [n2, m2] = size (s);
+    [srows, scols] = size (s);  % [n2, m2]
 
-    if (n2 != n || m2 != m)
+    if (srows != brows || scols != bcols)
       error ("dare: s (%dx%d) must be identically dimensioned with b (%dx%d)",
-              n2, m2, n, m);
+              srows, scols, brows, bcols);
     endif
 
     ao = a - (b/r)*s';
     qo = q - (s/r)*s';
   endif
-
+  
   ## check stabilizability
   if (! isstabilizable (ao, b, [], 1))
     error ("dare: a and b not stabilizable");
@@ -144,35 +90,35 @@
     error ("dare: (a,q) has non-minimal modes near unit circle");
   endif
 
-  ## Checking positive definiteness
+  ## to allow lqe design, don't force
+  ## qo to be positive semi-definite
+
+  ## checking positive definiteness
   if (isdefinite (r) <= 0)
     error ("dare: r must be positive definite");
   endif
 
-  ## if (isdefinite (qo) < 0)
-  ##   error ("dare: q not positive semidefinite");
-  ## endif
-
   ## solve the riccati equation
-  s1 = [ ao, zeros(n);
-        -qo,  eye(n) ];
-  s2 = [ eye(n),  (b/r)*b';
-        zeros(n),    ao'  ];
-
-  [c, d, s1, s2] = balance (s1, s2, opt);
-  [aa, bb, u, lam] = qz (s1, s2, "S");
-
-  u = d*u;
-  n1 = n+1;
-  n2 = 2*n;
-
-  ## unique stabilizing solution
-  x = u(n1:n2, 1:n) / u(1:n, 1:n);
-
-  ## corresponding gain matrix
-  g = (r+b'*x*b) \ (b'*x*a + s');
+  if (isempty (s))
+    ## unique stabilizing solution
+    x = slsb02od (a, b, q, r, b, true, false);
+    
+    ## corresponding gain matrix
+    g = (r+b'*x*b) \ (b'*x*a);
+  else
+    ## unique stabilizing solution
+    x = slsb02od (a, b, q, r, s, true, true);
+    
+    ## corresponding gain matrix
+    g = (r+b'*x*b) \ (b'*x*a + s');
+  endif
 
   ## closed-loop poles
   l = eig (a - b*g);
+  
+  ## TODO: use alphar, alphai and beta from SB02OD
 
 endfunction
+
+
+## TODO: add a test
\ No newline at end of file

--- a/extra/control-oo/src/Makefile	Fri Feb 12 22:52:54 2010 +0000
+++ b/extra/control-oo/src/Makefile	Sat Feb 13 13:06:39 2010 +0000
@@ -1,6 +1,6 @@
 all: slab08nd.oct slab13dd.oct slsb10hd.oct slsb10ed.oct slab13bd.oct \
      slsb01bd.oct slsb10fd.oct slsb10dd.oct slsb03md.oct slsb04md.oct \
-     slsb04qd.oct slsg03ad.oct
+     slsb04qd.oct slsg03ad.oct slsb02od.oct
 
 # transmission zeros of state-space models
 slab08nd.oct: slab08nd.cc
@@ -91,5 +91,11 @@
               SG03AD.f MB01RW.f MB01RD.f SG03AX.f SG03AY.f \
               MB02UU.f MB02UV.f
 
+# algebraic Riccati equations
+slsb02od.oct: slsb02od.cc
+	mkoctfile slsb02od.cc \
+              SB02OD.f SB02OU.f SB02OV.f SB02OW.f SB02OY.f \
+              SB02MR.f SB02MV.f
+
 clean:
 	rm *.o core octave-core *.oct *~

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/extra/control-oo/src/slsb02od.cc	Sat Feb 13 13:06:39 2010 +0000
@@ -0,0 +1,203 @@
+/*
+
+Copyright (C) 2010   Lukas F. Reichlin
+
+This file is part of LTI Syncope.
+
+LTI Syncope 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 3 of the License, or
+(at your option) any later version.
+
+LTI Syncope 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 this program. If not, see <http://www.gnu.org/licenses/>.
+
+Solve algebraic Riccati equation.
+Uses SLICOT SB02OD by courtesy of NICONET e.V.
+<http://www.slicot.org>
+
+Author: Lukas Reichlin <lukas.reichlin@gmail.com>
+Created: February 2010
+Version: 0.1
+
+*/
+
+#include <octave/oct.h>
+#include <f77-fcn.h>
+
+extern "C"
+{ 
+    int F77_FUNC (sb02od, SB02OD)
+                 (char& DICO, char& JOBB,
+                  char& FACT, char& UPLO,
+                  char& JOBL, char& SORT,
+                  int& N, int& M, int& P,
+                  double* A, int& LDA,
+                  double* B, int& LDB,
+                  double* Q, int& LDQ,
+                  double* R, int& LDR,
+                  double* L, int& LDL,
+                  double& RCOND,
+                  double* X, int& LDX,
+                  double* ALFAR, double* ALFAI,
+                  double* BETA,
+                  double* S, int& LDS,
+                  double* T, int& LDT,
+                  double* U, int& LDU,
+                  double& TOL,
+                  int* IWORK,
+                  double* DWORK, int& LDWORK,
+                  bool* BWORK,
+                  int& INFO);
+}
+
+int max (int a, int b)
+{
+    if (a > b)
+        return a;
+    else
+        return b;
+}
+
+int max (int a, int b, int c)
+{
+    int d = max (a, b);
+    
+    return max (c, d);
+}
+
+int max (int a, int b, int c, int d)
+{
+    int e = max (a, b);
+    int f = max (c, d);
+    
+    return max (e, f);
+}
+     
+DEFUN_DLD (slsb02od, args, nargout, "Slicot SB02OD Release 5.0")
+{
+    int nargin = args.length ();
+    octave_value_list retval;
+    
+    if (nargin != 7)
+    {
+        print_usage ();
+    }
+    else
+    {
+        // arguments in
+        char dico;
+        char jobb = 'B';
+        char fact = 'N';
+        char uplo = 'U';
+        char jobl;
+        char sort = 'S';
+        
+        NDArray a = args(0).array_value ();
+        NDArray b = args(1).array_value ();
+        NDArray q = args(2).array_value ();
+        NDArray r = args(3).array_value ();
+        NDArray l = args(4).array_value ();
+        int digital = args(5).int_value ();
+        int ijobl = args(6).int_value ();
+
+        if (digital == 0)
+            dico = 'C';
+        else
+            dico = 'D';
+            
+        if (ijobl == 0)
+            jobl = 'Z';
+        else
+            jobl = 'N';
+        
+        int n = a.rows ();      // n: number of states
+        int m = b.columns ();   // m: number of inputs
+        int p = 0;              // p: number of outputs, not used because FACT = 'N'
+        
+        int lda = max (1, n);
+        int ldb = max (1, n);
+        int ldq = max (1, n);
+        int ldr = max (1, m);
+        int ldl = max (1, n);
+        
+        // arguments out
+        double rcond;
+        
+        int ldx = max (1, n);
+        dim_vector dv (2);
+        dv(0) = ldx;
+        dv(1) = n;
+        NDArray x (dv);
+        
+        // unused output arguments
+        OCTAVE_LOCAL_BUFFER (double, alfar, 2*n);
+        OCTAVE_LOCAL_BUFFER (double, alfai, 2*n);
+        OCTAVE_LOCAL_BUFFER (double, beta, 2*n);
+        
+        int lds = max (1, 2*n + m);
+        OCTAVE_LOCAL_BUFFER (double, s, lds*lds);
+        
+        int ldt = max (1, 2*n + m);
+        OCTAVE_LOCAL_BUFFER (double, t, ldt * 2*n);
+        
+        int ldu = max (1, 2*n);
+        OCTAVE_LOCAL_BUFFER (double, u, ldu * 2*n);
+        
+        // tolerance
+        double tol = 0;  // use default value
+        
+        // workspace
+        int liwork = max (1, m, 2*n);
+        OCTAVE_LOCAL_BUFFER (int, iwork, liwork);
+        
+        int ldwork = max (7*(2*n + 1) + 16, 16*n, 2*n + m, 3*m);
+        OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
+        
+        OCTAVE_LOCAL_BUFFER (bool, bwork, 2*n);
+        
+        // error indicator
+        int info;
+
+
+        // SLICOT routine SB02OD
+        F77_XFCN (sb02od, SB02OD,
+                 (dico, jobb,
+                  fact, uplo,
+                  jobl, sort,
+                  n, m, p,
+                  a.fortran_vec (), lda,
+                  b.fortran_vec (), ldb,
+                  q.fortran_vec (), ldq,
+                  r.fortran_vec (), ldr,
+                  l.fortran_vec (), ldl,
+                  rcond,
+                  x.fortran_vec (), ldx,
+                  alfar, alfai,
+                  beta,
+                  s, lds,
+                  t, ldt,
+                  u, ldu,
+                  tol,
+                  iwork,
+                  dwork, ldwork,
+                  bwork,
+                  info));
+
+        if (f77_exception_encountered)
+            error ("lti: norm: slsb02od: exception in SLICOT subroutine SB02OD");
+            
+        if (info != 0)
+            error ("lti: norm: slsb02od: SB02OD returned info = %d", info);
+        
+        // return value
+        retval(0) = x;
+    }
+    
+    return retval;
+}