GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
schur.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-2015 John W. Eaton
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 
27 #include <string>
28 
29 #include "CmplxSCHUR.h"
30 #include "dbleSCHUR.h"
31 #include "fCmplxSCHUR.h"
32 #include "floatSCHUR.h"
33 
34 #include "defun.h"
35 #include "error.h"
36 #include "gripes.h"
37 #include "oct-obj.h"
38 #include "utils.h"
39 
40 template <class Matrix>
41 static octave_value
43 {
44  octave_value retval = a;
45 
46  octave_idx_type n = a.rows ();
47  assert (a.columns () == n);
48 
49  const typename Matrix::element_type zero = typename Matrix::element_type ();
50 
51  for (octave_idx_type i = 0; i < n; i++)
52  if (a(i,i) == zero)
53  return retval;
54 
56 
57  return retval;
58 }
59 
60 DEFUN (schur, args, nargout,
61  "-*- texinfo -*-\n\
62 @deftypefn {Built-in Function} {@var{S} =} schur (@var{A})\n\
63 @deftypefnx {Built-in Function} {@var{S} =} schur (@var{A}, \"real\")\n\
64 @deftypefnx {Built-in Function} {@var{S} =} schur (@var{A}, \"complex\")\n\
65 @deftypefnx {Built-in Function} {@var{S} =} schur (@var{A}, @var{opt})\n\
66 @deftypefnx {Built-in Function} {[@var{U}, @var{S}] =} schur (@dots{})\n\
67 @cindex Schur decomposition\n\
68 Compute the Schur@tie{}decomposition of @var{A}.\n\
69 \n\
70 The Schur@tie{}decomposition is defined as\n\
71 @tex\n\
72 $$\n\
73  S = U^T A U\n\
74 $$\n\
75 @end tex\n\
76 @ifnottex\n\
77 \n\
78 @example\n\
79 @code{@var{S} = @var{U}' * @var{A} * @var{U}}\n\
80 @end example\n\
81 \n\
82 @end ifnottex\n\
83 where @var{U} is a unitary matrix\n\
84 @tex\n\
85 ($U^T U$ is identity)\n\
86 @end tex\n\
87 @ifnottex\n\
88 (@code{@var{U}'* @var{U}} is identity)\n\
89 @end ifnottex\n\
90 and @var{S} is upper triangular. The eigenvalues of @var{A} (and @var{S})\n\
91 are the diagonal elements of @var{S}. If the matrix @var{A} is real, then\n\
92 the real Schur@tie{}decomposition is computed, in which the matrix @var{U}\n\
93 is orthogonal and @var{S} is block upper triangular with blocks of size at\n\
94 most\n\
95 @tex\n\
96 $2 \\times 2$\n\
97 @end tex\n\
98 @ifnottex\n\
99 @code{2 x 2}\n\
100 @end ifnottex\n\
101 along the diagonal. The diagonal elements of @var{S}\n\
102 (or the eigenvalues of the\n\
103 @tex\n\
104 $2 \\times 2$\n\
105 @end tex\n\
106 @ifnottex\n\
107 @code{2 x 2}\n\
108 @end ifnottex\n\
109 blocks, when appropriate) are the eigenvalues of @var{A} and @var{S}.\n\
110 \n\
111 The default for real matrices is a real Schur@tie{}decomposition.\n\
112 A complex decomposition may be forced by passing the flag\n\
113 @qcode{\"complex\"}.\n\
114 \n\
115 The eigenvalues are optionally ordered along the diagonal according to the\n\
116 value of @var{opt}. @code{@var{opt} = \"a\"} indicates that all eigenvalues\n\
117 with negative real parts should be moved to the leading block of @var{S}\n\
118 (used in @code{are}), @code{@var{opt} = \"d\"} indicates that all\n\
119 eigenvalues with magnitude less than one should be moved to the leading\n\
120 block of @var{S} (used in @code{dare}), and @code{@var{opt} = \"u\"}, the\n\
121 default, indicates that no ordering of eigenvalues should occur. The\n\
122 leading @var{k} columns of @var{U} always span the @var{A}-invariant\n\
123 subspace corresponding to the @var{k} leading eigenvalues of @var{S}.\n\
124 \n\
125 The Schur@tie{}decomposition is used to compute eigenvalues of a square\n\
126 matrix, and has applications in the solution of algebraic Riccati equations\n\
127 in control (see @code{are} and @code{dare}).\n\
128 @seealso{rsf2csf, ordschur, lu, chol, hess, qr, qz, svd}\n\
129 @end deftypefn")
130 {
131  octave_value_list retval;
132 
133  int nargin = args.length ();
134 
135  if (nargin < 1 || nargin > 2 || nargout > 2)
136  {
137  print_usage ();
138  return retval;
139  }
140 
141  octave_value arg = args(0);
142 
143  std::string ord;
144 
145  if (nargin == 2)
146  {
147  if (args(1).is_string ())
148  ord = args(1).string_value ();
149  else
150  {
151  error ("schur: second argument must be a string");
152  return retval;
153  }
154  }
155 
156  bool force_complex = false;
157 
158  if (ord == "real")
159  {
160  ord = std::string ();
161  }
162  else if (ord == "complex")
163  {
164  force_complex = true;
165  ord = std::string ();
166  }
167  else
168  {
169  char ord_char = ord.empty () ? 'U' : ord[0];
170 
171  if (ord_char != 'U' && ord_char != 'A' && ord_char != 'D'
172  && ord_char != 'u' && ord_char != 'a' && ord_char != 'd')
173  {
174  warning ("schur: incorrect ordered schur argument '%s'",
175  ord.c_str ());
176  return retval;
177  }
178  }
179 
180  octave_idx_type nr = arg.rows ();
181  octave_idx_type nc = arg.columns ();
182 
183  if (nr != nc)
184  {
186  return retval;
187  }
188 
189  if (! arg.is_numeric_type ())
190  gripe_wrong_type_arg ("schur", arg);
191  else if (arg.is_single_type ())
192  {
193  if (! force_complex && arg.is_real_type ())
194  {
195  FloatMatrix tmp = arg.float_matrix_value ();
196 
197  if (! error_state)
198  {
199  if (nargout == 0 || nargout == 1)
200  {
201  FloatSCHUR result (tmp, ord, false);
202  retval(0) = result.schur_matrix ();
203  }
204  else
205  {
206  FloatSCHUR result (tmp, ord, true);
207  retval(1) = result.schur_matrix ();
208  retval(0) = result.unitary_matrix ();
209  }
210  }
211  }
212  else
213  {
215 
216  if (! error_state)
217  {
218 
219  if (nargout == 0 || nargout == 1)
220  {
221  FloatComplexSCHUR result (ctmp, ord, false);
222  retval(0) = mark_upper_triangular (result.schur_matrix ());
223  }
224  else
225  {
226  FloatComplexSCHUR result (ctmp, ord, true);
227  retval(1) = mark_upper_triangular (result.schur_matrix ());
228  retval(0) = result.unitary_matrix ();
229  }
230  }
231  }
232  }
233  else
234  {
235  if (! force_complex && arg.is_real_type ())
236  {
237  Matrix tmp = arg.matrix_value ();
238 
239  if (! error_state)
240  {
241  if (nargout == 0 || nargout == 1)
242  {
243  SCHUR result (tmp, ord, false);
244  retval(0) = result.schur_matrix ();
245  }
246  else
247  {
248  SCHUR result (tmp, ord, true);
249  retval(1) = result.schur_matrix ();
250  retval(0) = result.unitary_matrix ();
251  }
252  }
253  }
254  else
255  {
256  ComplexMatrix ctmp = arg.complex_matrix_value ();
257 
258  if (! error_state)
259  {
260 
261  if (nargout == 0 || nargout == 1)
262  {
263  ComplexSCHUR result (ctmp, ord, false);
264  retval(0) = mark_upper_triangular (result.schur_matrix ());
265  }
266  else
267  {
268  ComplexSCHUR result (ctmp, ord, true);
269  retval(1) = mark_upper_triangular (result.schur_matrix ());
270  retval(0) = result.unitary_matrix ();
271  }
272  }
273  }
274  }
275 
276  return retval;
277 }
278 
279 /*
280 %!test
281 %! a = [1, 2, 3; 4, 5, 9; 7, 8, 6];
282 %! [u, s] = schur (a);
283 %! assert (u' * a * u, s, sqrt (eps));
284 
285 %!test
286 %! a = single ([1, 2, 3; 4, 5, 9; 7, 8, 6]);
287 %! [u, s] = schur (a);
288 %! assert (u' * a * u, s, sqrt (eps ("single")));
289 
290 %!error schur ()
291 %!error schur (1,2,3)
292 %!error [a,b,c] = schur (1)
293 %!error <argument must be a square matrix> schur ([1, 2, 3; 4, 5, 6])
294 %!error <wrong type argument 'cell'> schur ({1})
295 %!warning <incorrect ordered schur argument> schur ([1, 2; 3, 4], "bad_opt");
296 
297 */
298 
299 DEFUN (rsf2csf, args, nargout,
300  "-*- texinfo -*-\n\
301 @deftypefn {Function File} {[@var{U}, @var{T}] =} rsf2csf (@var{UR}, @var{TR})\n\
302 Convert a real, upper quasi-triangular Schur@tie{}form @var{TR} to a complex,\n\
303 upper triangular Schur@tie{}form @var{T}.\n\
304 \n\
305 Note that the following relations hold:\n\
306 \n\
307 @tex\n\
308 $UR \\cdot TR \\cdot {UR}^T = U T U^{\\dagger}$ and\n\
309 $U^{\\dagger} U$ is the identity matrix I.\n\
310 @end tex\n\
311 @ifnottex\n\
312 @tcode{@var{UR} * @var{TR} * @var{UR}' = @var{U} * @var{T} * @var{U}'} and\n\
313 @code{@var{U}' * @var{U}} is the identity matrix I.\n\
314 @end ifnottex\n\
315 \n\
316 Note also that @var{U} and @var{T} are not unique.\n\
317 @seealso{schur}\n\
318 @end deftypefn")
319 {
320  octave_value_list retval;
321 
322  if (args.length () == 2 && nargout <= 2)
323  {
324  if (! args(0).is_numeric_type ())
325  gripe_wrong_type_arg ("rsf2csf", args(0));
326  else if (! args(1).is_numeric_type ())
327  gripe_wrong_type_arg ("rsf2csf", args(1));
328  else if (args(0).is_complex_type () || args(1).is_complex_type ())
329  error ("rsf2csf: UR and TR must be real matrices");
330  else
331  {
332 
333  if (args(0).is_single_type () || args(1).is_single_type ())
334  {
335  FloatMatrix u = args(0).float_matrix_value ();
336  FloatMatrix t = args(1).float_matrix_value ();
337  if (! error_state)
338  {
339  FloatComplexSCHUR cs (FloatSCHUR (t, u));
340 
341  retval(1) = cs.schur_matrix ();
342  retval(0) = cs.unitary_matrix ();
343  }
344  }
345  else
346  {
347  Matrix u = args(0).matrix_value ();
348  Matrix t = args(1).matrix_value ();
349  if (! error_state)
350  {
351  ComplexSCHUR cs (SCHUR (t, u));
352 
353  retval(1) = cs.schur_matrix ();
354  retval(0) = cs.unitary_matrix ();
355  }
356  }
357  }
358  }
359  else
360  print_usage ();
361 
362  return retval;
363 }
364 
365 /*
366 %!test
367 %! A = [1, 1, 1, 2; 1, 2, 1, 1; 1, 1, 3, 1; -2, 1, 1, 1];
368 %! [u, t] = schur (A);
369 %! [U, T] = rsf2csf (u, t);
370 %! assert (norm (u * t * u' - U * T * U'), 0, 1e-12);
371 %! assert (norm (A - U * T * U'), 0, 1e-12);
372 
373 %!test
374 %! A = rand (10);
375 %! [u, t] = schur (A);
376 %! [U, T] = rsf2csf (u, t);
377 %! assert (norm (tril (T, -1)), 0);
378 %! assert (norm (U * U'), 1, 1e-14);
379 
380 %!test
381 %! A = [0, 1;-1, 0];
382 %! [u, t] = schur (A);
383 %! [U, T] = rsf2csf (u,t);
384 %! assert (U * T * U', A, 1e-14);
385 */
static octave_value mark_upper_triangular(const Matrix &a)
Definition: schur.cc:42
bool is_real_type(void) const
Definition: ov.h:651
void gripe_wrong_type_arg(const char *name, const char *s, bool is_error)
Definition: gripes.cc:135
octave_idx_type rows(void) const
Definition: ov.h:473
FloatComplexMatrix float_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:795
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
octave_idx_type length(void) const
Definition: oct-obj.h:89
bool is_numeric_type(void) const
Definition: ov.h:663
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
void error(const char *fmt,...)
Definition: error.cc:476
ComplexMatrix schur_matrix(void) const
Definition: CmplxSCHUR.h:75
void gripe_square_matrix_required(const char *name)
Definition: gripes.cc:81
FloatComplexMatrix unitary_matrix(void) const
Definition: fCmplxSCHUR.h:76
octave_idx_type rows(void) const
Definition: Array.h:313
octave_idx_type columns(void) const
Definition: ov.h:475
int error_state
Definition: error.cc:101
FloatComplexMatrix schur_matrix(void) const
Definition: fCmplxSCHUR.h:74
Definition: dMatrix.h:35
ComplexMatrix unitary_matrix(void) const
Definition: CmplxSCHUR.h:77
Matrix matrix_value(bool frc_str_conv=false) const
Definition: ov.h:773
double arg(double x)
Definition: lo-mappers.h:37
MatrixType matrix_type(void) const
Definition: ov.h:510
void warning(const char *fmt,...)
Definition: error.cc:681
FloatMatrix schur_matrix(void) const
Definition: floatSCHUR.h:71
FloatMatrix unitary_matrix(void) const
Definition: floatSCHUR.h:73
double element_type
Definition: Array.h:118
ComplexMatrix complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:791
FloatMatrix float_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:776
bool is_single_type(void) const
Definition: ov.h:611
Matrix unitary_matrix(void) const
Definition: dbleSCHUR.h:72
octave_idx_type columns(void) const
Definition: Array.h:322
Matrix schur_matrix(void) const
Definition: dbleSCHUR.h:70