GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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