Mercurial > octave-nkf
changeset 19781:56157a7505ed
Add new ordschur function.
* libinterp/corefcn/ordschur.cc: New file.
* libinterp/corefcn/module.mk: Include it in the list of source files.
* scripts/help/__unimplemented__.m: Remove ordschur from the list of unimplemented functions.
* doc/interpreter/linalg.txi: Add it to the interpreter manual.
* NEWS: Mention it.
* libinterp/corefcn/schur.cc (schur): Reference it from the documentation of the schur function.
Thanks to Carnë Draug for improving the original patch, and to Mike Miller for reviewing it and suggesting improvements.
| author | Sébastien Villemot <sebastien@debian.org> |
|---|---|
| date | Sat, 07 Feb 2015 21:51:20 +0100 |
| parents | 97690ea6f57a |
| children | 3fc946d5e91f |
| files | NEWS doc/interpreter/linalg.txi libinterp/corefcn/module.mk libinterp/corefcn/ordschur.cc libinterp/corefcn/schur.cc scripts/help/__unimplemented__.m |
| diffstat | 6 files changed, 262 insertions(+), 2 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS Wed Feb 18 09:20:59 2015 -0800 +++ b/NEWS Sat Feb 07 21:51:20 2015 +0100 @@ -175,6 +175,7 @@ linkaxes lscov numfields + ordschur qmr rotate sylvester
--- a/doc/interpreter/linalg.txi Wed Feb 18 09:20:59 2015 -0800 +++ b/doc/interpreter/linalg.txi Sat Feb 07 21:51:20 2015 +0100 @@ -175,6 +175,8 @@ @DOCSTRING(rsf2csf) +@DOCSTRING(ordschur) + @DOCSTRING(subspace) @DOCSTRING(svd)
--- a/libinterp/corefcn/module.mk Wed Feb 18 09:20:59 2015 -0800 +++ b/libinterp/corefcn/module.mk Sat Feb 07 21:51:20 2015 +0100 @@ -220,6 +220,7 @@ corefcn/oct-stream.cc \ corefcn/oct-strstrm.cc \ corefcn/octave-link.cc \ + corefcn/ordschur.cc \ corefcn/pager.cc \ corefcn/pinv.cc \ corefcn/pr-output.cc \
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libinterp/corefcn/ordschur.cc Sat Feb 07 21:51:20 2015 +0100 @@ -0,0 +1,257 @@ +/* + +Copyright (C) 2015 Sébastien Villemot <sebastien@debian.org> + +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 "defun.h" +#include "error.h" +#include "oct-obj.h" +#include "f77-fcn.h" + +extern "C" +{ + F77_RET_T + F77_FUNC (dtrsen, DTRSEN) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type*, const octave_idx_type&, + double*, const octave_idx_type&, double*, const octave_idx_type&, + double*, double*, octave_idx_type&, double&, double&, double*, + const octave_idx_type&, octave_idx_type*, + const octave_idx_type&, octave_idx_type&); + + F77_RET_T + F77_FUNC (ztrsen, ZTRSEN) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type*, const octave_idx_type&, + Complex*, const octave_idx_type&, Complex*, const octave_idx_type&, + Complex*, octave_idx_type&, double&, double&, Complex*, + const octave_idx_type&, octave_idx_type &); + + F77_RET_T + F77_FUNC (strsen, STRSEN) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type*, const octave_idx_type&, + float*, const octave_idx_type&, float*, const octave_idx_type&, + float*, float*, octave_idx_type&, float&, float&, float*, + const octave_idx_type&, octave_idx_type*, + const octave_idx_type&, octave_idx_type&); + + F77_RET_T + F77_FUNC (ctrsen, CTRSEN) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type*, const octave_idx_type&, + FloatComplex*, const octave_idx_type&, FloatComplex*, const octave_idx_type&, + FloatComplex*, octave_idx_type&, float&, float&, FloatComplex*, + const octave_idx_type&, octave_idx_type &); +} + +DEFUN (ordschur, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {[@var{UR}, @var{SR}] =} ordschur (@var{U}, @var{S}, @var{select})\n\ +Reorders the real Schur factorization (@var{U},@var{S}) obtained with the\n\ +@code{schur} function, so that selected eigenvalues appear in the upper left\n\ +diagonal blocks of the quasi triangular Schur matrix.\n\ +The logical vector @var{select} specifies the selected eigenvalues as they\n\ +appear along @var{S}'s diagonal.\n\ +\n\ +For example, given the matrix @code{@var{A} = [1, 2; 3, 4]}, and its Schur\n\ +decomposition\n\ +\n\ +@example\n\ +[@var{U}, @var{S}] = schur (@var{A})\n\ +@end example\n\ +\n\ +which returns\n\ +\n\ +@example\n\ +@group\n\ +@var{U} =\n\ +\n\ + -0.82456 -0.56577\n\ + 0.56577 -0.82456\n\ +\n\ +@var{S} =\n\ +\n\ + -0.37228 -1.00000\n\ + 0.00000 5.37228\n\ +\n\ +@end group\n\ +@end example\n\ +\n\ +It is possible to reorder the decomposition so that the positive eigenvalue is\n\ +in the upper left corner, by doing:\n\ +\n\ +@example\n\ +[@var{U}, @var{S}] = ordschur (@var{U}, @var{S}, [0,1])\n\ +@end example\n\ +\n\ +@seealso{schur}\n\ +@end deftypefn") +{ + const octave_idx_type nargin = args.length (); + octave_value_list retval; + + if (nargin != 3) + { + print_usage (); + return retval; + } + + const Array<octave_idx_type> sel = args(2).octave_idx_type_vector_value (); + if (error_state) + { + error ("ordschur: SELECT must be an array of integers"); + return retval; + } + const octave_idx_type n = sel.numel (); + + const dim_vector dimU = args(0).dims (); + const dim_vector dimS = args(1).dims (); + if (n != dimU(0)) + { + error ("ordschur: SELECT must have same length as the sides of U and S"); + return retval; + } + else if (n != dimU(0) || n != dimS(0) || n != dimU(1) || n != dimS(1)) + { + error ("ordschur: U and S must be square and of equal sizes"); + return retval; + } + + const bool double_type = args(0).is_double_type () + || args(1).is_double_type (); + const bool complex_type = args(0).is_complex_type () + || args(1).is_complex_type (); + +#define PREPARE_ARGS(TYPE, TYPE_M, TYPE_COND) \ + TYPE ## Matrix U = args(0).TYPE_M ## _value (); \ + TYPE ## Matrix S = args(1).TYPE_M ## _value (); \ + if (error_state) \ + { \ + error ("ordschur: U and S must be real or complex floating point matrices"); \ + return retval; \ + } \ + TYPE ## Matrix w (dim_vector (n, 1)); \ + TYPE ## Matrix work (dim_vector (n, 1)); \ + octave_idx_type m; \ + octave_idx_type info; \ + TYPE_COND cond1, cond2; + +#define PREPARE_OUTPUT()\ + if (info != 0) \ + { \ + error ("ordschur: trsen failed"); \ + return retval; \ + } \ + retval(0) = U; \ + retval(1) = S; \ + + if (double_type) + { + if (complex_type) + { + PREPARE_ARGS (Complex, complex_matrix, double) + + F77_XFCN (ztrsen, ztrsen, + (F77_CONST_CHAR_ARG ("N"), F77_CONST_CHAR_ARG ("V"), + sel.data (), n, S.fortran_vec (), n, U.fortran_vec (), n, + w.fortran_vec (), m, cond1, cond2, work.fortran_vec (), n, + info)); + PREPARE_OUTPUT() + } + else + { + PREPARE_ARGS (, matrix, double) + Matrix wi (dim_vector (n, 1)); + Array<octave_idx_type> iwork (dim_vector (n, 1)); + + F77_XFCN (dtrsen, dtrsen, + (F77_CONST_CHAR_ARG ("N"), F77_CONST_CHAR_ARG ("V"), + sel.data (), n, S.fortran_vec (), n, U.fortran_vec (), n, + w.fortran_vec (), wi.fortran_vec (), m, cond1, cond2, + work.fortran_vec (), n, iwork.fortran_vec (), n, info)); + PREPARE_OUTPUT () + } + } + else + { + if (complex_type) + { + PREPARE_ARGS (FloatComplex, float_complex_matrix, float) + + F77_XFCN (ctrsen, ctrsen, + (F77_CONST_CHAR_ARG ("N"), F77_CONST_CHAR_ARG ("V"), + sel.data (), n, S.fortran_vec (), n, U.fortran_vec (), n, + w.fortran_vec (), m, cond1, cond2, work.fortran_vec (), n, + info)); + PREPARE_OUTPUT () + } + else + { + PREPARE_ARGS (Float, float_matrix, float) + FloatMatrix wi (dim_vector (n, 1)); + Array<octave_idx_type> iwork (dim_vector (n, 1)); + + F77_XFCN (strsen, strsen, + (F77_CONST_CHAR_ARG ("N"), F77_CONST_CHAR_ARG ("V"), + sel.data (), n, S.fortran_vec (), n, U.fortran_vec (), n, + w.fortran_vec (), wi.fortran_vec (), m, cond1, cond2, + work.fortran_vec (), n, iwork.fortran_vec (), n, info)); + PREPARE_OUTPUT () + } + } + +#undef PREPARE_ARGS +#undef PREPARE_OUTPUT + + return retval; +} + +/* + +%!test +%! A = [1, 2, 3, -2; 4, 5, 6, -5 ; 7, 8, 9, -5; 10, 11, 12, 4 ]; +%! [U, T] = schur (A); +%! [US, TS] = ordschur (U, T, [ 0, 0, 1, 1 ]); +%! assert (US*TS*US', A, sqrt (eps)) +%! assert (diag (T)(3:4), diag (TS)(1:2), sqrt (eps)) + +%!test +%! A = [1, 2, 3, -2; 4, 5, 6, -5 ; 7, 8, 9, -5; 10, 11, 12, 4 ]; +%! [U, T] = schur (A); +%! [US, TS] = ordschur (single (U), single (T), [ 0, 0, 1, 1 ]); +%! assert (US*TS*US', A, sqrt (eps ("single"))) +%! assert (diag (T)(3:4), diag (TS)(1:2), sqrt (eps ("single"))) + +%!test +%! A = [1, 2, 3, -2; 4, 5, 6, -5 ; 7, 8, 9, -5; 10, 11, 12, 4+3i ]; +%! [U, T] = schur (A); +%! [US, TS] = ordschur (U, T, [ 0, 0, 1, 1 ]); +%! assert (US*TS*US', A, sqrt (eps)) +%! assert (diag (T)(3:4), diag (TS)(1:2), sqrt (eps)) + +%!test +%! A = [1, 2, 3, -2; 4, 5, 6, -5 ; 7, 8, 9, -5; 10, 11, 12, 4+3i ]; +%! [U, T] = schur (A); +%! [US, TS] = ordschur (single (U), single (T), [ 0, 0, 1, 1 ]); +%! assert (US*TS*US', A, sqrt (eps ("single"))) +%! assert (diag (T)(3:4), diag (TS)(1:2), sqrt (eps ("single"))) + +*/
--- a/libinterp/corefcn/schur.cc Wed Feb 18 09:20:59 2015 -0800 +++ b/libinterp/corefcn/schur.cc Sat Feb 07 21:51:20 2015 +0100 @@ -124,7 +124,7 @@ 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, lu, chol, hess, qr, qz, svd}\n\ +@seealso{rsf2csf, ordschur, lu, chol, hess, qr, qz, svd}\n\ @end deftypefn") { octave_value_list retval;