Mercurial > octave
view libinterp/corefcn/schur.cc @ 21139:538b57866b90
consistently use "typename" intead of "class" in template declarations
* Object.h, QtHandlesUtils.cc, QtHandlesUtils.h, ToolBarButton.cc,
ToolBarButton.h, Cell.h, __lin_interpn__.cc, bitfcns.cc, bsxfun.cc,
cellfun.cc, data.cc, filter.cc, gcd.cc, graphics.cc, help.cc, kron.cc,
lookup.cc, ls-mat5.cc, ls-oct-text.h, lu.cc, max.cc, mgorth.cc,
oct-map.cc, oct-map.h, oct-stream.cc, oct-stream.h, octave-link.h,
pr-output.cc, profiler.h, schur.cc, sparse-xdiv.cc, sparse-xpow.cc,
sqrtm.cc, symtab.h, tril.cc, typecast.cc, variables.cc, xdiv.cc,
zfstream.h, __init_fltk__.cc, __magick_read__.cc, chol.cc, qr.cc,
ov-base-diag.cc, ov-base-diag.h, ov-base-int.cc, ov-base-int.h,
ov-base-mat.cc, ov-base-mat.h, ov-base-scalar.cc, ov-base-scalar.h,
ov-base-sparse.cc, ov-base-sparse.h, ov-base.h, ov-classdef.cc,
ov-int-traits.h, ov-java.h, ov-usr-fcn.h, ov.cc, ov.h,
op-dms-template.cc, oct-parse.in.yy, parse.h, pt-mat.cc, Array-b.cc,
Array.cc, Array.h, CDiagMatrix.h, CMatrix.h, CNDArray.h,
DiagArray2.cc, DiagArray2.h, MArray.cc, MArray.h, MDiagArray2.cc,
MDiagArray2.h, MSparse.cc, MSparse.h, MatrixType.cc, Sparse.cc,
Sparse.h, dDiagMatrix.h, dMatrix.h, dNDArray.h, fCDiagMatrix.h,
fCMatrix.h, fCNDArray.h, fDiagMatrix.h, fMatrix.h, fNDArray.h,
idx-vector.cc, idx-vector.h, intNDArray.cc, intNDArray.h, DET.h,
base-aepbal.h, base-lu.cc, base-lu.h, base-qr.cc, base-qr.h,
bsxfun-defs.cc, eigs-base.cc, lo-mappers.h, lo-specfun.cc,
lo-specfun.h, oct-convn.cc, oct-fftw.cc, oct-norm.cc,
sparse-base-chol.cc, sparse-base-chol.h, sparse-base-lu.cc,
sparse-base-lu.h, sparse-dmsolve.cc, mx-inlines.cc,
action-container.h, base-list.h, lo-traits.h, lo-utils.h,
oct-base64.h, oct-binmap.h, oct-cmplx.h, oct-inttypes.cc,
oct-inttypes.h, oct-locbuf.h, oct-refcount.h, oct-sort.cc, oct-sort.h:
Use "typename" instead of "class" in template declarations.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Sun, 24 Jan 2016 13:50:04 -0500 |
parents | 499b851fbfae |
children | fcac5dbbf9ed |
line wrap: on
line source
/* Copyright (C) 1996-2015 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 "errwarn.h" #include "ovl.h" #include "utils.h" template <typename 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 {} {@var{S} =} schur (@var{A})\n\ @deftypefnx {} {@var{S} =} schur (@var{A}, \"real\")\n\ @deftypefnx {} {@var{S} =} schur (@var{A}, \"complex\")\n\ @deftypefnx {} {@var{S} =} schur (@var{A}, @var{opt})\n\ @deftypefnx {} {[@var{U}, @var{S}] =} schur (@dots{})\n\ @cindex Schur decomposition\n\ Compute the Schur@tie{}decomposition of @var{A}.\n\ \n\ The Schur@tie{}decomposition is defined as\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} is real, then\n\ the real Schur@tie{}decomposition is computed, in which the matrix @var{U}\n\ is orthogonal and @var{S} is block upper triangular with blocks of size at\n\ 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\n\ @qcode{\"complex\"}.\n\ \n\ The eigenvalues are optionally ordered along the diagonal according to the\n\ value of @var{opt}. @code{@var{opt} = \"a\"} indicates that all eigenvalues\n\ with negative real parts should be moved to the leading block of @var{S}\n\ (used in @code{are}), @code{@var{opt} = \"d\"} indicates that all\n\ eigenvalues with magnitude less than one should be moved to the leading\n\ block of @var{S} (used in @code{dare}), and @code{@var{opt} = \"u\"}, the\n\ default, indicates that no ordering of eigenvalues should occur. The\n\ leading @var{k} 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 square\n\ matrix, and has applications in the solution of algebraic Riccati equations\n\ in control (see @code{are} and @code{dare}).\n\ @seealso{rsf2csf, ordschur, lu, chol, hess, qr, qz, svd}\n\ @end deftypefn") { int nargin = args.length (); if (nargin < 1 || nargin > 2 || nargout > 2) print_usage (); octave_value arg = args(0); std::string ord; if (nargin == 2) ord = args(1).xstring_value ("schur: second argument must be a string"); bool force_complex = false; if (ord == "real") { ord = ""; } else if (ord == "complex") { force_complex = true; ord = ""; } 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 '%s'", ord.c_str ()); return ovl (); } } octave_idx_type nr = arg.rows (); octave_idx_type nc = arg.columns (); if (nr != nc) err_square_matrix_required ("schur", "A"); if (! arg.is_numeric_type ()) err_wrong_type_arg ("schur", arg); octave_value_list retval; if (arg.is_single_type ()) { if (! force_complex && arg.is_real_type ()) { FloatMatrix tmp = arg.float_matrix_value (); if (nargout <= 1) { FloatSCHUR result (tmp, ord, false); retval = ovl (result.schur_matrix ()); } else { FloatSCHUR result (tmp, ord, true); retval = ovl (result.unitary_matrix (), result.schur_matrix ()); } } else { FloatComplexMatrix ctmp = arg.float_complex_matrix_value (); if (nargout <= 1) { FloatComplexSCHUR result (ctmp, ord, false); retval = ovl (mark_upper_triangular (result.schur_matrix ())); } else { FloatComplexSCHUR result (ctmp, ord, true); retval = ovl (result.unitary_matrix (), mark_upper_triangular (result.schur_matrix ())); } } } else { if (! force_complex && arg.is_real_type ()) { Matrix tmp = arg.matrix_value (); if (nargout <= 1) { SCHUR result (tmp, ord, false); retval = ovl (result.schur_matrix ()); } else { SCHUR result (tmp, ord, true); retval = ovl (result.unitary_matrix (), result.schur_matrix ()); } } else { ComplexMatrix ctmp = arg.complex_matrix_value (); if (nargout <= 1) { ComplexSCHUR result (ctmp, ord, false); retval = ovl (mark_upper_triangular (result.schur_matrix ())); } else { ComplexSCHUR result (ctmp, ord, true); retval = ovl (result.unitary_matrix (), mark_upper_triangular (result.schur_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"))); %!error schur () %!error schur (1,2,3) %!error [a,b,c] = schur (1) %!error <must be a square matrix> schur ([1, 2, 3; 4, 5, 6]) %!error <wrong type argument 'cell'> schur ({1}) %!warning <incorrect ordered schur argument> schur ([1, 2; 3, 4], "bad_opt"); */ DEFUN (rsf2csf, args, nargout, "-*- texinfo -*-\n\ @deftypefn {} {[@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\ @tcode{@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") { if (args.length () != 2 || nargout > 2) print_usage (); if (! args(0).is_numeric_type ()) err_wrong_type_arg ("rsf2csf", args(0)); if (! args(1).is_numeric_type ()) err_wrong_type_arg ("rsf2csf", args(1)); if (args(0).is_complex_type () || args(1).is_complex_type ()) error ("rsf2csf: UR and TR must be real matrices"); 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 (); FloatComplexSCHUR cs (FloatSCHUR (t, u)); return ovl (cs.unitary_matrix (), cs.schur_matrix ()); } else { Matrix u = args(0).matrix_value (); Matrix t = args(1).matrix_value (); ComplexSCHUR cs (SCHUR (t, u)); return ovl (cs.unitary_matrix (), cs.schur_matrix ()); } } /* %!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); */