Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/schur.cc @ 10607:f7501986e42d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
make schur more Matlab compatible
author Jaroslav Hajek <highegg@gmail.com>
date Thu, 06 May 2010 09:49:36 +0200
parents d0ce5e973937
children 9c9e07f5eb1c
comparison
equal deleted inserted replaced
10606:ec34c7acd057 10607:f7501986e42d
39 #include "utils.h" 39 #include "utils.h"
40 40
41 DEFUN_DLD (schur, args, nargout, 41 DEFUN_DLD (schur, args, nargout,
42 "-*- texinfo -*-\n\ 42 "-*- texinfo -*-\n\
43 @deftypefn {Loadable Function} {@var{s} =} schur (@var{a})\n\ 43 @deftypefn {Loadable Function} {@var{s} =} schur (@var{a})\n\
44 @deftypefnx {Loadable Function} {@var{s} =} schur (@var{a}, \"complex\")\n\
44 @deftypefnx {Loadable Function} {[@var{u}, @var{s}] =} schur (@var{a}, @var{opt})\n\ 45 @deftypefnx {Loadable Function} {[@var{u}, @var{s}] =} schur (@var{a}, @var{opt})\n\
45 @cindex Schur decomposition\n\ 46 @cindex Schur decomposition\n\
46 The Schur decomposition is used to compute eigenvalues of a\n\ 47 The Schur decomposition is used to compute eigenvalues of a\n\
47 square matrix, and has applications in the solution of algebraic\n\ 48 square matrix, and has applications in the solution of algebraic\n\
48 Riccati equations in control (see @code{are} and @code{dare}).\n\ 49 Riccati equations in control (see @code{are} and @code{dare}).\n\
146 @end tex\n\ 147 @end tex\n\
147 @ifnottex\n\ 148 @ifnottex\n\
148 @code{s}.\n\ 149 @code{s}.\n\
149 @end ifnottex\n\ 150 @end ifnottex\n\
150 \n\ 151 \n\
152 A complex decomposition may be forced by passing \"complex\" as @var{opt}.\n\
153 \n\
151 The eigenvalues are optionally ordered along the diagonal according to\n\ 154 The eigenvalues are optionally ordered along the diagonal according to\n\
152 the value of @code{opt}. @code{opt = \"a\"} indicates that all\n\ 155 the value of @code{opt}. @code{opt = \"a\"} indicates that all\n\
153 eigenvalues with negative real parts should be moved to the leading\n\ 156 eigenvalues with negative real parts should be moved to the leading\n\
154 block of\n\ 157 block of\n\
155 @tex\n\ 158 @tex\n\
227 error ("schur: expecting string as second argument"); 230 error ("schur: expecting string as second argument");
228 return retval; 231 return retval;
229 } 232 }
230 } 233 }
231 234
232 char ord_char = ord.empty () ? 'U' : ord[0]; 235 bool force_complex = false;
233 236
234 if (ord_char != 'U' && ord_char != 'A' && ord_char != 'D' 237 if (ord == "complex")
235 && ord_char != 'u' && ord_char != 'a' && ord_char != 'd') 238 {
236 { 239 force_complex = true;
237 warning ("schur: incorrect ordered schur argument `%c'", 240 ord = std::string ();
238 ord.c_str ()); 241 }
239 return retval; 242 else
243 {
244 char ord_char = ord.empty () ? 'U' : ord[0];
245
246 if (ord_char != 'U' && ord_char != 'A' && ord_char != 'D'
247 && ord_char != 'u' && ord_char != 'a' && ord_char != 'd')
248 {
249 warning ("schur: incorrect ordered schur argument `%c'",
250 ord.c_str ());
251 return retval;
252 }
240 } 253 }
241 254
242 octave_idx_type nr = arg.rows (); 255 octave_idx_type nr = arg.rows ();
243 octave_idx_type nc = arg.columns (); 256 octave_idx_type nc = arg.columns ();
244 257
245 int arg_is_empty = empty_arg ("schur", nr, nc);
246
247 if (arg_is_empty < 0)
248 return retval;
249 else if (arg_is_empty > 0)
250 return octave_value_list (2, Matrix ());
251
252 if (nr != nc) 258 if (nr != nc)
253 { 259 {
254 gripe_square_matrix_required ("schur"); 260 gripe_square_matrix_required ("schur");
255 return retval; 261 return retval;
256 } 262 }
257 263
258 if (arg.is_single_type ()) 264 if (! arg.is_numeric_type ())
259 { 265 gripe_wrong_type_arg ("schur", arg);
260 if (arg.is_real_type ()) 266 else if (arg.is_single_type ())
267 {
268 if (! force_complex && arg.is_real_type ())
261 { 269 {
262 FloatMatrix tmp = arg.float_matrix_value (); 270 FloatMatrix tmp = arg.float_matrix_value ();
263 271
264 if (! error_state) 272 if (! error_state)
265 { 273 {
274 retval(1) = result.schur_matrix (); 282 retval(1) = result.schur_matrix ();
275 retval(0) = result.unitary_matrix (); 283 retval(0) = result.unitary_matrix ();
276 } 284 }
277 } 285 }
278 } 286 }
279 else if (arg.is_complex_type ()) 287 else
280 { 288 {
281 FloatComplexMatrix ctmp = arg.float_complex_matrix_value (); 289 FloatComplexMatrix ctmp = arg.float_complex_matrix_value ();
282 290
283 if (! error_state) 291 if (! error_state)
284 { 292 {
297 } 305 }
298 } 306 }
299 } 307 }
300 else 308 else
301 { 309 {
302 if (arg.is_real_type ()) 310 if (! force_complex && arg.is_real_type ())
303 { 311 {
304 Matrix tmp = arg.matrix_value (); 312 Matrix tmp = arg.matrix_value ();
305 313
306 if (! error_state) 314 if (! error_state)
307 { 315 {
316 retval(1) = result.schur_matrix (); 324 retval(1) = result.schur_matrix ();
317 retval(0) = result.unitary_matrix (); 325 retval(0) = result.unitary_matrix ();
318 } 326 }
319 } 327 }
320 } 328 }
321 else if (arg.is_complex_type ()) 329 else
322 { 330 {
323 ComplexMatrix ctmp = arg.complex_matrix_value (); 331 ComplexMatrix ctmp = arg.complex_matrix_value ();
324 332
325 if (! error_state) 333 if (! error_state)
326 { 334 {
336 retval(1) = result.schur_matrix (); 344 retval(1) = result.schur_matrix ();
337 retval(0) = result.unitary_matrix (); 345 retval(0) = result.unitary_matrix ();
338 } 346 }
339 } 347 }
340 } 348 }
341 else
342 {
343 gripe_wrong_type_arg ("schur", arg);
344 }
345 } 349 }
346 350
347 return retval; 351 return retval;
348 } 352 }
349 353