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