GNU Octave  4.2.1
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
eig.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-2017 John W. Eaton
4 Copyright (C) 2016 Barbara Lócsi
5 
6 This file is part of Octave.
7 
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #if defined (HAVE_CONFIG_H)
25 # include "config.h"
26 #endif
27 
28 #include "defun.h"
29 #include "error.h"
30 #include "errwarn.h"
31 #include "ovl.h"
32 
33 #include "EIG.h"
34 #include "fEIG.h"
35 #include "oct-string.h"
36 
37 DEFUN (eig, args, nargout,
38  doc: /* -*- texinfo -*-
39 @deftypefn {} {@var{lambda} =} eig (@var{A})
40 @deftypefnx {} {@var{lambda} =} eig (@var{A}, @var{B})
41 @deftypefnx {} {[@var{V}, @var{lambda}] =} eig (@var{A})
42 @deftypefnx {} {[@var{V}, @var{lambda}] =} eig (@var{A}, @var{B})
43 @deftypefnx {} {[@var{V}, @var{lambda}, @var{W}] =} eig (@var{A})
44 @deftypefnx {} {[@var{V}, @var{lambda}, @var{W}] =} eig (@var{A}, @var{B})
45 @deftypefnx {} {[@dots{}] =} eig (@var{A}, @var{balanceOption})
46 @deftypefnx {} {[@dots{}] =} eig (@var{A}, @var{B}, @var{algorithm})
47 @deftypefnx {} {[@dots{}] =} eig (@dots{}, @var{eigvalOption})
48 Compute the right eigenvalues(V) and optionally the eigenvectors(lambda) and
49 the left eigenvalues(W) of a matrix or a pair of matrices.
50 
51 The flag @var{balanceOption} can be one of:
52 
53 @table @asis
54 @item @qcode{"balance"}
55 Preliminary balancing is on. (default)
56 
57 @item @qcode{"nobalance"}
58 Disables preliminary balancing.
59 @end table
60 
61 The flag @var{eigvalOption} can be one of:
62 
63 @table @asis
64 @item @qcode{"matrix"}
65 Return the eigenvalues in a diagonal matrix. (default if 2 or 3 outputs
66 are specified)
67 
68 @item @qcode{"vector"}
69 Return the eigenvalues in a column vector. (default if 1 output is
70 specified, e.g. @var{lambda} = eig (@var{A}))
71 @end table
72 
73 The flag @var{algorithm} can be one of:
74 
75 @table @asis
76 @item @qcode{"chol"}
77 Uses the Cholesky factorization of B. (default if A is symmetric (Hermitian)
78 and B is symmetric (Hermitian) positive definite)
79 
80 @item @qcode{"qz"}
81 Uses the QZ algorithm. (When A or B are not symmetric always the
82 QZ algorithm will be used)
83 @end table
84 
85 @multitable @columnfractions .31 .23 .23 .23
86 @headitem @tab no flag @tab chol @tab qz
87 @item both are symmetric
88 @tab @qcode{"chol"}
89 @tab @qcode{"chol"}
90 @tab @qcode{"qz"}
91 @item at least one is not symmetric
92 @tab @qcode{"qz"}
93 @tab @qcode{"qz"}
94 @tab @qcode{"qz"}
95 @end multitable
96 
97 The eigenvalues returned by @code{eig} are not ordered.
98 @seealso{eigs, svd}
99 @end deftypefn */)
100 {
101  int nargin = args.length ();
102 
103  if (nargin > 4 || nargin == 0)
104  print_usage ();
105 
107 
109 
110  arg_a = args(0);
111 
112  if (arg_a.is_empty ())
113  return octave_value_list (2, Matrix ());
114 
115  if (! arg_a.is_float_type ())
116  err_wrong_type_arg ("eig", arg_a);
117 
118  if (arg_a.rows () != arg_a.columns ())
119  err_square_matrix_required ("eig", "A");
120 
121  // determine if it's AEP or GEP
122  bool AEPcase = nargin == 1 || args(1).is_string ();
123 
124  if (! AEPcase)
125  {
126  arg_b = args(1);
127 
128  if (arg_b.is_empty ())
129  return octave_value_list (2, Matrix ());
130 
131  if (! arg_b.is_float_type ())
132  err_wrong_type_arg ("eig", arg_b);
133 
134  if (arg_b.rows () != arg_b.columns ())
135  err_square_matrix_required ("eig", "B");
136  }
137 
138  bool qz_flag = false;
139  bool chol_flag = false;
140  bool balance_flag = false;
141  bool no_balance_flag = false;
142  bool matrix_flag = false;
143  bool vector_flag = false;
144 
145  for (int i = (AEPcase ? 1 : 2); i < args.length (); ++i)
146  {
147  if (! args(i).is_string ())
148  err_wrong_type_arg ("eig", args(i));
149 
150  std::string arg_i = args(i).string_value ();
151  if (octave::string::strcmpi (arg_i, "qz"))
152  qz_flag = true;
153  else if (octave::string::strcmpi (arg_i, "chol"))
154  chol_flag = true;
155  else if (octave::string::strcmpi (arg_i, "balance"))
156  balance_flag = true;
157  else if (octave::string::strcmpi (arg_i, "nobalance"))
158  no_balance_flag = true;
159  else if (octave::string::strcmpi (arg_i, "matrix"))
160  matrix_flag = true;
161  else if (octave::string::strcmpi (arg_i, "vector"))
162  vector_flag = true;
163  else
164  error ("eig: invalid option \"%s\"", arg_i.c_str ());
165  }
166 
167  if (balance_flag && no_balance_flag)
168  error ("eig: \"balance\" and \"nobalance\" options are mutually exclusive");
169  if (vector_flag && matrix_flag)
170  error ("eig: \"vector\" and \"matrix\" options are mutually exclusive");
171  if (qz_flag && chol_flag)
172  error ("eig: \"qz\" and \"chol\" options are mutually exclusive");
173 
174  if (AEPcase)
175  {
176  if (qz_flag)
177  error ("eig: invalid \"qz\" option for algebraic eigenvalue problem");
178  if (chol_flag)
179  error ("eig: invalid \"chol\" option for algebraic eigenvalue problem");
180  }
181  else
182  {
183  if (balance_flag)
184  error ("eig: invalid \"balance\" option for generalized eigenvalue problem");
185  if (no_balance_flag)
186  error ("eig: invalid \"nobalance\" option for generalized eigenvalue problem");
187  }
188 
189  // Default is to balance
190  const bool balance = no_balance_flag ? false : true;
191  const bool force_qz = qz_flag;
192 
193 
194  Matrix tmp_a, tmp_b;
195  ComplexMatrix ctmp_a, ctmp_b;
196  FloatMatrix ftmp_a, ftmp_b;
197  FloatComplexMatrix fctmp_a, fctmp_b;
198 
199  if (arg_a.is_single_type ())
200  {
202  if (AEPcase)
203  {
204  if (arg_a.is_real_type ())
205  {
206  ftmp_a = arg_a.float_matrix_value ();
207 
208  result = FloatEIG (ftmp_a, nargout > 1, nargout > 2, balance);
209  }
210  else
211  {
212  fctmp_a = arg_a.float_complex_matrix_value ();
213 
214  result = FloatEIG (fctmp_a, nargout > 1, nargout > 2, balance);
215  }
216  }
217  else
218  {
219  if (arg_a.is_real_type () && arg_b.is_real_type ())
220  {
221  ftmp_a = arg_a.float_matrix_value ();
222  ftmp_b = arg_b.float_matrix_value ();
223 
224  result = FloatEIG (ftmp_a, ftmp_b, nargout > 1, nargout > 2,
225  force_qz);
226  }
227  else
228  {
229  fctmp_a = arg_a.float_complex_matrix_value ();
230  fctmp_b = arg_b.float_complex_matrix_value ();
231 
232  result = FloatEIG (fctmp_a, fctmp_b, nargout > 1, nargout > 2,
233  force_qz);
234  }
235  }
236 
237  if (nargout == 0 || nargout == 1)
238  {
239  if (matrix_flag)
240  retval = ovl (FloatComplexDiagMatrix (result.eigenvalues ()));
241  else
242  retval = ovl (result.eigenvalues ());
243  }
244  else if (nargout == 2)
245  {
246  if (vector_flag)
247  retval = ovl (result.right_eigenvectors (), result.eigenvalues ());
248  else
249  retval = ovl (result.right_eigenvectors (),
250  FloatComplexDiagMatrix (result.eigenvalues ()));
251  }
252  else
253  {
254  if (vector_flag)
255  retval = ovl (result.right_eigenvectors (),
256  result.eigenvalues (),
257  result.left_eigenvectors ());
258  else
259  retval = ovl (result.right_eigenvectors (),
261  result.left_eigenvectors ());
262  }
263  }
264  else
265  {
266  EIG result;
267 
268  if (AEPcase)
269  {
270  if (arg_a.is_real_type ())
271  {
272  tmp_a = arg_a.matrix_value ();
273 
274  result = EIG (tmp_a, nargout > 1, nargout > 2, balance);
275  }
276  else
277  {
278  ctmp_a = arg_a.complex_matrix_value ();
279 
280  result = EIG (ctmp_a, nargout > 1, nargout > 2, balance);
281  }
282  }
283  else
284  {
285  if (arg_a.is_real_type () && arg_b.is_real_type ())
286  {
287  tmp_a = arg_a.matrix_value ();
288  tmp_b = arg_b.matrix_value ();
289 
290  result = EIG (tmp_a, tmp_b, nargout > 1, nargout > 2, force_qz);
291  }
292  else
293  {
294  ctmp_a = arg_a.complex_matrix_value ();
295  ctmp_b = arg_b.complex_matrix_value ();
296 
297  result = EIG (ctmp_a, ctmp_b, nargout > 1, nargout > 2, force_qz);
298  }
299  }
300 
301  if (nargout == 0 || nargout == 1)
302  {
303  if (matrix_flag)
304  retval = ovl (ComplexDiagMatrix (result.eigenvalues ()));
305  else
306  retval = ovl (result.eigenvalues ());
307  }
308  else if (nargout == 2)
309  {
310  if (vector_flag)
311  retval = ovl (result.right_eigenvectors (), result.eigenvalues ());
312  else
313  retval = ovl (result.right_eigenvectors (),
314  ComplexDiagMatrix (result.eigenvalues ()));
315  }
316  else
317  {
318  if (vector_flag)
319  retval = ovl (result.right_eigenvectors (),
320  result.eigenvalues (),
321  result.left_eigenvectors ());
322  else
323  retval = ovl (result.right_eigenvectors (),
324  ComplexDiagMatrix (result.eigenvalues ()),
325  result.left_eigenvectors ());
326  }
327  }
328 
329  return retval;
330 }
331 
332 /*
333 %!assert (eig ([1, 2; 2, 1]), [-1; 3], sqrt (eps))
334 
335 %!test
336 %! [v, d] = eig ([1, 2; 2, 1]);
337 %! x = 1 / sqrt (2);
338 %! assert (d, [-1, 0; 0, 3], sqrt (eps))
339 %! assert (v, [-x, x; x, x], sqrt (eps))
340 
341 %!test
342 %! [v, d, w] = eig ([1, 2; 2, 1]);
343 %! x = 1 / sqrt (2);
344 %! assert (w, [-x, x; x, x], sqrt (eps))
345 
346 %!test
347 %! [v, d] = eig ([1, 2; 2, 1], "balance");
348 %! x = 1 / sqrt (2);
349 %! assert (d, [-1, 0; 0, 3], sqrt (eps))
350 %! assert (v, [-x, x; x, x], sqrt (eps))
351 
352 %!test
353 %! [v, d, w] = eig ([1, 2; 2, 1], "balance");
354 %! x = 1 / sqrt (2);
355 %! assert (w, [-x, x; x, x], sqrt (eps));
356 
357 %!assert (eig (single ([1, 2; 2, 1])), single ([-1; 3]), sqrt (eps ("single")))
358 
359 %!assert (eig (single ([1, 2; 2, 1]), "balance"),
360 %! single ([-1; 3]), sqrt (eps ("single")))
361 
362 %!test
363 %! [v, d] = eig (single ([1, 2; 2, 1]));
364 %! x = single (1 / sqrt (2));
365 %! assert (d, single ([-1, 0; 0, 3]), sqrt (eps ("single")))
366 %! assert (v, [-x, x; x, x], sqrt (eps ("single")))
367 
368 %!test
369 %! [v, d, w] = eig (single ([1, 2; 2, 1]));
370 %! x = single (1 / sqrt (2));
371 %! assert (w, [-x, x; x, x], sqrt (eps ("single")))
372 
373 %!test
374 %! [v, d] = eig (single ([1, 2; 2, 1]), "balance");
375 %! x = single (1 / sqrt (2));
376 %! assert (d, single ([-1, 0; 0, 3]), sqrt (eps ("single")));
377 %! assert (v, [-x, x; x, x], sqrt (eps ("single")))
378 
379 %!test
380 %! [v, d, w] = eig (single ([1, 2; 2, 1]), "balance");
381 %! x = single (1 / sqrt (2));
382 %! assert (w, [-x, x; x, x], sqrt (eps ("single")))
383 
384 
385 ## If (at least one of) the matrices are non-symmetric,
386 ## regardless the algorithm flag the qz algorithm should be used.
387 ## So the results without algorithm flag, with "qz" and with "chol"
388 ## should be the same.
389 %!function nonsym_chol_2_output (A, B, res = sqrt (eps))
390 %! [v, d] = eig (A, B);
391 %! [v2, d2] = eig (A, B, "qz");
392 %! [v3, d3] = eig (A, B, "chol");
393 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), res)
394 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), res)
395 %! assert (v, v2)
396 %! assert (v, v3)
397 %! assert (d, d2)
398 %! assert (d, d3)
399 %!endfunction
400 
401 %!test nonsym_chol_2_output ([1, 2; -1, 1], [3, 3; 1, 2])
402 %!test nonsym_chol_2_output ([1+3i, 2+3i; 3-8i, 8+3i], [8+i, 3+i; 4-9i, 3+i])
403 %!test nonsym_chol_2_output ([1, 2; 3, 8], [8, 3; 4, 3])
404 
405 %!test nonsym_chol_2_output (single ([1, 2; -1, 1]),
406 %! single ([3, 3; 1, 2]), sqrt (eps ("single")))
407 %!test nonsym_chol_2_output (single ([1+3i, 2+3i; 3-8i, 8+3i]),
408 %! single ([8+i, 3+i; 4-9i, 3+i]),
409 %! sqrt (eps ("single")))
410 
411 %!function nonsym_chol_3_output (A, B, res = sqrt (eps))
412 %! [v, d, w] = eig (A, B);
413 %! [v2, d2, w2] = eig (A, B, "qz");
414 %! [v3, d3, w3] = eig (A, B, "chol");
415 %! wt = w';
416 %! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, res)
417 %! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, res)
418 %! assert (v, v2)
419 %! assert (v, v3)
420 %! assert (d, d2)
421 %! assert (d, d3)
422 %! assert (w, w2)
423 %! assert (w, w3)
424 %!endfunction
425 
426 %!test nonsym_chol_3_output ([1, 2; -1, 1], [3, 3; 1, 2])
427 %!test nonsym_chol_3_output ([1+3i, 2+3i; 3-8i, 8+3i], [8+i, 3+i; 4-9i, 3+i])
428 %!test nonsym_chol_3_output ([1, 2; 3, 8], [8, 3; 4, 3])
429 
430 %!test nonsym_chol_3_output (single ([1, 2; -1, 1]),
431 %! single ([3, 3; 1, 2]), sqrt (eps ("single")))
432 %!test nonsym_chol_3_output (single ([1+3i, 2+3i; 3-8i, 8+3i]),
433 %! single ([8+i, 3+i; 4-9i, 3+i]),
434 %! sqrt (eps ("single")))
435 
436 ## If the matrices are symmetric,
437 ## then the chol method is default.
438 ## So the results without algorithm flag and with "chol" should be the same.
439 %!function sym_chol_2_input (A, B, res = sqrt (eps))
440 %! [v, d] = eig (A, B);
441 %! [v2, d2] = eig (A, B, "chol");
442 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), res)
443 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), res)
444 %! assert (v, v2)
445 %! assert (d, d2)
446 %!endfunction
447 
448 %!test sym_chol_2_input ([1, 2; 2, 1], [3, -2; -2, 3])
449 %!test sym_chol_2_input ([1+3i, 2+i; 2-i, 1+3i], [5+9i, 2+i; 2-i, 5+9i])
450 %!test sym_chol_2_input ([1, 1+i; 1-i, 1], [2, 0; 0, 2])
451 
452 %!test sym_chol_2_input (single ([1, 2; 2, 1]), single ([3, -2; -2, 3]),
453 %! sqrt (eps ("single")))
454 %!test sym_chol_2_input (single ([1+3i, 2+i; 2-i, 1+3i]),
455 %! single ([5+9i, 2+i; 2-i, 5+9i]), sqrt (eps ("single")))
456 %!test sym_chol_2_input (single ([1, 1+i; 1-i, 1]), single ([2, 0; 0, 2]),
457 %! sqrt (eps ("single")))
458 
459 %!function sym_chol_3_input (A, B, res = sqrt (eps))
460 %! [v, d, w] = eig (A, B);
461 %! [v2, d2, w2] = eig (A, B, "chol");
462 %! wt = w';
463 %! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, res)
464 %! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, res)
465 %! assert (v, v2)
466 %! assert (d, d2)
467 %! assert (w, w2)
468 %!endfunction
469 
470 %!test sym_chol_3_input ([1, 2; 2, 1], [3, -2; -2, 3])
471 %!test sym_chol_3_input ([1+3i, 2+i; 2-i, 1+3i], [5+9i, 2+i; 2-i, 5+9i])
472 %!test sym_chol_3_input ([1, 1+i; 1-i, 1], [2, 0; 0, 2])
473 
474 %!test sym_chol_3_input (single ([1, 2; 2, 1]), single ([3, -2; -2, 3]),
475 %! sqrt (eps ("single")))
476 %!test sym_chol_3_input (single ([1+3i, 2+i; 2-i, 1+3i]),
477 %! single ([5+9i, 2+i; 2-i, 5+9i]), sqrt (eps ("single")))
478 %!test sym_chol_3_input (single ([1, 1+i; 1-i, 1]), single ([2, 0; 0, 2]),
479 %! sqrt (eps ("single")))
480 
481 ## "balance" is always default
482 ## so the results with and without "balance" should be the same
483 ## while in this case "nobalance" should produce different result
484 %!test
485 %! A = [3 -2 -0.9 0; -2 4 1 -0; -0 0 -1 0; -0.5 -0.5 0.1 1];
486 %! [V1, D1] = eig (A);
487 %! [V2, D2] = eig (A, "balance");
488 %! [V3, D3] = eig (A, "nobalance");
489 %! assert (V1, V2)
490 %! assert (D1, D2)
491 %! assert (isequal (V2, V3), false)
492 
493 ## Testing the flags in all combination.
494 ## If 2 flags are on, than the result should be the same regardless
495 ## of the flags order.
496 ## option1 represents the first order while option2 represents the other order.
497 ## d and d2 should be a diagonal matrix if "matrix" flag is on while
498 ## these should be column vectors if the "vector" flag is on.
499 %!function test_eig_args (args, options1, options2, testd = @() true)
500 %! [v, d, w] = eig (args{:}, options1{:});
501 %! [v2, d2, w2] = eig (args{:}, options2{:});
502 %! assert (testd (d))
503 %! assert (testd (d2))
504 %! assert (v, v2)
505 %! assert (d, d2)
506 %! assert (w, w2)
507 %!endfunction
508 
509 %!function qz_chol_with_shapes (A, B)
510 %! for shapes = struct ("name", {"vector", "matrix"},
511 %! "test", {@isvector, @isdiag})
512 %! test_eig_args ({A, B}, {"qz", shapes.name},
513 %! {shapes.name, "qz"}, shapes.test);
514 %! test_eig_args ({A, B}, {"chol", shapes.name},
515 %! {shapes.name, "chol"}, shapes.test);
516 %! endfor
517 %!endfunction
518 
519 %!function balance_nobalance_with_shapes (A)
520 %! for shapes = struct ("name", {"vector", "matrix"},
521 %! "test", {@isvector, @isdiag})
522 %! test_eig_args ({A}, {"balance", shapes.name},
523 %! {shapes.name, "balance"}, shapes.test);
524 %! test_eig_args ({A}, {"nobalance", shapes.name},
525 %! {shapes.name, "nobalance"}, shapes.test);
526 %! endfor
527 %!endfunction
528 
529 ## Default return format:
530 ## diagonal matrix if 2 or 3 outputs are specified
531 ## column vector if 1 output is specified
532 %!function test_shapes (args)
533 %! d = eig (args{:});
534 %! assert (isvector(d))
535 %! d2 = eig (args{:}, "vector");
536 %! assert (isvector(d2))
537 %! [v, d3] = eig (args{:});
538 %! assert (isdiag(d3))
539 %! d4 = eig (args{:}, "matrix");
540 %! assert (isdiag(d4))
541 %! [v, d5, w] = eig (args{:});
542 %! assert (isdiag(d5))
543 %! d6 = eig (args{:}, "matrix");
544 %! assert (isdiag(d6))
545 %! assert (d, d2)
546 %! assert (d3, d4)
547 %! assert (d5, d6)
548 %! assert (d, diag(d3))
549 %! assert (d, diag(d5))
550 %!endfunction
551 
552 %!function shapes_AEP (A)
553 %! test_shapes({A});
554 %!endfunction
555 
556 %!function shapes_GEP (A, B)
557 %! test_shapes({A, B});
558 %!endfunction
559 
560 %!test balance_nobalance_with_shapes ([1, 2; 2, 1]);
561 %!test balance_nobalance_with_shapes (single ([1, 2; 2, 1]));
562 
563 %!test shapes_AEP ([1, 2; 2, 1]);
564 %!test shapes_AEP (single ([1, 2; 2, 1]));
565 
566 %!test qz_chol_with_shapes ([1, 1+i; 1-i, 1], [2, 0; 0, 2]);
567 %!test qz_chol_with_shapes ([1, 2; 3, 8], [8, 3; 4, 3]);
568 %!test qz_chol_with_shapes ([1, 2; -1, 1], [3, 3; 1, 2]);
569 
570 %!test qz_chol_with_shapes (single ([1, 1+i; 1-i, 1]), single ([2, 0; 0, 2]));
571 %!test qz_chol_with_shapes (single ([1, 2; 3, 8]), single ([8, 3; 4, 3]));
572 %!test qz_chol_with_shapes (single ([1, 2; -1, 1]), single ([3, 3; 1, 2]));
573 
574 %!test shapes_GEP ([1, 1+i; 1-i, 1], [2, 0; 0, 2]);
575 %!test shapes_GEP ([1, 2; 3, 8], [8, 3; 4, 3]);
576 %!test shapes_GEP ([1, 2; -1, 1], [3, 3; 1, 2]);
577 
578 %!test shapes_GEP (single ([1, 1+i; 1-i, 1]), single ([2, 0; 0, 2]));
579 %!test shapes_GEP (single ([1, 2; 3, 8]), single ([8, 3; 4, 3]));
580 %!test shapes_GEP (single ([1, 2; -1, 1]), single ([3, 3; 1, 2]));
581 
582 %!function chol_qz_accuracy (A, B, is_qz_accurate, is_chol_accurate)
583 %! [V1,D1] = eig (A,B, 'qz');
584 %! [V2,D2] = eig (A,B); #default is chol
585 %! assert (isequal (A*V1,A*V1*D1), is_qz_accurate)
586 %! assert (isequal (A*V2, A*V2*D2), is_chol_accurate)
587 %!endfunction
588 
589 %!test
590 %! minij_100 = gallery('minij',100);
591 %! chol_qz_accuracy (minij_100, minij_100, false, true);
592 
593 %!test
594 %! moler_100 = gallery('moler',100);
595 %! chol_qz_accuracy (moler_100, moler_100, false, true);
596 
597 %!test
598 %! A = diag([10^-16, 10^-15]);
599 %! chol_qz_accuracy (A, A, true, false);
600 
601 %!error eig ()
602 %!error eig (false)
603 %!error eig ([1, 2; 3, 4], [4, 3; 2, 1], 1)
604 
605 %!error <EIG requires same size matrices>
606 %! eig ([1, 2; 3, 4], 2)
607 %!error <must be a square matrix>
608 %! eig ([1, 2; 3, 4; 5, 6])
609 %!error <wrong type argument>
610 %! eig ("abcd")
611 %!error <invalid option "abcd">
612 %! eig ([1 2 ; 2 3], "abcd")
613 %!error <invalid "chol" option for algebraic eigenvalue problem>
614 %! eig ([1 2 ; 2 3], "chol")
615 %!error <invalid "qz" option for algebraic eigenvalue problem>
616 %! eig ([1 2 ; 2 3], "qz")
617 %!error <wrong type argument>
618 %! eig (false, [1 2 ; 2 3])
619 %!error <invalid option "abcd">
620 %! eig ([1 2 ; 2 3], [1 2 ; 2 3], "abcd")
621 %!error <invalid "qz" option for algebraic eigenvalue problem>
622 %! eig ([1 2 ; 2 3], "balance", "qz")
623 %!error <invalid option "abcd">
624 %! eig ([1 2 ; 2 3], [1 2 ; 2 3], "vector", "abcd")
625 %!error <invalid option "abcd">
626 %! eig ([1 2 ; 2 3], "balance", "matrix", "abcd")
627 %!error <"balance" and "nobalance" options are mutually exclusive>
628 %! eig ([1 2 ; 2 3], "balance", "nobalance")
629 %!error <"balance" and "nobalance" options are mutually exclusive>
630 %! eig ([1 2 ; 2 3], "nobalance", "balance")
631 %!error <"vector" and "matrix" options are mutually exclusive>
632 %! eig ([1 2 ; 2 3], "matrix", "vector")
633 %!error <"vector" and "matrix" options are mutually exclusive>
634 %! eig ([1 2 ; 2 3], "vector", "matrix")
635 %!error <"vector" and "matrix" options are mutually exclusive>
636 %! eig ([1 2 ; 2 3], [1 2 ; 2 3], "matrix", "vector")
637 %!error <"vector" and "matrix" options are mutually exclusive>
638 %! eig ([1 2 ; 2 3], [1 2 ; 2 3], "vector", "matrix")
639 %!error <wrong type argument>
640 %! eig ([1 2 ; 2 3], [1 2 ; 2 3], false)
641 %!error <wrong type argument>
642 %! eig ([1 2 ; 2 3], [1 2 ; 2 3], [1 2 ; 2 3])
643 */
bool is_real_type(void) const
Definition: ov.h:667
octave_idx_type rows(void) const
Definition: ov.h:489
FloatComplexMatrix left_eigenvectors(void) const
Definition: fEIG.h:119
FloatComplexMatrix float_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:809
OCTAVE_EXPORT octave_value_list isa nd deftypefn *return ovl(args(0).is_integer_type())
OCTINTERP_API void print_usage(void)
Definition: defun.cc:52
std::string doc
Definition: ov-fcn.h:217
Definition: EIG.h:34
ComplexMatrix left_eigenvectors(void) const
Definition: EIG.h:119
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:46
void error(const char *fmt,...)
Definition: error.cc:570
void err_square_matrix_required(const char *fcn, const char *name)
Definition: errwarn.cc:112
FloatComplexColumnVector eigenvalues(void) const
Definition: fEIG.h:117
ComplexColumnVector eigenvalues(void) const
Definition: EIG.h:117
octave_value arg_b
Definition: sylvester.cc:68
bool is_float_type(void) const
Definition: ov.h:630
JNIEnv void * args
Definition: ov-java.cc:67
octave_idx_type columns(void) const
Definition: ov.h:491
OCTAVE_EXPORT octave_value_list return the number of command line arguments passed to Octave If called with the optional argument the function xample nargout(@histc)
Definition: ov-usr-fcn.cc:935
Definition: fEIG.h:34
int nargin
Definition: graphics.cc:10115
octave_value retval
Definition: data.cc:6294
Definition: dMatrix.h:37
Matrix matrix_value(bool frc_str_conv=false) const
Definition: ov.h:787
void err_wrong_type_arg(const char *name, const char *s)
Definition: errwarn.cc:156
With real return the complex result
Definition: data.cc:3375
bool is_empty(void) const
Definition: ov.h:542
bool strcmpi(const T &str_a, const T &str_b)
True if strings are the same, ignoring case.
Definition: oct-string.cc:129
octave_value arg_a
Definition: sylvester.cc:67
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
ComplexMatrix complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:805
FloatComplexMatrix right_eigenvectors(void) const
Definition: fEIG.h:118
FloatMatrix float_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:790
ComplexMatrix right_eigenvectors(void) const
Definition: EIG.h:118
bool is_single_type(void) const
Definition: ov.h:627
If this string is the system will ring the terminal sometimes it is useful to be able to print the original representation of the string
Definition: utils.cc:854