--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/corefcn/schur.cc	Fri Jul 27 15:35:00 2012 -0700
@@ -0,0 +1,381 @@
+/*
+
+Copyright (C) 1996-2012 John W. Eaton
+
+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 3 of the License, 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, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <string>
+
+#include "CmplxSCHUR.h"
+#include "dbleSCHUR.h"
+#include "fCmplxSCHUR.h"
+#include "floatSCHUR.h"
+
+#include "defun.h"
+#include "error.h"
+#include "gripes.h"
+#include "oct-obj.h"
+#include "utils.h"
+
+template <class Matrix>
+static octave_value
+mark_upper_triangular (const Matrix& a)
+{
+  octave_value retval = a;
+
+  octave_idx_type n = a.rows ();
+  assert (a.columns () == n);
+
+  const typename Matrix::element_type zero = typename Matrix::element_type ();
+
+  for (octave_idx_type i = 0; i < n; i++)
+    if (a(i,i) == zero)
+      return retval;
+
+  retval.matrix_type (MatrixType::Upper);
+
+  return retval;
+}
+
+DEFUN (schur, args, nargout,
+  "-*- texinfo -*-\n\
+@deftypefn  {Built-in Function} {@var{S} =} schur (@var{A})\n\
+@deftypefnx {Built-in Function} {@var{S} =} schur (@var{A}, \"real\")\n\
+@deftypefnx {Built-in Function} {@var{S} =} schur (@var{A}, \"complex\")\n\
+@deftypefnx {Built-in Function} {@var{S} =} schur (@var{A}, @var{opt})\n\
+@deftypefnx {Built-in Function} {[@var{U}, @var{S}] =} schur (@var{A}, @dots{})\n\
+@cindex Schur decomposition\n\
+Compute the Schur@tie{}decomposition of @var{A}\n\
+@tex\n\
+$$\n\
+ S = U^T A U\n\
+$$\n\
+@end tex\n\
+@ifnottex\n\
+\n\
+@example\n\
+@code{@var{S} = @var{U}' * @var{A} * @var{U}}\n\
+@end example\n\
+\n\
+@end ifnottex\n\
+where @var{U} is a unitary matrix\n\
+@tex\n\
+($U^T U$ is identity)\n\
+@end tex\n\
+@ifnottex\n\
+(@code{@var{U}'* @var{U}} is identity)\n\
+@end ifnottex\n\
+and @var{S} is upper triangular.  The eigenvalues of @var{A} (and @var{S})\n\
+are the diagonal elements of @var{S}.  If the matrix @var{A}\n\
+is real, then the real Schur@tie{}decomposition is computed, in which the\n\
+matrix @var{U} is orthogonal and @var{S} is block upper triangular\n\
+with blocks of size at most\n\
+@tex\n\
+$2 \\times 2$\n\
+@end tex\n\
+@ifnottex\n\
+@code{2 x 2}\n\
+@end ifnottex\n\
+along the diagonal.  The diagonal elements of @var{S}\n\
+(or the eigenvalues of the\n\
+@tex\n\
+$2 \\times 2$\n\
+@end tex\n\
+@ifnottex\n\
+@code{2 x 2}\n\
+@end ifnottex\n\
+blocks, when appropriate) are the eigenvalues of @var{A} and @var{S}.\n\
+\n\
+The default for real matrices is a real Schur@tie{}decomposition.\n\
+A complex decomposition may be forced by passing the flag \"complex\".\n\
+\n\
+The eigenvalues are optionally ordered along the diagonal according to\n\
+the value of @var{opt}.  @code{@var{opt} = \"a\"} indicates that all\n\
+eigenvalues with negative real parts should be moved to the leading\n\
+block of @var{S}\n\
+(used in @code{are}), @code{@var{opt} = \"d\"} indicates that all eigenvalues\n\
+with magnitude less than one should be moved to the leading block of @var{S}\n\
+(used in @code{dare}), and @code{@var{opt} = \"u\"}, the default, indicates\n\
+that no ordering of eigenvalues should occur.  The leading @var{k}\n\
+columns of @var{U} always span the @var{A}-invariant\n\
+subspace corresponding to the @var{k} leading eigenvalues of @var{S}.\n\
+\n\
+The Schur@tie{}decomposition is used to compute eigenvalues of a\n\
+square matrix, and has applications in the solution of algebraic\n\
+Riccati equations in control (see @code{are} and @code{dare}).\n\
+@seealso{rsf2csf}\n\
+@end deftypefn")
+{
+  octave_value_list retval;
+
+  int nargin = args.length ();
+
+  if (nargin < 1 || nargin > 2 || nargout > 2)
+    {
+      print_usage ();
+      return retval;
+    }
+
+  octave_value arg = args(0);
+
+  std::string ord;
+
+  if (nargin == 2)
+    {
+      ord = args(1).string_value ();
+
+      if (error_state)
+        {
+          error ("schur: second argument must be a string");
+          return retval;
+        }
+    }
+
+  bool force_complex = false;
+
+  if (ord == "real")
+    {
+      ord = std::string ();
+    }
+  else if (ord == "complex")
+    {
+      force_complex = true;
+      ord = std::string ();
+    }
+  else
+    {
+      char ord_char = ord.empty () ? 'U' : ord[0];
+
+      if (ord_char != 'U' && ord_char != 'A' && ord_char != 'D'
+          && ord_char != 'u' && ord_char != 'a' && ord_char != 'd')
+        {
+          warning ("schur: incorrect ordered schur argument `%c'",
+                   ord.c_str ());
+          return retval;
+        }
+    }
+
+  octave_idx_type nr = arg.rows ();
+  octave_idx_type nc = arg.columns ();
+
+  if (nr != nc)
+    {
+      gripe_square_matrix_required ("schur");
+      return retval;
+    }
+
+  if (! arg.is_numeric_type ())
+    gripe_wrong_type_arg ("schur", arg);
+  else if (arg.is_single_type ())
+    {
+      if (! force_complex && arg.is_real_type ())
+        {
+          FloatMatrix tmp = arg.float_matrix_value ();
+
+          if (! error_state)
+            {
+              if (nargout == 0 || nargout == 1)
+                {
+                  FloatSCHUR result (tmp, ord, false);
+                  retval(0) = result.schur_matrix ();
+                }
+              else
+                {
+                  FloatSCHUR result (tmp, ord, true);
+                  retval(1) = result.schur_matrix ();
+                  retval(0) = result.unitary_matrix ();
+                }
+            }
+        }
+      else
+        {
+          FloatComplexMatrix ctmp = arg.float_complex_matrix_value ();
+
+          if (! error_state)
+            {
+
+              if (nargout == 0 || nargout == 1)
+                {
+                  FloatComplexSCHUR result (ctmp, ord, false);
+                  retval(0) = mark_upper_triangular (result.schur_matrix ());
+                }
+              else
+                {
+                  FloatComplexSCHUR result (ctmp, ord, true);
+                  retval(1) = mark_upper_triangular (result.schur_matrix ());
+                  retval(0) = result.unitary_matrix ();
+                }
+            }
+        }
+    }
+  else
+    {
+      if (! force_complex && arg.is_real_type ())
+        {
+          Matrix tmp = arg.matrix_value ();
+
+          if (! error_state)
+            {
+              if (nargout == 0 || nargout == 1)
+                {
+                  SCHUR result (tmp, ord, false);
+                  retval(0) = result.schur_matrix ();
+                }
+              else
+                {
+                  SCHUR result (tmp, ord, true);
+                  retval(1) = result.schur_matrix ();
+                  retval(0) = result.unitary_matrix ();
+                }
+            }
+        }
+      else
+        {
+          ComplexMatrix ctmp = arg.complex_matrix_value ();
+
+          if (! error_state)
+            {
+
+              if (nargout == 0 || nargout == 1)
+                {
+                  ComplexSCHUR result (ctmp, ord, false);
+                  retval(0) = mark_upper_triangular (result.schur_matrix ());
+                }
+              else
+                {
+                  ComplexSCHUR result (ctmp, ord, true);
+                  retval(1) = mark_upper_triangular (result.schur_matrix ());
+                  retval(0) = result.unitary_matrix ();
+                }
+            }
+        }
+    }
+
+  return retval;
+}
+
+/*
+%!test
+%! a = [1, 2, 3; 4, 5, 9; 7, 8, 6];
+%! [u, s] = schur (a);
+%! assert (u' * a * u, s, sqrt (eps));
+
+%!test
+%! a = single ([1, 2, 3; 4, 5, 9; 7, 8, 6]);
+%! [u, s] = schur (a);
+%! assert (u' * a * u, s, sqrt (eps ("single")));
+
+%!test
+%! fail ("schur ([1, 2; 3, 4], 2)", "warning");
+
+%!error schur ()
+%!error <argument must be a square matrix> schur ([1, 2, 3; 4, 5, 6])
+*/
+
+DEFUN (rsf2csf, args, nargout,
+  "-*- texinfo -*-\n\
+@deftypefn {Function File} {[@var{U}, @var{T}] =} rsf2csf (@var{UR}, @var{TR})\n\
+Convert a real, upper quasi-triangular Schur@tie{}form @var{TR} to a complex,\n\
+upper triangular Schur@tie{}form @var{T}.\n\
+\n\
+Note that the following relations hold:\n\
+\n\
+@tex\n\
+$UR \\cdot TR \\cdot {UR}^T = U T U^{\\dagger}$ and\n\
+$U^{\\dagger} U$ is the identity matrix I.\n\
+@end tex\n\
+@ifnottex\n\
+@xcode{@var{UR} * @var{TR} * @var{UR}' = @var{U} * @var{T} * @var{U}'} and\n\
+@code{@var{U}' * @var{U}} is the identity matrix I.\n\
+@end ifnottex\n\
+\n\
+Note also that @var{U} and @var{T} are not unique.\n\
+@seealso{schur}\n\
+@end deftypefn")
+{
+  octave_value_list retval;
+
+  if (args.length () == 2 && nargout <= 2)
+    {
+      if (! args(0).is_numeric_type ())
+        gripe_wrong_type_arg ("rsf2csf", args(0));
+      else if (! args(1).is_numeric_type ())
+        gripe_wrong_type_arg ("rsf2csf", args(1));
+      else if (args(0).is_complex_type () || args(1).is_complex_type ())
+        error ("rsf2csf: UR and TR must be real matrices");
+      else
+        {
+
+          if (args(0).is_single_type () || args(1).is_single_type ())
+            {
+              FloatMatrix u = args(0).float_matrix_value ();
+              FloatMatrix t = args(1).float_matrix_value ();
+              if (! error_state)
+                {
+                  FloatComplexSCHUR cs (FloatSCHUR (t, u));
+
+                  retval(1) = cs.schur_matrix ();
+                  retval(0) = cs.unitary_matrix ();
+                }
+            }
+          else
+            {
+              Matrix u = args(0).matrix_value ();
+              Matrix t = args(1).matrix_value ();
+              if (! error_state)
+                {
+                  ComplexSCHUR cs (SCHUR (t, u));
+
+                  retval(1) = cs.schur_matrix ();
+                  retval(0) = cs.unitary_matrix ();
+                }
+            }
+        }
+    }
+  else
+    print_usage ();
+
+  return retval;
+}
+
+/*
+%!test
+%! A = [1, 1, 1, 2; 1, 2, 1, 1; 1, 1, 3, 1; -2, 1, 1, 1];
+%! [u, t] = schur (A);
+%! [U, T] = rsf2csf (u, t);
+%! assert (norm (u * t * u' - U * T * U'), 0, 1e-12);
+%! assert (norm (A - U * T * U'), 0, 1e-12);
+
+%!test
+%! A = rand (10);
+%! [u, t] = schur (A);
+%! [U, T] = rsf2csf (u, t);
+%! assert (norm (tril (T, -1)), 0);
+%! assert (norm (U * U'), 1, 1e-14);
+
+%!test
+%! A = [0, 1;-1, 0];
+%! [u, t] = schur (A);
+%! [U, T] = rsf2csf (u,t);
+%! assert (U * T * U', A, 1e-14);
+*/