Mercurial > octave
changeset 22305:510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
* libinterp/corefcn/eig.cc: add preliminary balancing option, computation of
left eigenvectors as a third output, choosing among generalized eigenvalue
algorithms (chol or qz), and choosing among return value formats of the
eigenvalues (vector or matrix). Expand documentation for new options and
add several new tests (and remove duplicated code in existing tests).
* liboctave/numeric/EIG.cc, liboctave/numeric/fEIG.cc: change dgeev, zgeev,
sgeev, and cgeev, to dgeevx, zgeevx, sgeevx, and cgeevx respectively which
allow for more control over thr solution process. Add new flags to the
functions to support the new options added to the interpreter's eig.
* liboctave/numeric/EIG.h, liboctave/numeric/fEIG.h: fix function declaration
to include the new options.
| author | Barbara Locsi <locsi.barbara@gmail.com> |
|---|---|
| date | Tue, 16 Aug 2016 02:07:58 +0100 |
| parents | c0cdf1c92086 |
| children | 2cd1afd0f12f |
| files | libinterp/corefcn/__qp__.cc libinterp/corefcn/eig.cc libinterp/corefcn/xpow.cc liboctave/numeric/EIG.cc liboctave/numeric/EIG.h liboctave/numeric/fEIG.cc liboctave/numeric/fEIG.h |
| diffstat | 7 files changed, 1069 insertions(+), 459 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/__qp__.cc Mon Aug 15 19:39:24 2016 -0400 +++ b/libinterp/corefcn/__qp__.cc Tue Aug 16 02:07:58 2016 +0100 @@ -148,7 +148,7 @@ } ColumnVector eigenvalH = real (eigH.eigenvalues ()); - Matrix eigenvecH = real (eigH.eigenvectors ()); + Matrix eigenvecH = real (eigH.right_eigenvectors ()); double minReal = eigenvalH.min (); octave_idx_type indminR = 0; for (octave_idx_type i = 0; i < n; i++) @@ -296,7 +296,7 @@ } ColumnVector eigenvalrH = real (eigrH.eigenvalues ()); - Matrix eigenvecrH = real (eigrH.eigenvectors ()); + Matrix eigenvecrH = real (eigrH.right_eigenvectors ()); double mRrH = eigenvalrH.min (); indminR = 0; for (octave_idx_type i = 0; i < n; i++)
--- a/libinterp/corefcn/eig.cc Mon Aug 15 19:39:24 2016 -0400 +++ b/libinterp/corefcn/eig.cc Tue Aug 16 02:07:58 2016 +0100 @@ -1,6 +1,7 @@ /* Copyright (C) 1996-2015 John W. Eaton +Copyright (C) 2016 Barbara Lócsi This file is part of Octave. @@ -24,14 +25,14 @@ # include "config.h" #endif -#include "EIG.h" -#include "fEIG.h" - #include "defun.h" #include "error.h" #include "errwarn.h" #include "ovl.h" -#include "utils.h" + +#include "EIG.h" +#include "fEIG.h" +#include "oct-string.h" DEFUN (eig, args, nargout, doc: /* -*- texinfo -*- @@ -39,12 +40,59 @@ @deftypefnx {} {@var{lambda} =} eig (@var{A}, @var{B}) @deftypefnx {} {[@var{V}, @var{lambda}] =} eig (@var{A}) @deftypefnx {} {[@var{V}, @var{lambda}] =} eig (@var{A}, @var{B}) -Compute the eigenvalues (and optionally the eigenvectors) of a matrix -or a pair of matrices +@deftypefnx {} {[@var{V}, @var{lambda}, @var{W}] =} eig (@var{A}) +@deftypefnx {} {[@var{V}, @var{lambda}, @var{W}] =} eig (@var{A}, @var{B}) +@deftypefnx {} {[@dots{}] =} eig (@var{A}, @var{balanceOption}) +@deftypefnx {} {[@dots{}] =} eig (@var{A}, @var{B}, @var{algorithm}) +@deftypefnx {} {[@dots{}] =} eig (@dots{}, @var{eigvalOption}) +Compute the right eigenvalues(V) and optionally the eigenvectors(lambda) and +the left eigenvalues(W) of a matrix or a pair of matrices. + +The flag @var{balanceOption} can be one of: + +@table @asis +@item @qcode{"balance"} +Preliminary balancing is on. (default) + +@item @qcode{"nobalance"} +Disables preliminary balancing. +@end table + +The flag @var{eigvalOption} can be one of: + +@table @asis +@item @qcode{"matrix"} +Return the eigenvalues in a diagonal matrix. (default if 2 or 3 outputs +are specified) -The algorithm used depends on whether there are one or two input -matrices, if they are real or complex, and if they are symmetric -(Hermitian if complex) or non-symmetric. +@item @qcode{"vector"} +Return the eigenvalues in a column vector. (default if 1 output is +specified, e.g. @var{lambda} = eig (@var{A})) +@end table + +The flag @var{algorithm} can be one of: + +@table @asis +@item @qcode{"chol"} +Uses the Cholesky factorization of B. (default if A is symmetric (Hermitian) +and B is symmetric (Hermitian) positive definite) + +@item @qcode{"qz"} +Uses the QZ algorithm. (When A or B are not symmetric always the +QZ algorithm will be used) +@end table + +@multitable @columnfractions .31 .23 .23 .23 +@headitem @tab no flag @tab chol @tab qz +@item both are symmetric +@tab @qcode{"chol"} +@tab @qcode{"chol"} +@tab @qcode{"qz"} +@item at least one is not symmetric +@tab @qcode{"qz"} +@tab @qcode{"qz"} +@tab @qcode{"qz"} +@end multitable The eigenvalues returned by @code{eig} are not ordered. @seealso{eigs, svd} @@ -52,50 +100,96 @@ { int nargin = args.length (); - if (nargin > 2 || nargin == 0) + if (nargin > 4 || nargin == 0) print_usage (); octave_value_list retval; octave_value arg_a, arg_b; - octave_idx_type nr_a, nr_b, nc_a, nc_b; - nr_a = nr_b = nc_a = nc_b = 0; + arg_a = args(0); - arg_a = args(0); - nr_a = arg_a.rows (); - nc_a = arg_a.columns (); - - int arg_is_empty = empty_arg ("eig", nr_a, nc_a); - if (arg_is_empty < 0) - return retval; - else if (arg_is_empty > 0) + if (arg_a.is_empty ()) return octave_value_list (2, Matrix ()); - if (! arg_a.is_double_type () && ! arg_a.is_single_type ()) + if (! arg_a.is_float_type ()) err_wrong_type_arg ("eig", arg_a); - if (nargin == 2) - { - arg_b = args(1); - nr_b = arg_b.rows (); - nc_b = arg_b.columns (); + if (arg_a.rows () != arg_a.columns ()) + err_square_matrix_required ("eig", "A"); + + // determine if it's AEP or GEP + bool AEPcase = nargin == 1 || args(1).is_string (); + + if (! AEPcase) + { + arg_b = args(1); + + if (arg_b.is_empty ()) + return octave_value_list (2, Matrix ()); + + if (! arg_b.is_float_type ()) + err_wrong_type_arg ("eig", arg_b); + + if (arg_b.rows () != arg_b.columns ()) + err_square_matrix_required ("eig", "B"); + } - arg_is_empty = empty_arg ("eig", nr_b, nc_b); - if (arg_is_empty < 0) - return retval; - else if (arg_is_empty > 0) - return ovl (2, Matrix ()); + bool qz_flag = false; + bool chol_flag = false; + bool balance_flag = false; + bool no_balance_flag = false; + bool matrix_flag = false; + bool vector_flag = false; + + for (int i = (AEPcase ? 1 : 2); i < args.length (); ++i) + { + if (! args(i).is_string ()) + err_wrong_type_arg ("eig", args(i)); - if (! arg_b.is_single_type () && ! arg_b.is_double_type ()) - err_wrong_type_arg ("eig", arg_b); + std::string arg_i = args(i).string_value (); + if (octave::string::strcmpi (arg_i, "qz")) + qz_flag = true; + else if (octave::string::strcmpi (arg_i, "chol")) + chol_flag = true; + else if (octave::string::strcmpi (arg_i, "balance")) + balance_flag = true; + else if (octave::string::strcmpi (arg_i, "nobalance")) + no_balance_flag = true; + else if (octave::string::strcmpi (arg_i, "matrix")) + matrix_flag = true; + else if (octave::string::strcmpi (arg_i, "vector")) + vector_flag = true; + else + error ("eig: invalid option \"%s\"", arg_i.c_str ()); } - if (nr_a != nc_a) - err_square_matrix_required ("eig", "A"); + if (balance_flag && no_balance_flag) + error ("eig: \"balance\" and \"nobalance\" options are mutually exclusive"); + if (vector_flag && matrix_flag) + error ("eig: \"vector\" and \"matrix\" options are mutually exclusive"); + if (qz_flag && chol_flag) + error ("eig: \"qz\" and \"chol\" options are mutually exclusive"); - if (nargin == 2 && nr_b != nc_b) - err_square_matrix_required ("eig", "B"); + if (AEPcase) + { + if (qz_flag) + error ("eig: invalid \"qz\" option for algebraic eigenvalue problem"); + if (chol_flag) + error ("eig: invalid \"chol\" option for algebraic eigenvalue problem"); + } + else + { + if (balance_flag) + error ("eig: invalid \"balance\" option for generalized eigenvalue problem"); + if (no_balance_flag) + error ("eig: invalid \"nobalance\" option for generalized eigenvalue problem"); + } + + // Default is to balance + const bool balance = no_balance_flag ? false : true; + const bool force_qz = qz_flag; + Matrix tmp_a, tmp_b; ComplexMatrix ctmp_a, ctmp_b; @@ -105,99 +199,130 @@ if (arg_a.is_single_type ()) { FloatEIG result; - - if (nargin == 1) + if (AEPcase) { if (arg_a.is_real_type ()) { ftmp_a = arg_a.float_matrix_value (); - result = FloatEIG (ftmp_a, nargout > 1); + result = FloatEIG (ftmp_a, nargout > 1, nargout > 2, balance); } else { fctmp_a = arg_a.float_complex_matrix_value (); - result = FloatEIG (fctmp_a, nargout > 1); + result = FloatEIG (fctmp_a, nargout > 1, nargout > 2, balance); } } - else if (nargin == 2) + else { if (arg_a.is_real_type () && arg_b.is_real_type ()) { ftmp_a = arg_a.float_matrix_value (); ftmp_b = arg_b.float_matrix_value (); - result = FloatEIG (ftmp_a, ftmp_b, nargout > 1); + result = FloatEIG (ftmp_a, ftmp_b, nargout > 1, nargout > 2, + force_qz); } else { fctmp_a = arg_a.float_complex_matrix_value (); fctmp_b = arg_b.float_complex_matrix_value (); - result = FloatEIG (fctmp_a, fctmp_b, nargout > 1); + result = FloatEIG (fctmp_a, fctmp_b, nargout > 1, nargout > 2, + force_qz); } } if (nargout == 0 || nargout == 1) { - retval = ovl (result.eigenvalues ()); + if (matrix_flag) + retval = ovl (FloatComplexDiagMatrix (result.eigenvalues ())); + else + retval = ovl (result.eigenvalues ()); + } + else if (nargout == 2) + { + if (vector_flag) + retval = ovl (result.right_eigenvectors (), result.eigenvalues ()); + else + retval = ovl (result.right_eigenvectors (), + FloatComplexDiagMatrix (result.eigenvalues ())); } else { - // Blame it on Matlab. - FloatComplexDiagMatrix d (result.eigenvalues ()); - - retval = ovl (result.eigenvectors (), d); + if (vector_flag) + retval = ovl (result.right_eigenvectors (), + result.eigenvalues (), + result.left_eigenvectors ()); + else + retval = ovl (result.right_eigenvectors (), + FloatComplexDiagMatrix (result.eigenvalues ()), + result.left_eigenvectors ()); } } else { EIG result; - if (nargin == 1) + if (AEPcase) { if (arg_a.is_real_type ()) { tmp_a = arg_a.matrix_value (); - result = EIG (tmp_a, nargout > 1); + result = EIG (tmp_a, nargout > 1, nargout > 2, balance); } else { ctmp_a = arg_a.complex_matrix_value (); - result = EIG (ctmp_a, nargout > 1); + result = EIG (ctmp_a, nargout > 1, nargout > 2, balance); } } - else if (nargin == 2) + else { if (arg_a.is_real_type () && arg_b.is_real_type ()) { tmp_a = arg_a.matrix_value (); tmp_b = arg_b.matrix_value (); - result = EIG (tmp_a, tmp_b, nargout > 1); + result = EIG (tmp_a, tmp_b, nargout > 1, nargout > 2, force_qz); } else { ctmp_a = arg_a.complex_matrix_value (); ctmp_b = arg_b.complex_matrix_value (); - result = EIG (ctmp_a, ctmp_b, nargout > 1); + result = EIG (ctmp_a, ctmp_b, nargout > 1, nargout > 2, force_qz); } } if (nargout == 0 || nargout == 1) { - retval = ovl (result.eigenvalues ()); + if (matrix_flag) + retval = ovl (ComplexDiagMatrix (result.eigenvalues ())); + else + retval = ovl (result.eigenvalues ()); + } + else if (nargout == 2) + { + if (vector_flag) + retval = ovl (result.right_eigenvectors (), result.eigenvalues ()); + else + retval = ovl (result.right_eigenvectors (), + ComplexDiagMatrix (result.eigenvalues ())); } else { - // Blame it on Matlab. - ComplexDiagMatrix d (result.eigenvalues ()); - - retval = ovl (result.eigenvectors (), d); + if (vector_flag) + retval = ovl (result.right_eigenvectors (), + result.eigenvalues (), + result.left_eigenvectors ()); + else + retval = ovl (result.right_eigenvectors (), + ComplexDiagMatrix (result.eigenvalues ()), + result.left_eigenvectors ()); } } @@ -210,88 +335,309 @@ %!test %! [v, d] = eig ([1, 2; 2, 1]); %! x = 1 / sqrt (2); -%! assert (d, [-1, 0; 0, 3], sqrt (eps)); -%! assert (v, [-x, x; x, x], sqrt (eps)); +%! assert (d, [-1, 0; 0, 3], sqrt (eps)) +%! assert (v, [-x, x; x, x], sqrt (eps)) + +%!test +%! [v, d, w] = eig ([1, 2; 2, 1]); +%! x = 1 / sqrt (2); +%! assert (w, [-x, x; x, x], sqrt (eps)) + +%!test +%! [v, d] = eig ([1, 2; 2, 1], "balance"); +%! x = 1 / sqrt (2); +%! assert (d, [-1, 0; 0, 3], sqrt (eps)) +%! assert (v, [-x, x; x, x], sqrt (eps)) + +%!test +%! [v, d, w] = eig ([1, 2; 2, 1], "balance"); +%! x = 1 / sqrt (2); +%! assert (w, [-x, x; x, x], sqrt (eps)); %!assert (eig (single ([1, 2; 2, 1])), single ([-1; 3]), sqrt (eps ("single"))) +%!assert (eig (single ([1, 2; 2, 1]), "balance"), +%! single ([-1; 3]), sqrt (eps ("single"))) + %!test %! [v, d] = eig (single ([1, 2; 2, 1])); %! x = single (1 / sqrt (2)); -%! assert (d, single ([-1, 0; 0, 3]), sqrt (eps ("single"))); -%! assert (v, [-x, x; x, x], sqrt (eps ("single"))); - -%!test -%! A = [1, 2; -1, 1]; B = [3, 3; 1, 2]; -%! [v, d] = eig (A, B); -%! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps)); -%! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps)); - -%!test -%! A = single ([1, 2; -1, 1]); B = single ([3, 3; 1, 2]); -%! [v, d] = eig (A, B); -%! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps ("single"))); -%! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps ("single"))); +%! assert (d, single ([-1, 0; 0, 3]), sqrt (eps ("single"))) +%! assert (v, [-x, x; x, x], sqrt (eps ("single"))) %!test -%! A = [1, 2; 2, 1]; B = [3, -2; -2, 3]; -%! [v, d] = eig (A, B); -%! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps)); -%! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps)); +%! [v, d, w] = eig (single ([1, 2; 2, 1])); +%! x = single (1 / sqrt (2)); +%! assert (w, [-x, x; x, x], sqrt (eps ("single"))) %!test -%! A = single ([1, 2; 2, 1]); B = single ([3, -2; -2, 3]); -%! [v, d] = eig (A, B); -%! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps ("single"))); -%! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps ("single"))); - -%!test -%! A = [1+3i, 2+i; 2-i, 1+3i]; B = [5+9i, 2+i; 2-i, 5+9i]; -%! [v, d] = eig (A, B); -%! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps)); -%! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps)); +%! [v, d] = eig (single ([1, 2; 2, 1]), "balance"); +%! x = single (1 / sqrt (2)); +%! assert (d, single ([-1, 0; 0, 3]), sqrt (eps ("single"))); +%! assert (v, [-x, x; x, x], sqrt (eps ("single"))) %!test -%! A = single ([1+3i, 2+i; 2-i, 1+3i]); B = single ([5+9i, 2+i; 2-i, 5+9i]); -%! [v, d] = eig (A, B); -%! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps ("single"))); -%! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps ("single"))); +%! [v, d, w] = eig (single ([1, 2; 2, 1]), "balance"); +%! x = single (1 / sqrt (2)); +%! assert (w, [-x, x; x, x], sqrt (eps ("single"))) + + +## If (at least one of) the matrices are non-symmetric, +## regardless the algorithm flag the qz algorithm should be used. +## So the results without algorithm flag, with "qz" and with "chol" +## should be the same. +%!function nonsym_chol_2_output (A, B, res = sqrt (eps)) +%! [v, d] = eig (A, B); +%! [v2, d2] = eig (A, B, "qz"); +%! [v3, d3] = eig (A, B, "chol"); +%! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), res) +%! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), res) +%! assert (v, v2) +%! assert (v, v3) +%! assert (d, d2) +%! assert (d, d3) +%!endfunction + +%!test nonsym_chol_2_output ([1, 2; -1, 1], [3, 3; 1, 2]) +%!test nonsym_chol_2_output ([1+3i, 2+3i; 3-8i, 8+3i], [8+i, 3+i; 4-9i, 3+i]) +%!test nonsym_chol_2_output ([1, 2; 3, 8], [8, 3; 4, 3]) + +%!test nonsym_chol_2_output (single ([1, 2; -1, 1]), +%! single ([3, 3; 1, 2]), sqrt (eps ("single"))) +%!test nonsym_chol_2_output (single ([1+3i, 2+3i; 3-8i, 8+3i]), +%! single ([8+i, 3+i; 4-9i, 3+i]), +%! sqrt (eps ("single"))) + +%!function nonsym_chol_3_output (A, B, res = sqrt (eps)) +%! [v, d, w] = eig (A, B); +%! [v2, d2, w2] = eig (A, B, "qz"); +%! [v3, d3, w3] = eig (A, B, "chol"); +%! wt = w'; +%! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, res) +%! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, res) +%! assert (v, v2) +%! assert (v, v3) +%! assert (d, d2) +%! assert (d, d3) +%! assert (w, w2) +%! assert (w, w3) +%!endfunction + +%!test nonsym_chol_3_output ([1, 2; -1, 1], [3, 3; 1, 2]) +%!test nonsym_chol_3_output ([1+3i, 2+3i; 3-8i, 8+3i], [8+i, 3+i; 4-9i, 3+i]) +%!test nonsym_chol_3_output ([1, 2; 3, 8], [8, 3; 4, 3]) +%!test nonsym_chol_3_output (single ([1, 2; -1, 1]), +%! single ([3, 3; 1, 2]), sqrt (eps ("single"))) +%!test nonsym_chol_3_output (single ([1+3i, 2+3i; 3-8i, 8+3i]), +%! single ([8+i, 3+i; 4-9i, 3+i]), +%! sqrt (eps ("single"))) + +## If the matrices are symmetric, +## then the chol method is default. +## So the results without algorithm flag and with "chol" should be the same. +%!function sym_chol_2_input (A, B, res = sqrt (eps)) +%! [v, d] = eig (A, B); +%! [v2, d2] = eig (A, B, "chol"); +%! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), res) +%! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), res) +%! assert (v, v2) +%! assert (d, d2) +%!endfunction + +%!test sym_chol_2_input ([1, 2; 2, 1], [3, -2; -2, 3]) +%!test sym_chol_2_input ([1+3i, 2+i; 2-i, 1+3i], [5+9i, 2+i; 2-i, 5+9i]) +%!test sym_chol_2_input ([1, 1+i; 1-i, 1], [2, 0; 0, 2]) + +%!test sym_chol_2_input (single ([1, 2; 2, 1]), single ([3, -2; -2, 3]), +%! sqrt (eps ("single"))) +%!test sym_chol_2_input (single ([1+3i, 2+i; 2-i, 1+3i]), +%! single ([5+9i, 2+i; 2-i, 5+9i]), sqrt (eps ("single"))) +%!test sym_chol_2_input (single ([1, 1+i; 1-i, 1]), single ([2, 0; 0, 2]), +%! sqrt (eps ("single"))) + +%!function sym_chol_3_input (A, B, res = sqrt (eps)) +%! [v, d, w] = eig (A, B); +%! [v2, d2, w2] = eig (A, B, "chol"); +%! wt = w'; +%! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, res) +%! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, res) +%! assert (v, v2) +%! assert (d, d2) +%! assert (w, w2) +%!endfunction + +%!test sym_chol_3_input ([1, 2; 2, 1], [3, -2; -2, 3]) +%!test sym_chol_3_input ([1+3i, 2+i; 2-i, 1+3i], [5+9i, 2+i; 2-i, 5+9i]) +%!test sym_chol_3_input ([1, 1+i; 1-i, 1], [2, 0; 0, 2]) + +%!test sym_chol_3_input (single ([1, 2; 2, 1]), single ([3, -2; -2, 3]), +%! sqrt (eps ("single"))) +%!test sym_chol_3_input (single ([1+3i, 2+i; 2-i, 1+3i]), +%! single ([5+9i, 2+i; 2-i, 5+9i]), sqrt (eps ("single"))) +%!test sym_chol_3_input (single ([1, 1+i; 1-i, 1]), single ([2, 0; 0, 2]), +%! sqrt (eps ("single"))) + +## "balance" is always default +## so the results with and without "balance" should be the same +## while in this case "nobalance" should produce different result %!test -%! A = [1+3i, 2+3i; 3-8i, 8+3i]; B = [8+i, 3+i; 4-9i, 3+i]; -%! [v, d] = eig (A, B); -%! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps)); -%! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps)); +%! A = [3 -2 -0.9 0; -2 4 1 -0; -0 0 -1 0; -0.5 -0.5 0.1 1]; +%! [V1, D1] = eig (A); +%! [V2, D2] = eig (A, "balance"); +%! [V3, D3] = eig (A, "nobalance"); +%! assert (V1, V2) +%! assert (D1, D2) +%! assert (isequal (V2, V3), false) + +## Testing the flags in all combination. +## If 2 flags are on, than the result should be the same regardless +## of the flags order. +## option1 represents the first order while option2 represents the other order. +## d and d2 should be a diagonal matrix if "matrix" flag is on while +## these should be column vectors if the "vector" flag is on. +%!function test_eig_args (args, options1, options2, testd = @() true) +%! [v, d, w] = eig (args{:}, options1{:}); +%! [v2, d2, w2] = eig (args{:}, options2{:}); +%! assert (testd (d)) +%! assert (testd (d2)) +%! assert (v, v2) +%! assert (d, d2) +%! assert (w, w2) +%!endfunction + +%!function qz_chol_with_shapes (A, B) +%! for shapes = struct ("name", {"vector", "matrix"}, +%! "test", {@isvector, @isdiag}) +%! test_eig_args ({A, B}, {"qz", shapes.name}, +%! {shapes.name, "qz"}, shapes.test); +%! test_eig_args ({A, B}, {"chol", shapes.name}, +%! {shapes.name, "chol"}, shapes.test); +%! endfor +%!endfunction + +%!function balance_nobalance_with_shapes (A) +%! for shapes = struct ("name", {"vector", "matrix"}, +%! "test", {@isvector, @isdiag}) +%! test_eig_args ({A}, {"balance", shapes.name}, +%! {shapes.name, "balance"}, shapes.test); +%! test_eig_args ({A}, {"nobalance", shapes.name}, +%! {shapes.name, "nobalance"}, shapes.test); +%! endfor +%!endfunction -%!test -%! A = single ([1+3i, 2+3i; 3-8i, 8+3i]); B = single ([8+i, 3+i; 4-9i, 3+i]); -%! [v, d] = eig (A, B); -%! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps ("single"))); -%! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps ("single"))); +## Default return format: +## diagonal matrix if 2 or 3 outputs are specified +## column vector if 1 output is specified +%!function test_shapes (args) +%! d = eig (args{:}); +%! assert (isvector(d)) +%! d2 = eig (args{:}, "vector"); +%! assert (isvector(d2)) +%! [v, d3] = eig (args{:}); +%! assert (isdiag(d3)) +%! d4 = eig (args{:}, "matrix"); +%! assert (isdiag(d4)) +%! [v, d5, w] = eig (args{:}); +%! assert (isdiag(d5)) +%! d6 = eig (args{:}, "matrix"); +%! assert (isdiag(d6)) +%! assert (d, d2) +%! assert (d3, d4) +%! assert (d5, d6) +%! assert (d, diag(d3)) +%! assert (d, diag(d5)) +%!endfunction + +%!function shapes_AEP (A) +%! test_shapes({A}); +%!endfunction + +%!function shapes_GEP (A, B) +%! test_shapes({A, B}); +%!endfunction + +%!test balance_nobalance_with_shapes ([1, 2; 2, 1]); +%!test balance_nobalance_with_shapes (single ([1, 2; 2, 1])); + +%!test shapes_AEP ([1, 2; 2, 1]); +%!test shapes_AEP (single ([1, 2; 2, 1])); + +%!test qz_chol_with_shapes ([1, 1+i; 1-i, 1], [2, 0; 0, 2]); +%!test qz_chol_with_shapes ([1, 2; 3, 8], [8, 3; 4, 3]); +%!test qz_chol_with_shapes ([1, 2; -1, 1], [3, 3; 1, 2]); + +%!test qz_chol_with_shapes (single ([1, 1+i; 1-i, 1]), single ([2, 0; 0, 2])); +%!test qz_chol_with_shapes (single ([1, 2; 3, 8]), single ([8, 3; 4, 3])); +%!test qz_chol_with_shapes (single ([1, 2; -1, 1]), single ([3, 3; 1, 2])); + +%!test shapes_GEP ([1, 1+i; 1-i, 1], [2, 0; 0, 2]); +%!test shapes_GEP ([1, 2; 3, 8], [8, 3; 4, 3]); +%!test shapes_GEP ([1, 2; -1, 1], [3, 3; 1, 2]); + +%!test shapes_GEP (single ([1, 1+i; 1-i, 1]), single ([2, 0; 0, 2])); +%!test shapes_GEP (single ([1, 2; 3, 8]), single ([8, 3; 4, 3])); +%!test shapes_GEP (single ([1, 2; -1, 1]), single ([3, 3; 1, 2])); + +%!function chol_qz_accuracy (A, B, is_qz_accurate, is_chol_accurate) +%! [V1,D1] = eig (A,B, 'qz'); +%! [V2,D2] = eig (A,B); #default is chol +%! assert (isequal (A*V1,A*V1*D1), is_qz_accurate) +%! assert (isequal (A*V2, A*V2*D2), is_chol_accurate) +%!endfunction %!test -%! A = [1, 2; 3, 8]; B = [8, 3; 4, 3]; -%! [v, d] = eig (A, B); -%! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps)); -%! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps)); +%! minij_100 = gallery('minij',100); +%! chol_qz_accuracy (minij_100, minij_100, false, true); %!test -%! A = [1, 1+i; 1-i, 1]; B = [2, 0; 0, 2]; -%! [v, d] = eig (A, B); -%! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps)); -%! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps)); +%! moler_100 = gallery('moler',100); +%! chol_qz_accuracy (moler_100, moler_100, false, true); %!test -%! A = single ([1, 1+i; 1-i, 1]); B = single ([2, 0; 0, 2]); -%! [v, d] = eig (A, B); -%! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps ("single"))); -%! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps ("single"))); +%! A = diag([10^-16, 10^-15]); +%! chol_qz_accuracy (A, A, true, false); %!error eig () +%!error eig (false) %!error eig ([1, 2; 3, 4], [4, 3; 2, 1], 1) -%!error <EIG requires same size matrices> eig ([1, 2; 3, 4], 2) -%!error <must be a square matrix> eig ([1, 2; 3, 4; 5, 6]) -%!error <wrong type argument> eig ("abcd") -%!error <wrong type argument> eig ([1 2 ; 2 3], "abcd") -%!error <wrong type argument> eig (false, [1 2 ; 2 3]) + +%!error <EIG requires same size matrices> +%! eig ([1, 2; 3, 4], 2) +%!error <must be a square matrix> +%! eig ([1, 2; 3, 4; 5, 6]) +%!error <wrong type argument> +%! eig ("abcd") +%!error <invalid option "abcd"> +%! eig ([1 2 ; 2 3], "abcd") +%!error <invalid "chol" option for algebraic eigenvalue problem> +%! eig ([1 2 ; 2 3], "chol") +%!error <invalid "qz" option for algebraic eigenvalue problem> +%! eig ([1 2 ; 2 3], "qz") +%!error <wrong type argument> +%! eig (false, [1 2 ; 2 3]) +%!error <invalid option "abcd"> +%! eig ([1 2 ; 2 3], [1 2 ; 2 3], "abcd") +%!error <invalid "qz" option for algebraic eigenvalue problem> +%! eig ([1 2 ; 2 3], "balance", "qz") +%!error <invalid option "abcd"> +%! eig ([1 2 ; 2 3], [1 2 ; 2 3], "vector", "abcd") +%!error <invalid option "abcd"> +%! eig ([1 2 ; 2 3], "balance", "matrix", "abcd") +%!error <"balance" and "nobalance" options are mutually exclusive> +%! eig ([1 2 ; 2 3], "balance", "nobalance") +%!error <"balance" and "nobalance" options are mutually exclusive> +%! eig ([1 2 ; 2 3], "nobalance", "balance") +%!error <"vector" and "matrix" options are mutually exclusive> +%! eig ([1 2 ; 2 3], "matrix", "vector") +%!error <"vector" and "matrix" options are mutually exclusive> +%! eig ([1 2 ; 2 3], "vector", "matrix") +%!error <"vector" and "matrix" options are mutually exclusive> +%! eig ([1 2 ; 2 3], [1 2 ; 2 3], "matrix", "vector") +%!error <"vector" and "matrix" options are mutually exclusive> +%! eig ([1 2 ; 2 3], [1 2 ; 2 3], "vector", "matrix") +%!error <wrong type argument> +%! eig ([1 2 ; 2 3], [1 2 ; 2 3], false) +%!error <wrong type argument> +%! eig ([1 2 ; 2 3], [1 2 ; 2 3], [1 2 ; 2 3]) */
--- a/libinterp/corefcn/xpow.cc Mon Aug 15 19:39:24 2016 -0400 +++ b/libinterp/corefcn/xpow.cc Tue Aug 16 02:07:58 2016 +0100 @@ -123,7 +123,7 @@ EIG b_eig (b); ComplexColumnVector lambda (b_eig.eigenvalues ()); - ComplexMatrix Q (b_eig.eigenvectors ()); + ComplexMatrix Q (b_eig.right_eigenvectors ()); for (octave_idx_type i = 0; i < nr; i++) { @@ -174,7 +174,7 @@ try { ComplexColumnVector lambda (b_eig.eigenvalues ()); - ComplexMatrix Q (b_eig.eigenvectors ()); + ComplexMatrix Q (b_eig.right_eigenvectors ()); for (octave_idx_type i = 0; i < nr; i++) { @@ -262,7 +262,7 @@ try { ComplexColumnVector lambda (a_eig.eigenvalues ()); - ComplexMatrix Q (a_eig.eigenvectors ()); + ComplexMatrix Q (a_eig.right_eigenvectors ()); for (octave_idx_type i = 0; i < nr; i++) lambda(i) = std::pow (lambda(i), b); @@ -338,7 +338,7 @@ try { ComplexColumnVector lambda (a_eig.eigenvalues ()); - ComplexMatrix Q (a_eig.eigenvectors ()); + ComplexMatrix Q (a_eig.right_eigenvectors ()); for (octave_idx_type i = 0; i < nr; i++) lambda(i) = std::pow (lambda(i), b); @@ -386,7 +386,7 @@ try { ComplexColumnVector lambda (b_eig.eigenvalues ()); - ComplexMatrix Q (b_eig.eigenvectors ()); + ComplexMatrix Q (b_eig.right_eigenvectors ()); for (octave_idx_type i = 0; i < nr; i++) { @@ -434,7 +434,7 @@ try { ComplexColumnVector lambda (b_eig.eigenvalues ()); - ComplexMatrix Q (b_eig.eigenvectors ()); + ComplexMatrix Q (b_eig.right_eigenvectors ()); for (octave_idx_type i = 0; i < nr; i++) { @@ -522,7 +522,7 @@ try { ComplexColumnVector lambda (a_eig.eigenvalues ()); - ComplexMatrix Q (a_eig.eigenvectors ()); + ComplexMatrix Q (a_eig.right_eigenvectors ()); for (octave_idx_type i = 0; i < nr; i++) lambda(i) = std::pow (lambda(i), b); @@ -557,7 +557,7 @@ try { ComplexColumnVector lambda (a_eig.eigenvalues ()); - ComplexMatrix Q (a_eig.eigenvectors ()); + ComplexMatrix Q (a_eig.right_eigenvectors ()); for (octave_idx_type i = 0; i < nr; i++) lambda(i) = std::pow (lambda(i), b); @@ -1527,7 +1527,7 @@ try { FloatComplexColumnVector lambda (b_eig.eigenvalues ()); - FloatComplexMatrix Q (b_eig.eigenvectors ()); + FloatComplexMatrix Q (b_eig.right_eigenvectors ()); for (octave_idx_type i = 0; i < nr; i++) { @@ -1579,7 +1579,7 @@ try { FloatComplexColumnVector lambda (b_eig.eigenvalues ()); - FloatComplexMatrix Q (b_eig.eigenvectors ()); + FloatComplexMatrix Q (b_eig.right_eigenvectors ()); for (octave_idx_type i = 0; i < nr; i++) { @@ -1667,7 +1667,7 @@ try { FloatComplexColumnVector lambda (a_eig.eigenvalues ()); - FloatComplexMatrix Q (a_eig.eigenvectors ()); + FloatComplexMatrix Q (a_eig.right_eigenvectors ()); for (octave_idx_type i = 0; i < nr; i++) lambda(i) = std::pow (lambda(i), b); @@ -1733,7 +1733,7 @@ try { FloatComplexColumnVector lambda (a_eig.eigenvalues ()); - FloatComplexMatrix Q (a_eig.eigenvectors ()); + FloatComplexMatrix Q (a_eig.right_eigenvectors ()); for (octave_idx_type i = 0; i < nr; i++) lambda(i) = std::pow (lambda(i), b); @@ -1781,7 +1781,7 @@ try { FloatComplexColumnVector lambda (b_eig.eigenvalues ()); - FloatComplexMatrix Q (b_eig.eigenvectors ()); + FloatComplexMatrix Q (b_eig.right_eigenvectors ()); for (octave_idx_type i = 0; i < nr; i++) { @@ -1829,7 +1829,7 @@ try { FloatComplexColumnVector lambda (b_eig.eigenvalues ()); - FloatComplexMatrix Q (b_eig.eigenvectors ()); + FloatComplexMatrix Q (b_eig.right_eigenvectors ()); for (octave_idx_type i = 0; i < nr; i++) { @@ -1917,7 +1917,7 @@ try { FloatComplexColumnVector lambda (a_eig.eigenvalues ()); - FloatComplexMatrix Q (a_eig.eigenvectors ()); + FloatComplexMatrix Q (a_eig.right_eigenvectors ()); for (octave_idx_type i = 0; i < nr; i++) lambda(i) = std::pow (lambda(i), b); @@ -1952,7 +1952,7 @@ try { FloatComplexColumnVector lambda (a_eig.eigenvalues ()); - FloatComplexMatrix Q (a_eig.eigenvectors ()); + FloatComplexMatrix Q (a_eig.right_eigenvectors ()); for (octave_idx_type i = 0; i < nr; i++) lambda(i) = std::pow (lambda(i), b);
--- a/liboctave/numeric/EIG.cc Mon Aug 15 19:39:24 2016 -0400 +++ b/liboctave/numeric/EIG.cc Tue Aug 16 02:07:58 2016 +0100 @@ -1,6 +1,7 @@ /* Copyright (C) 1994-2015 John W. Eaton +Copyright (C) 2016 Barbara Lócsi This file is part of Octave. @@ -32,26 +33,40 @@ extern "C" { F77_RET_T - F77_FUNC (dgeev, DGEEV) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - const F77_INT&, F77_DBLE*, - const F77_INT&, F77_DBLE*, F77_DBLE*, - F77_DBLE*, const F77_INT&, F77_DBLE*, - const F77_INT&, F77_DBLE*, - const F77_INT&, F77_INT& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); + F77_FUNC (dgeevx, DGEEVX) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const F77_INT&, F77_DBLE*, + const F77_INT&, F77_DBLE*, F77_DBLE*, + F77_DBLE*, const F77_INT&, F77_DBLE*, + const F77_INT&, F77_INT&, + F77_INT&, F77_DBLE*, F77_DBLE&, + F77_DBLE*, F77_DBLE*, F77_DBLE*, + const F77_INT&, F77_INT*, + F77_INT& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); F77_RET_T - F77_FUNC (zgeev, ZGEEV) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - const F77_INT&, F77_DBLE_CMPLX*, - const F77_INT&, F77_DBLE_CMPLX*, - F77_DBLE_CMPLX*, const F77_INT&, F77_DBLE_CMPLX*, - const F77_INT&, F77_DBLE_CMPLX*, - const F77_INT&, F77_DBLE*, F77_INT& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); + F77_FUNC (zgeevx, ZGEEVX) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const F77_INT&, F77_DBLE_CMPLX*, + const F77_INT&, F77_DBLE_CMPLX*, + F77_DBLE_CMPLX*, const F77_INT&, F77_DBLE_CMPLX*, + const F77_INT&, F77_INT&, + F77_INT&, F77_DBLE*, F77_DBLE&, + F77_DBLE*, F77_DBLE*, F77_DBLE_CMPLX*, + const F77_INT&, F77_DBLE*, + F77_INT& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (dsyev, DSYEV) (F77_CONST_CHAR_ARG_DECL, @@ -135,14 +150,14 @@ } octave_idx_type -EIG::init (const Matrix& a, bool calc_ev) +EIG::init (const Matrix& a, bool calc_rev, bool calc_lev, bool balance) { if (a.any_element_is_inf_or_nan ()) (*current_liboctave_error_handler) ("EIG: matrix contains Inf or NaN values"); if (a.is_symmetric ()) - return symmetric_init (a, calc_ev); + return symmetric_init (a, calc_rev, calc_lev); octave_idx_type n = a.rows (); @@ -160,46 +175,77 @@ Array<double> wi (dim_vector (n, 1)); double *pwi = wi.fortran_vec (); - octave_idx_type tnvr = calc_ev ? n : 0; + octave_idx_type tnvr = calc_rev ? n : 0; Matrix vr (tnvr, tnvr); double *pvr = vr.fortran_vec (); + octave_idx_type tnvl = calc_lev ? n : 0; + Matrix vl (tnvl, tnvl); + double *pvl = vl.fortran_vec (); + octave_idx_type lwork = -1; double dummy_work; - double *dummy = 0; - octave_idx_type idummy = 1; + octave_idx_type ilo; + octave_idx_type ihi; + + Array<double> scale (dim_vector (n, 1)); + double *pscale = scale.fortran_vec (); + + double abnrm; + + Array<double> rconde (dim_vector (n, 1)); + double *prconde = rconde.fortran_vec (); + + Array<double> rcondv (dim_vector (n, 1)); + double *prcondv = rcondv.fortran_vec (); - F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), - F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), - n, tmp_data, n, pwr, pwi, dummy, - idummy, pvr, n, &dummy_work, lwork, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); + octave_idx_type dummy_iwork; + + F77_XFCN (dgeevx, DGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1), + F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 ("N", 1), + n, tmp_data, n, pwr, pwi, pvl, + n, pvr, n, ilo, ihi, pscale, + abnrm, prconde, prcondv, &dummy_work, + lwork, &dummy_iwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); if (info != 0) - (*current_liboctave_error_handler) ("dgeev workspace query failed"); + (*current_liboctave_error_handler) ("dgeevx workspace query failed"); lwork = static_cast<octave_idx_type> (dummy_work); Array<double> work (dim_vector (lwork, 1)); double *pwork = work.fortran_vec (); - F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), - F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), - n, tmp_data, n, pwr, pwi, dummy, - idummy, pvr, n, pwork, lwork, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); + F77_XFCN (dgeevx, DGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1), + F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 ("N", 1), + n, tmp_data, n, pwr, pwi, pvl, + n, pvr, n, ilo, ihi, pscale, + abnrm, prconde, prcondv, pwork, + lwork, &dummy_iwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); if (info < 0) - (*current_liboctave_error_handler) ("unrecoverable error in dgeev"); + (*current_liboctave_error_handler) ("unrecoverable error in dgeevx"); if (info > 0) - (*current_liboctave_error_handler) ("dgeev failed to converge"); + (*current_liboctave_error_handler) ("dgeevx failed to converge"); lambda.resize (n); - octave_idx_type nvr = calc_ev ? n : 0; + octave_idx_type nvr = calc_rev ? n : 0; v.resize (nvr, nvr); + octave_idx_type nvl = calc_lev ? n : 0; + w.resize (nvl, nvl); for (octave_idx_type j = 0; j < n; j++) { @@ -208,6 +254,9 @@ lambda.elem (j) = Complex (wr.elem (j)); for (octave_idx_type i = 0; i < nvr; i++) v.elem (i, j) = vr.elem (i, j); + + for (octave_idx_type i = 0; i < nvl; i++) + w.elem (i, j) = vl.elem (i, j); } else { @@ -224,6 +273,14 @@ v.elem (i, j) = Complex (real_part, imag_part); v.elem (i, j+1) = Complex (real_part, -imag_part); } + + for (octave_idx_type i = 0; i < nvl; i++) + { + double real_part = vl.elem (i, j); + double imag_part = vl.elem (i, j+1); + w.elem (i, j) = Complex (real_part, imag_part); + w.elem (i, j+1) = Complex (real_part, -imag_part); + } j++; } } @@ -232,7 +289,7 @@ } octave_idx_type -EIG::symmetric_init (const Matrix& a, bool calc_ev) +EIG::symmetric_init (const Matrix& a, bool calc_rev, bool calc_lev) { octave_idx_type n = a.rows (); @@ -250,7 +307,7 @@ octave_idx_type lwork = -1; double dummy_work; - F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), F77_CONST_CHAR_ARG2 ("U", 1), n, tmp_data, n, pwr, &dummy_work, lwork, info F77_CHAR_ARG_LEN (1) @@ -263,7 +320,7 @@ Array<double> work (dim_vector (lwork, 1)); double *pwork = work.fortran_vec (); - F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), F77_CONST_CHAR_ARG2 ("U", 1), n, tmp_data, n, pwr, pwork, lwork, info F77_CHAR_ARG_LEN (1) @@ -276,20 +333,21 @@ (*current_liboctave_error_handler) ("dsyev failed to converge"); lambda = ComplexColumnVector (wr); - v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix (); + v = calc_rev ? ComplexMatrix (atmp) : ComplexMatrix (); + w = calc_lev ? ComplexMatrix (atmp) : ComplexMatrix (); return info; } octave_idx_type -EIG::init (const ComplexMatrix& a, bool calc_ev) +EIG::init (const ComplexMatrix& a, bool calc_rev, bool calc_lev, bool balance) { if (a.any_element_is_inf_or_nan ()) (*current_liboctave_error_handler) ("EIG: matrix contains Inf or NaN values"); if (a.is_hermitian ()) - return hermitian_init (a, calc_ev); + return hermitian_init (a, calc_rev, calc_lev); octave_idx_type n = a.rows (); @@ -301,12 +359,16 @@ ComplexMatrix atmp = a; Complex *tmp_data = atmp.fortran_vec (); - ComplexColumnVector w (n); - Complex *pw = w.fortran_vec (); + ComplexColumnVector wr (n); + Complex *pw = wr.fortran_vec (); - octave_idx_type nvr = calc_ev ? n : 0; - ComplexMatrix vtmp (nvr, nvr); - Complex *pv = vtmp.fortran_vec (); + octave_idx_type nvr = calc_rev ? n : 0; + ComplexMatrix vrtmp (nvr, nvr); + Complex *pvr = vrtmp.fortran_vec (); + + octave_idx_type nvl = calc_lev ? n : 0; + ComplexMatrix vltmp (nvl, nvl); + Complex *pvl = vltmp.fortran_vec (); octave_idx_type lwork = -1; Complex dummy_work; @@ -315,44 +377,70 @@ Array<double> rwork (dim_vector (lrwork, 1)); double *prwork = rwork.fortran_vec (); - Complex *dummy = 0; - octave_idx_type idummy = 1; + octave_idx_type ilo; + octave_idx_type ihi; + + Array<double> scale (dim_vector (n, 1)); + double *pscale = scale.fortran_vec (); + + double abnrm; + + Array<double> rconde (dim_vector (n, 1)); + double *prconde = rconde.fortran_vec (); + + Array<double> rcondv (dim_vector (n, 1)); + double *prcondv = rcondv.fortran_vec (); - F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), - F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), - n, F77_DBLE_CMPLX_ARG (tmp_data), n, F77_DBLE_CMPLX_ARG (pw), F77_DBLE_CMPLX_ARG (dummy), idummy, - F77_DBLE_CMPLX_ARG (pv), n, F77_DBLE_CMPLX_ARG (&dummy_work), lwork, prwork, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); + F77_XFCN (zgeevx, ZGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1), + F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 ("N", 1), + n, F77_DBLE_CMPLX_ARG (tmp_data), n, + F77_DBLE_CMPLX_ARG (pw), F77_DBLE_CMPLX_ARG (pvl), n, + F77_DBLE_CMPLX_ARG (pvr), n, ilo, ihi, pscale, abnrm, + prconde, prcondv, + F77_DBLE_CMPLX_ARG (&dummy_work), lwork, prwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); if (info != 0) - (*current_liboctave_error_handler) ("zgeev workspace query failed"); + (*current_liboctave_error_handler) ("zgeevx workspace query failed"); lwork = static_cast<octave_idx_type> (dummy_work.real ()); Array<Complex> work (dim_vector (lwork, 1)); Complex *pwork = work.fortran_vec (); - F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), - F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), - n, F77_DBLE_CMPLX_ARG (tmp_data), n, F77_DBLE_CMPLX_ARG (pw), F77_DBLE_CMPLX_ARG (dummy), idummy, - F77_DBLE_CMPLX_ARG (pv), n, F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); + F77_XFCN (zgeevx, ZGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1), + F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 ("N", 1), + n, F77_DBLE_CMPLX_ARG (tmp_data), n, + F77_DBLE_CMPLX_ARG (pw), F77_DBLE_CMPLX_ARG (pvl), n, + F77_DBLE_CMPLX_ARG (pvr), n, ilo, ihi, pscale, abnrm, + prconde, prcondv, + F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); if (info < 0) - (*current_liboctave_error_handler) ("unrecoverable error in zgeev"); + (*current_liboctave_error_handler) ("unrecoverable error in zgeevx"); if (info > 0) - (*current_liboctave_error_handler) ("zgeev failed to converge"); + (*current_liboctave_error_handler) ("zgeevx failed to converge"); - lambda = w; - v = vtmp; + lambda = wr; + v = vrtmp; + w = vltmp; return info; } octave_idx_type -EIG::hermitian_init (const ComplexMatrix& a, bool calc_ev) +EIG::hermitian_init (const ComplexMatrix& a, bool calc_rev, bool calc_lev) { octave_idx_type n = a.rows (); @@ -374,9 +462,10 @@ Array<double> rwork (dim_vector (lrwork, 1)); double *prwork = rwork.fortran_vec (); - F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), F77_CONST_CHAR_ARG2 ("U", 1), - n, F77_DBLE_CMPLX_ARG (tmp_data), n, pwr, F77_DBLE_CMPLX_ARG (&dummy_work), lwork, + n, F77_DBLE_CMPLX_ARG (tmp_data), n, pwr, + F77_DBLE_CMPLX_ARG (&dummy_work), lwork, prwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); @@ -388,9 +477,10 @@ Array<Complex> work (dim_vector (lwork, 1)); Complex *pwork = work.fortran_vec (); - F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), F77_CONST_CHAR_ARG2 ("U", 1), - n, F77_DBLE_CMPLX_ARG (tmp_data), n, pwr, F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, info + n, F77_DBLE_CMPLX_ARG (tmp_data), n, pwr, + F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); @@ -401,13 +491,15 @@ (*current_liboctave_error_handler) ("zheev failed to converge"); lambda = ComplexColumnVector (wr); - v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix (); + v = calc_rev ? ComplexMatrix (atmp) : ComplexMatrix (); + w = calc_lev ? ComplexMatrix (atmp) : ComplexMatrix (); return info; } octave_idx_type -EIG::init (const Matrix& a, const Matrix& b, bool calc_ev) +EIG::init (const Matrix& a, const Matrix& b, bool calc_rev, bool calc_lev, + bool force_qz) { if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ()) (*current_liboctave_error_handler) @@ -427,13 +519,16 @@ Matrix tmp = b; double *tmp_data = tmp.fortran_vec (); - F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), - n, tmp_data, n, - info - F77_CHAR_ARG_LEN (1))); + if (! force_qz) + { + F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), + n, tmp_data, n, + info + F77_CHAR_ARG_LEN (1))); - if (a.is_symmetric () && b.is_symmetric () && info == 0) - return symmetric_init (a, b, calc_ev); + if (a.is_symmetric () && b.is_symmetric () && info == 0) + return symmetric_init (a, b, calc_rev, calc_lev); + } Matrix atmp = a; double *atmp_data = atmp.fortran_vec (); @@ -450,21 +545,22 @@ Array<double> beta (dim_vector (n, 1)); double *pbeta = beta.fortran_vec (); - octave_idx_type tnvr = calc_ev ? n : 0; + octave_idx_type tnvr = calc_rev ? n : 0; Matrix vr (tnvr, tnvr); double *pvr = vr.fortran_vec (); + octave_idx_type tnvl = calc_lev ? n : 0; + Matrix vl (tnvl, tnvl); + double *pvl = vl.fortran_vec (); + octave_idx_type lwork = -1; double dummy_work; - double *dummy = 0; - octave_idx_type idummy = 1; - - F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), - F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), n, atmp_data, n, btmp_data, n, par, pai, pbeta, - dummy, idummy, pvr, n, + pvl, n, pvr, n, &dummy_work, lwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); @@ -476,11 +572,11 @@ Array<double> work (dim_vector (lwork, 1)); double *pwork = work.fortran_vec (); - F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), - F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), n, atmp_data, n, btmp_data, n, par, pai, pbeta, - dummy, idummy, pvr, n, + pvl, n, pvr, n, pwork, lwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); @@ -492,9 +588,12 @@ (*current_liboctave_error_handler) ("dggev failed to converge"); lambda.resize (n); - octave_idx_type nvr = calc_ev ? n : 0; + octave_idx_type nvr = calc_rev ? n : 0; v.resize (nvr, nvr); + octave_idx_type nvl = calc_lev ? n : 0; + w.resize (nvl, nvl); + for (octave_idx_type j = 0; j < n; j++) { if (ai.elem (j) == 0.0) @@ -502,6 +601,8 @@ lambda.elem (j) = Complex (ar.elem (j) / beta.elem (j)); for (octave_idx_type i = 0; i < nvr; i++) v.elem (i, j) = vr.elem (i, j); + for (octave_idx_type i = 0; i < nvl; i++) + w.elem (i, j) = vl.elem (i, j); } else { @@ -520,6 +621,13 @@ v.elem (i, j) = Complex (real_part, imag_part); v.elem (i, j+1) = Complex (real_part, -imag_part); } + for (octave_idx_type i = 0; i < nvl; i++) + { + double real_part = vl.elem (i, j); + double imag_part = vl.elem (i, j+1); + w.elem (i, j) = Complex (real_part, imag_part); + w.elem (i, j+1) = Complex (real_part, -imag_part); + } j++; } } @@ -528,7 +636,8 @@ } octave_idx_type -EIG::symmetric_init (const Matrix& a, const Matrix& b, bool calc_ev) +EIG::symmetric_init (const Matrix& a, const Matrix& b, bool calc_rev, + bool calc_lev) { octave_idx_type n = a.rows (); octave_idx_type nb = b.rows (); @@ -553,7 +662,7 @@ octave_idx_type lwork = -1; double dummy_work; - F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), F77_CONST_CHAR_ARG2 ("U", 1), n, atmp_data, n, btmp_data, n, @@ -568,7 +677,7 @@ Array<double> work (dim_vector (lwork, 1)); double *pwork = work.fortran_vec (); - F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), F77_CONST_CHAR_ARG2 ("U", 1), n, atmp_data, n, btmp_data, n, @@ -583,13 +692,15 @@ (*current_liboctave_error_handler) ("dsygv failed to converge"); lambda = ComplexColumnVector (wr); - v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix (); + v = calc_rev ? ComplexMatrix (atmp) : ComplexMatrix (); + w = calc_lev ? ComplexMatrix (atmp) : ComplexMatrix (); return info; } octave_idx_type -EIG::init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_ev) +EIG::init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_rev, + bool calc_lev, bool force_qz) { if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ()) (*current_liboctave_error_handler) @@ -609,13 +720,16 @@ ComplexMatrix tmp = b; Complex*tmp_data = tmp.fortran_vec (); - F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), - n, F77_DBLE_CMPLX_ARG (tmp_data), n, - info - F77_CHAR_ARG_LEN (1))); + if (! force_qz) + { + F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), + n, F77_DBLE_CMPLX_ARG (tmp_data), n, + info + F77_CHAR_ARG_LEN (1))); - if (a.is_hermitian () && b.is_hermitian () && info == 0) - return hermitian_init (a, b, calc_ev); + if (a.is_hermitian () && b.is_hermitian () && info == 0) + return hermitian_init (a, b, calc_rev, calc_lev); + } ComplexMatrix atmp = a; Complex *atmp_data = atmp.fortran_vec (); @@ -629,9 +743,13 @@ ComplexColumnVector beta (n); Complex *pbeta = beta.fortran_vec (); - octave_idx_type nvr = calc_ev ? n : 0; - ComplexMatrix vtmp (nvr, nvr); - Complex *pv = vtmp.fortran_vec (); + octave_idx_type nvr = calc_rev ? n : 0; + ComplexMatrix vrtmp (nvr, nvr); + Complex *pvr = vrtmp.fortran_vec (); + + octave_idx_type nvl = calc_lev ? n : 0; + ComplexMatrix vltmp (nvl, nvl); + Complex *pvl = vltmp.fortran_vec (); octave_idx_type lwork = -1; Complex dummy_work; @@ -640,14 +758,15 @@ Array<double> rwork (dim_vector (lrwork, 1)); double *prwork = rwork.fortran_vec (); - Complex *dummy = 0; - octave_idx_type idummy = 1; - F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), - F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), - n, F77_DBLE_CMPLX_ARG (atmp_data), n, F77_DBLE_CMPLX_ARG (btmp_data), n, - F77_DBLE_CMPLX_ARG (palpha), F77_DBLE_CMPLX_ARG (pbeta), F77_DBLE_CMPLX_ARG (dummy), idummy, - F77_DBLE_CMPLX_ARG (pv), n, F77_DBLE_CMPLX_ARG (&dummy_work), lwork, prwork, info + F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), + n, F77_DBLE_CMPLX_ARG (atmp_data), n, + F77_DBLE_CMPLX_ARG (btmp_data), n, + F77_DBLE_CMPLX_ARG (palpha), + F77_DBLE_CMPLX_ARG (pbeta), + F77_DBLE_CMPLX_ARG (pvl), n, + F77_DBLE_CMPLX_ARG (pvr), n, + F77_DBLE_CMPLX_ARG (&dummy_work), lwork, prwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); @@ -659,10 +778,14 @@ Complex *pwork = work.fortran_vec (); F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), - F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), - n, F77_DBLE_CMPLX_ARG (atmp_data), n, F77_DBLE_CMPLX_ARG (btmp_data), n, - F77_DBLE_CMPLX_ARG (palpha), F77_DBLE_CMPLX_ARG (pbeta), F77_DBLE_CMPLX_ARG (dummy), idummy, - F77_DBLE_CMPLX_ARG (pv), n, F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, info + F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), + n, F77_DBLE_CMPLX_ARG (atmp_data), n, + F77_DBLE_CMPLX_ARG (btmp_data), n, + F77_DBLE_CMPLX_ARG (palpha), + F77_DBLE_CMPLX_ARG (pbeta), + F77_DBLE_CMPLX_ARG (pvl), n, + F77_DBLE_CMPLX_ARG (pvr), n, + F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); @@ -677,14 +800,15 @@ for (octave_idx_type j = 0; j < n; j++) lambda.elem (j) = alpha.elem (j) / beta.elem (j); - v = vtmp; + v = vrtmp; + w = vltmp; return info; } octave_idx_type EIG::hermitian_init (const ComplexMatrix& a, const ComplexMatrix& b, - bool calc_ev) + bool calc_rev, bool calc_lev) { octave_idx_type n = a.rows (); octave_idx_type nb = b.rows (); @@ -713,7 +837,7 @@ Array<double> rwork (dim_vector (lrwork, 1)); double *prwork = rwork.fortran_vec (); - F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), F77_CONST_CHAR_ARG2 ("U", 1), n, F77_DBLE_CMPLX_ARG (atmp_data), n, F77_DBLE_CMPLX_ARG (btmp_data), n, @@ -729,7 +853,7 @@ Array<Complex> work (dim_vector (lwork, 1)); Complex *pwork = work.fortran_vec (); - F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), F77_CONST_CHAR_ARG2 ("U", 1), n, F77_DBLE_CMPLX_ARG (atmp_data), n, F77_DBLE_CMPLX_ARG (btmp_data), n, @@ -744,7 +868,8 @@ (*current_liboctave_error_handler) ("zhegv failed to converge"); lambda = ComplexColumnVector (wr); - v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix (); + v = calc_rev ? ComplexMatrix (atmp) : ComplexMatrix (); + w = calc_lev ? ComplexMatrix (atmp) : ComplexMatrix (); return info; }
--- a/liboctave/numeric/EIG.h Mon Aug 15 19:39:24 2016 -0400 +++ b/liboctave/numeric/EIG.h Tue Aug 16 02:07:58 2016 +0100 @@ -40,62 +40,66 @@ public: - EIG (void) : lambda (), v () { } + EIG (void) : lambda (), v (), w () { } - EIG (const Matrix& a, bool calc_eigenvectors = true) - : lambda (), v () + EIG (const Matrix& a, bool calc_rev = true, + bool calc_lev = true, bool balance = true) + : lambda (), v (), w () { - init (a, calc_eigenvectors); + init (a, calc_rev, calc_lev, balance); } - EIG (const Matrix& a, octave_idx_type& info, bool calc_eigenvectors = true) - : lambda (), v () + EIG (const Matrix& a, octave_idx_type& info, + bool calc_rev = true, bool calc_lev = true, bool balance = true) + : lambda (), v (), w () { - info = init (a, calc_eigenvectors); + info = init (a, calc_rev, calc_lev, balance); } - EIG (const Matrix& a, const Matrix& b, bool calc_eigenvectors = true) - : lambda (), v () + EIG (const Matrix& a, const Matrix& b, + bool calc_rev = true, bool calc_lev = true, bool force_qz = false) + : lambda (), v (), w () { - init (a, b, calc_eigenvectors); + init (a, b, calc_rev, calc_lev, force_qz); } EIG (const Matrix& a, const Matrix& b, octave_idx_type& info, - bool calc_eigenvectors = true) - : lambda (), v () + bool calc_rev = true, bool calc_lev = true, bool force_qz = false) + : lambda (), v (), w () { - info = init (a, b, calc_eigenvectors); + info = init (a, b, calc_rev, calc_lev, force_qz); } - EIG (const ComplexMatrix& a, bool calc_eigenvectors = true) - : lambda (), v () + EIG (const ComplexMatrix& a, bool calc_rev = true, + bool calc_lev = true, bool balance = true) + : lambda (), v (), w () { - init (a, calc_eigenvectors); + init (a, calc_rev, calc_lev, balance); } EIG (const ComplexMatrix& a, octave_idx_type& info, - bool calc_eigenvectors = true) - : lambda (), v () + bool calc_rev = true, bool calc_lev = true, bool balance = true) + : lambda (), v (), w () { - info = init (a, calc_eigenvectors); + info = init (a, calc_rev, calc_lev, balance); } EIG (const ComplexMatrix& a, const ComplexMatrix& b, - bool calc_eigenvectors = true) - : lambda (), v () + bool calc_rev = true, bool calc_lev = true, bool force_qz = false) + : lambda (), v (), w () { - init (a, b, calc_eigenvectors); + init (a, b, calc_rev, calc_lev, force_qz); } EIG (const ComplexMatrix& a, const ComplexMatrix& b, - octave_idx_type& info, bool calc_eigenvectors = true) - : lambda (), v () + octave_idx_type& info, bool calc_rev = true, bool calc_lev = true, + bool force_qz = false) + : lambda (), v (), w () { - info = init (a, b, calc_eigenvectors); + info = init (a, b, calc_rev, calc_lev, force_qz); } - EIG (const EIG& a) - : lambda (a.lambda), v (a.v) { } + EIG (const EIG& a) : lambda (a.lambda), v (a.v), w (a.w) { } EIG& operator = (const EIG& a) { @@ -103,6 +107,7 @@ { lambda = a.lambda; v = a.v; + w = a.w; } return *this; } @@ -110,8 +115,8 @@ ~EIG (void) { } ComplexColumnVector eigenvalues (void) const { return lambda; } - - ComplexMatrix eigenvectors (void) const { return v; } + ComplexMatrix right_eigenvectors (void) const { return v; } + ComplexMatrix left_eigenvectors (void) const { return w; } friend std::ostream& operator << (std::ostream& os, const EIG& a); @@ -119,28 +124,32 @@ ComplexColumnVector lambda; ComplexMatrix v; + ComplexMatrix w; - octave_idx_type init (const Matrix& a, bool calc_eigenvectors); + octave_idx_type init (const Matrix& a, bool calc_rev, bool calc_lev, + bool balance); octave_idx_type init (const Matrix& a, const Matrix& b, - bool calc_eigenvectors); + bool calc_rev, bool calc_lev, bool force_qz); - octave_idx_type init (const ComplexMatrix& a, bool calc_eigenvectors); + octave_idx_type init (const ComplexMatrix& a, bool calc_rev, + bool calc_lev, bool balance); octave_idx_type init (const ComplexMatrix& a, const ComplexMatrix& b, - bool calc_eigenvectors); + bool calc_rev, bool calc_lev, bool force_qz); - octave_idx_type symmetric_init (const Matrix& a, bool calc_eigenvectors); + octave_idx_type symmetric_init (const Matrix& a, bool calc_rev, + bool calc_lev); octave_idx_type symmetric_init (const Matrix& a, const Matrix& b, - bool calc_eigenvectors); + bool calc_rev, bool calc_lev); octave_idx_type hermitian_init (const ComplexMatrix& a, - bool calc_eigenvectors); + bool calc_rev, bool calc_lev); octave_idx_type hermitian_init (const ComplexMatrix& a, const ComplexMatrix& b, - bool calc_eigenvectors); + bool calc_rev, bool calc_lev); };
--- a/liboctave/numeric/fEIG.cc Mon Aug 15 19:39:24 2016 -0400 +++ b/liboctave/numeric/fEIG.cc Tue Aug 16 02:07:58 2016 +0100 @@ -1,6 +1,7 @@ /* Copyright (C) 1994-2015 John W. Eaton +Copyright (C) 2016 Barbara Lócsi This file is part of Octave. @@ -32,26 +33,38 @@ extern "C" { F77_RET_T - F77_FUNC (sgeev, SGEEV) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - const F77_INT&, F77_REAL*, - const F77_INT&, F77_REAL*, F77_REAL*, F77_REAL*, - const F77_INT&, F77_REAL*, - const F77_INT&, F77_REAL*, - const F77_INT&, F77_INT& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); + F77_FUNC (sgeevx, SGEEVX) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const F77_INT&, F77_REAL*, + const F77_INT&, F77_REAL*, F77_REAL*, F77_REAL*, + const F77_INT&, F77_REAL*, + const F77_INT&, F77_INT&, + F77_INT&, F77_REAL*, F77_REAL&, F77_REAL*, + F77_REAL*, F77_REAL*, const F77_INT&, + F77_INT*, F77_INT& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); F77_RET_T - F77_FUNC (cgeev, CGEEV) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - const F77_INT&, F77_CMPLX*, - const F77_INT&, F77_CMPLX*, F77_CMPLX*, - const F77_INT&, F77_CMPLX*, - const F77_INT&, F77_CMPLX*, - const F77_INT&, F77_REAL*, F77_INT& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); + F77_FUNC (cgeevx, CGEEVX) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const F77_INT&, F77_CMPLX*, + const F77_INT&, F77_CMPLX*, F77_CMPLX*, + const F77_INT&, F77_CMPLX*, + const F77_INT&, F77_INT&, + F77_INT&, F77_REAL*, F77_REAL&, F77_REAL*, + F77_REAL*, F77_CMPLX*, const F77_INT&, + F77_REAL*, F77_INT& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (ssyev, SSYEV) (F77_CONST_CHAR_ARG_DECL, @@ -132,14 +145,14 @@ } octave_idx_type -FloatEIG::init (const FloatMatrix& a, bool calc_ev) +FloatEIG::init (const FloatMatrix& a, bool calc_rev, bool calc_lev, bool balance) { if (a.any_element_is_inf_or_nan ()) (*current_liboctave_error_handler) ("EIG: matrix contains Inf or NaN values"); if (a.is_symmetric ()) - return symmetric_init (a, calc_ev); + return symmetric_init (a, calc_rev, calc_lev); octave_idx_type n = a.rows (); @@ -157,45 +170,75 @@ Array<float> wi (dim_vector (n, 1)); float *pwi = wi.fortran_vec (); - volatile octave_idx_type nvr = calc_ev ? n : 0; + volatile octave_idx_type nvr = calc_rev ? n : 0; FloatMatrix vr (nvr, nvr); float *pvr = vr.fortran_vec (); + volatile octave_idx_type nvl = calc_lev ? n : 0; + FloatMatrix vl (nvl, nvl); + float *pvl = vl.fortran_vec (); + octave_idx_type lwork = -1; float dummy_work; - float *dummy = 0; - octave_idx_type idummy = 1; + octave_idx_type ilo; + octave_idx_type ihi; + + Array<float> scale (dim_vector (n, 1)); + float *pscale = scale.fortran_vec (); + + float abnrm; + + Array<float> rconde (dim_vector (n, 1)); + float *prconde = rconde.fortran_vec (); + + Array<float> rcondv (dim_vector (n, 1)); + float *prcondv = rcondv.fortran_vec (); - F77_XFCN (sgeev, SGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), - F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), - n, tmp_data, n, pwr, pwi, dummy, - idummy, pvr, n, &dummy_work, lwork, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); + octave_idx_type dummy_iwork; + + F77_XFCN (sgeevx, SGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1), + F77_CONST_CHAR_ARG2 ("N", 1), + F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 ("N", 1), + n, tmp_data, n, pwr, pwi, + pvl, n, pvr, n, + ilo, ihi, pscale, abnrm, prconde, prcondv, + &dummy_work, lwork, &dummy_iwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); if (info != 0) - (*current_liboctave_error_handler) ("sgeev workspace query failed"); + (*current_liboctave_error_handler) ("sgeevx workspace query failed"); lwork = static_cast<octave_idx_type> (dummy_work); Array<float> work (dim_vector (lwork, 1)); float *pwork = work.fortran_vec (); - F77_XFCN (sgeev, SGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), - F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), - n, tmp_data, n, pwr, pwi, dummy, - idummy, pvr, n, pwork, lwork, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); + F77_XFCN (sgeevx, SGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1), + F77_CONST_CHAR_ARG2 ("N", 1), + F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 ("N", 1), + n, tmp_data, n, pwr, pwi, + pvl, n, pvr, n, + ilo, ihi, pscale, abnrm, prconde, prcondv, + pwork, lwork, &dummy_iwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); if (info < 0) - (*current_liboctave_error_handler) ("unrecoverable error in sgeev"); + (*current_liboctave_error_handler) ("unrecoverable error in sgeevx"); if (info > 0) - (*current_liboctave_error_handler) ("sgeev failed to converge"); + (*current_liboctave_error_handler) ("sgeevx failed to converge"); lambda.resize (n); v.resize (nvr, nvr); + w.resize (nvl, nvl); for (octave_idx_type j = 0; j < n; j++) { @@ -204,6 +247,9 @@ lambda.elem (j) = FloatComplex (wr.elem (j)); for (octave_idx_type i = 0; i < nvr; i++) v.elem (i, j) = vr.elem (i, j); + + for (octave_idx_type i = 0; i < nvl; i++) + w.elem (i, j) = vl.elem (i, j); } else { @@ -220,6 +266,13 @@ v.elem (i, j) = FloatComplex (real_part, imag_part); v.elem (i, j+1) = FloatComplex (real_part, -imag_part); } + for (octave_idx_type i = 0; i < nvl; i++) + { + float real_part = vl.elem (i, j); + float imag_part = vl.elem (i, j+1); + w.elem (i, j) = FloatComplex (real_part, imag_part); + w.elem (i, j+1) = FloatComplex (real_part, -imag_part); + } j++; } } @@ -228,7 +281,7 @@ } octave_idx_type -FloatEIG::symmetric_init (const FloatMatrix& a, bool calc_ev) +FloatEIG::symmetric_init (const FloatMatrix& a, bool calc_rev, bool calc_lev) { octave_idx_type n = a.rows (); @@ -246,7 +299,7 @@ octave_idx_type lwork = -1; float dummy_work; - F77_XFCN (ssyev, SSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (ssyev, SSYEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), F77_CONST_CHAR_ARG2 ("U", 1), n, tmp_data, n, pwr, &dummy_work, lwork, info F77_CHAR_ARG_LEN (1) @@ -259,7 +312,7 @@ Array<float> work (dim_vector (lwork, 1)); float *pwork = work.fortran_vec (); - F77_XFCN (ssyev, SSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (ssyev, SSYEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), F77_CONST_CHAR_ARG2 ("U", 1), n, tmp_data, n, pwr, pwork, lwork, info F77_CHAR_ARG_LEN (1) @@ -272,20 +325,22 @@ (*current_liboctave_error_handler) ("ssyev failed to converge"); lambda = FloatComplexColumnVector (wr); - v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix (); + v = calc_rev ? FloatComplexMatrix (atmp) : FloatComplexMatrix (); + w = calc_lev ? FloatComplexMatrix (atmp) : FloatComplexMatrix (); return info; } octave_idx_type -FloatEIG::init (const FloatComplexMatrix& a, bool calc_ev) +FloatEIG::init (const FloatComplexMatrix& a, bool calc_rev, bool calc_lev, + bool balance) { if (a.any_element_is_inf_or_nan ()) (*current_liboctave_error_handler) ("EIG: matrix contains Inf or NaN values"); if (a.is_hermitian ()) - return hermitian_init (a, calc_ev); + return hermitian_init (a, calc_rev, calc_lev); octave_idx_type n = a.rows (); @@ -297,12 +352,16 @@ FloatComplexMatrix atmp = a; FloatComplex *tmp_data = atmp.fortran_vec (); - FloatComplexColumnVector w (n); - FloatComplex *pw = w.fortran_vec (); + FloatComplexColumnVector wr (n); + FloatComplex *pw = wr.fortran_vec (); - octave_idx_type nvr = calc_ev ? n : 0; - FloatComplexMatrix vtmp (nvr, nvr); - FloatComplex *pv = vtmp.fortran_vec (); + octave_idx_type nvr = calc_rev ? n : 0; + FloatComplexMatrix vrtmp (nvr, nvr); + FloatComplex *pvr = vrtmp.fortran_vec (); + + octave_idx_type nvl = calc_lev ? n : 0; + FloatComplexMatrix vltmp (nvl, nvl); + FloatComplex *pvl = vltmp.fortran_vec (); octave_idx_type lwork = -1; FloatComplex dummy_work; @@ -311,44 +370,69 @@ Array<float> rwork (dim_vector (lrwork, 1)); float *prwork = rwork.fortran_vec (); - FloatComplex *dummy = 0; - octave_idx_type idummy = 1; + octave_idx_type ilo; + octave_idx_type ihi; + + Array<float> scale (dim_vector (n, 1)); + float *pscale = scale.fortran_vec (); + + float abnrm; + + Array<float> rconde (dim_vector (n, 1)); + float *prconde = rconde.fortran_vec (); - F77_XFCN (cgeev, CGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), - F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), - n, F77_CMPLX_ARG (tmp_data), n, F77_CMPLX_ARG (pw), F77_CMPLX_ARG (dummy), idummy, - F77_CMPLX_ARG (pv), n, F77_CMPLX_ARG (&dummy_work), lwork, prwork, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); + Array<float> rcondv (dim_vector (n, 1)); + float *prcondv = rcondv.fortran_vec (); + + F77_XFCN (cgeevx, CGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1), + F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 ("N", 1), + n, F77_CMPLX_ARG (tmp_data), n, F77_CMPLX_ARG (pw), + F77_CMPLX_ARG (pvl), n, F77_CMPLX_ARG (pvr), n, + ilo, ihi, pscale, abnrm, prconde, prcondv, + F77_CMPLX_ARG (&dummy_work), lwork, prwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); if (info != 0) - (*current_liboctave_error_handler) ("cgeev workspace query failed"); + (*current_liboctave_error_handler) ("cgeevx workspace query failed"); lwork = static_cast<octave_idx_type> (dummy_work.real ()); Array<FloatComplex> work (dim_vector (lwork, 1)); FloatComplex *pwork = work.fortran_vec (); - F77_XFCN (cgeev, CGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), - F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), - n, F77_CMPLX_ARG (tmp_data), n, F77_CMPLX_ARG (pw), F77_CMPLX_ARG (dummy), idummy, - F77_CMPLX_ARG (pv), n, F77_CMPLX_ARG (pwork), lwork, prwork, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); + F77_XFCN (cgeevx, CGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1), + F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 ("N", 1), + n, F77_CMPLX_ARG (tmp_data), n, F77_CMPLX_ARG (pw), + F77_CMPLX_ARG (pvl), n, F77_CMPLX_ARG (pvr), n, + ilo, ihi, pscale, abnrm, prconde, prcondv, + F77_CMPLX_ARG (pwork), lwork, prwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); if (info < 0) - (*current_liboctave_error_handler) ("unrecoverable error in cgeev"); + (*current_liboctave_error_handler) ("unrecoverable error in cgeevx"); if (info > 0) - (*current_liboctave_error_handler) ("cgeev failed to converge"); + (*current_liboctave_error_handler) ("cgeevx failed to converge"); - lambda = w; - v = vtmp; + lambda = wr; + v = vrtmp; + w = vltmp; return info; } octave_idx_type -FloatEIG::hermitian_init (const FloatComplexMatrix& a, bool calc_ev) +FloatEIG::hermitian_init (const FloatComplexMatrix& a, bool calc_rev, + bool calc_lev) { octave_idx_type n = a.rows (); @@ -370,9 +454,10 @@ Array<float> rwork (dim_vector (lrwork, 1)); float *prwork = rwork.fortran_vec (); - F77_XFCN (cheev, CHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (cheev, CHEEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), F77_CONST_CHAR_ARG2 ("U", 1), - n, F77_CMPLX_ARG (tmp_data), n, pwr, F77_CMPLX_ARG (&dummy_work), lwork, + n, F77_CMPLX_ARG (tmp_data), n, pwr, + F77_CMPLX_ARG (&dummy_work), lwork, prwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); @@ -384,9 +469,10 @@ Array<FloatComplex> work (dim_vector (lwork, 1)); FloatComplex *pwork = work.fortran_vec (); - F77_XFCN (cheev, CHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (cheev, CHEEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), F77_CONST_CHAR_ARG2 ("U", 1), - n, F77_CMPLX_ARG (tmp_data), n, pwr, F77_CMPLX_ARG (pwork), lwork, prwork, info + n, F77_CMPLX_ARG (tmp_data), n, pwr, + F77_CMPLX_ARG (pwork), lwork, prwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); @@ -397,13 +483,15 @@ (*current_liboctave_error_handler) ("cheev failed to converge"); lambda = FloatComplexColumnVector (wr); - v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix (); + v = calc_rev ? FloatComplexMatrix (atmp) : FloatComplexMatrix (); + w = calc_lev ? FloatComplexMatrix (atmp) : FloatComplexMatrix (); return info; } octave_idx_type -FloatEIG::init (const FloatMatrix& a, const FloatMatrix& b, bool calc_ev) +FloatEIG::init (const FloatMatrix& a, const FloatMatrix& b, bool calc_rev, + bool calc_lev, bool force_qz) { if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ()) (*current_liboctave_error_handler) @@ -422,14 +510,16 @@ FloatMatrix tmp = b; float *tmp_data = tmp.fortran_vec (); + if (! force_qz) + { + F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), + n, tmp_data, n, + info + F77_CHAR_ARG_LEN (1))); - F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), - n, tmp_data, n, - info - F77_CHAR_ARG_LEN (1))); - - if (a.is_symmetric () && b.is_symmetric () && info == 0) - return symmetric_init (a, b, calc_ev); + if (a.is_symmetric () && b.is_symmetric () && info == 0) + return symmetric_init (a, b, calc_rev, calc_lev); + } FloatMatrix atmp = a; float *atmp_data = atmp.fortran_vec (); @@ -446,21 +536,22 @@ Array<float> beta (dim_vector (n, 1)); float *pbeta = beta.fortran_vec (); - volatile octave_idx_type nvr = calc_ev ? n : 0; + volatile octave_idx_type nvr = calc_rev ? n : 0; FloatMatrix vr (nvr, nvr); float *pvr = vr.fortran_vec (); + volatile octave_idx_type nvl = calc_lev ? n : 0; + FloatMatrix vl (nvl, nvl); + float *pvl = vl.fortran_vec (); + octave_idx_type lwork = -1; float dummy_work; - float *dummy = 0; - octave_idx_type idummy = 1; - - F77_XFCN (sggev, SGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), - F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (sggev, SGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), n, atmp_data, n, btmp_data, n, par, pai, pbeta, - dummy, idummy, pvr, n, + pvl, n, pvr, n, &dummy_work, lwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); @@ -472,11 +563,11 @@ Array<float> work (dim_vector (lwork, 1)); float *pwork = work.fortran_vec (); - F77_XFCN (sggev, SGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), - F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (sggev, SGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), n, atmp_data, n, btmp_data, n, par, pai, pbeta, - dummy, idummy, pvr, n, + pvl, n, pvr, n, pwork, lwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); @@ -489,6 +580,8 @@ lambda.resize (n); v.resize (nvr, nvr); + w.resize (nvl, nvl); + for (octave_idx_type j = 0; j < n; j++) { @@ -497,6 +590,9 @@ lambda.elem (j) = FloatComplex (ar.elem (j) / beta.elem (j)); for (octave_idx_type i = 0; i < nvr; i++) v.elem (i, j) = vr.elem (i, j); + + for (octave_idx_type i = 0; i < nvl; i++) + w.elem (i, j) = vl.elem (i, j); } else { @@ -515,6 +611,13 @@ v.elem (i, j) = FloatComplex (real_part, imag_part); v.elem (i, j+1) = FloatComplex (real_part, -imag_part); } + for (octave_idx_type i = 0; i < nvl; i++) + { + float real_part = vl.elem (i, j); + float imag_part = vl.elem (i, j+1); + w.elem (i, j) = FloatComplex (real_part, imag_part); + w.elem (i, j+1) = FloatComplex (real_part, -imag_part); + } j++; } } @@ -524,7 +627,7 @@ octave_idx_type FloatEIG::symmetric_init (const FloatMatrix& a, const FloatMatrix& b, - bool calc_ev) + bool calc_rev, bool calc_lev) { octave_idx_type n = a.rows (); octave_idx_type nb = b.rows (); @@ -549,7 +652,7 @@ octave_idx_type lwork = -1; float dummy_work; - F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), F77_CONST_CHAR_ARG2 ("U", 1), n, atmp_data, n, btmp_data, n, @@ -564,7 +667,7 @@ Array<float> work (dim_vector (lwork, 1)); float *pwork = work.fortran_vec (); - F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), F77_CONST_CHAR_ARG2 ("U", 1), n, atmp_data, n, btmp_data, n, @@ -579,14 +682,15 @@ (*current_liboctave_error_handler) ("ssygv failed to converge"); lambda = FloatComplexColumnVector (wr); - v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix (); + v = calc_rev ? FloatComplexMatrix (atmp) : FloatComplexMatrix (); + w = calc_lev ? FloatComplexMatrix (atmp) : FloatComplexMatrix (); return info; } octave_idx_type FloatEIG::init (const FloatComplexMatrix& a, const FloatComplexMatrix& b, - bool calc_ev) + bool calc_rev, bool calc_lev, bool force_qz) { if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ()) (*current_liboctave_error_handler) @@ -606,13 +710,16 @@ FloatComplexMatrix tmp = b; FloatComplex *tmp_data = tmp.fortran_vec (); - F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), - n, F77_CMPLX_ARG (tmp_data), n, - info - F77_CHAR_ARG_LEN (1))); + if (! force_qz) + { + F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), + n, F77_CMPLX_ARG (tmp_data), n, + info + F77_CHAR_ARG_LEN (1))); - if (a.is_hermitian () && b.is_hermitian () && info == 0) - return hermitian_init (a, b, calc_ev); + if (a.is_hermitian () && b.is_hermitian () && info == 0) + return hermitian_init (a, b, calc_rev, calc_lev); + } FloatComplexMatrix atmp = a; FloatComplex *atmp_data = atmp.fortran_vec (); @@ -626,9 +733,13 @@ FloatComplexColumnVector beta (n); FloatComplex *pbeta = beta.fortran_vec (); - octave_idx_type nvr = calc_ev ? n : 0; - FloatComplexMatrix vtmp (nvr, nvr); - FloatComplex *pv = vtmp.fortran_vec (); + octave_idx_type nvr = calc_rev ? n : 0; + FloatComplexMatrix vrtmp (nvr, nvr); + FloatComplex *pvr = vrtmp.fortran_vec (); + + octave_idx_type nvl = calc_lev ? n : 0; + FloatComplexMatrix vltmp (nvl, nvl); + FloatComplex *pvl = vltmp.fortran_vec (); octave_idx_type lwork = -1; FloatComplex dummy_work; @@ -637,14 +748,13 @@ Array<float> rwork (dim_vector (lrwork, 1)); float *prwork = rwork.fortran_vec (); - FloatComplex *dummy = 0; - octave_idx_type idummy = 1; - - F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), - F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), - n, F77_CMPLX_ARG (atmp_data), n, F77_CMPLX_ARG (btmp_data), n, - F77_CMPLX_ARG (palpha), F77_CMPLX_ARG (pbeta), F77_CMPLX_ARG (dummy), idummy, - F77_CMPLX_ARG (pv), n, F77_CMPLX_ARG (&dummy_work), lwork, prwork, info + F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), + n, F77_CMPLX_ARG (atmp_data), n, + F77_CMPLX_ARG (btmp_data), n, + F77_CMPLX_ARG (palpha), F77_CMPLX_ARG (pbeta), + F77_CMPLX_ARG (pvl), n, F77_CMPLX_ARG (pvr), n, + F77_CMPLX_ARG (&dummy_work), lwork, prwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); @@ -655,11 +765,13 @@ Array<FloatComplex> work (dim_vector (lwork, 1)); FloatComplex *pwork = work.fortran_vec (); - F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), - F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), - n, F77_CMPLX_ARG (atmp_data), n, F77_CMPLX_ARG (btmp_data), n, - F77_CMPLX_ARG (palpha), F77_CMPLX_ARG (pbeta), F77_CMPLX_ARG (dummy), idummy, - F77_CMPLX_ARG (pv), n, F77_CMPLX_ARG (pwork), lwork, prwork, info + F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), + F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), + n, F77_CMPLX_ARG (atmp_data), n, + F77_CMPLX_ARG (btmp_data), n, + F77_CMPLX_ARG (palpha), F77_CMPLX_ARG (pbeta), + F77_CMPLX_ARG (pvl), n, F77_CMPLX_ARG (pvr), n, + F77_CMPLX_ARG (pwork), lwork, prwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); @@ -674,14 +786,16 @@ for (octave_idx_type j = 0; j < n; j++) lambda.elem (j) = alpha.elem (j) / beta.elem (j); - v = vtmp; + v = vrtmp; + w = vltmp; return info; } octave_idx_type FloatEIG::hermitian_init (const FloatComplexMatrix& a, - const FloatComplexMatrix& b, bool calc_ev) + const FloatComplexMatrix& b, + bool calc_rev, bool calc_lev) { octave_idx_type n = a.rows (); octave_idx_type nb = b.rows (); @@ -710,7 +824,7 @@ Array<float> rwork (dim_vector (lrwork, 1)); float *prwork = rwork.fortran_vec (); - F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), F77_CONST_CHAR_ARG2 ("U", 1), n, F77_CMPLX_ARG (atmp_data), n, F77_CMPLX_ARG (btmp_data), n, @@ -726,7 +840,7 @@ Array<FloatComplex> work (dim_vector (lwork, 1)); FloatComplex *pwork = work.fortran_vec (); - F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), + F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), F77_CONST_CHAR_ARG2 ("U", 1), n, F77_CMPLX_ARG (atmp_data), n, F77_CMPLX_ARG (btmp_data), n, @@ -741,7 +855,8 @@ (*current_liboctave_error_handler) ("zhegv failed to converge"); lambda = FloatComplexColumnVector (wr); - v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix (); + v = calc_rev ? FloatComplexMatrix (atmp) : FloatComplexMatrix (); + w = calc_lev ? FloatComplexMatrix (atmp) : FloatComplexMatrix (); return info; }
--- a/liboctave/numeric/fEIG.h Mon Aug 15 19:39:24 2016 -0400 +++ b/liboctave/numeric/fEIG.h Tue Aug 16 02:07:58 2016 +0100 @@ -40,64 +40,66 @@ public: - FloatEIG (void) - : lambda (), v () { } + FloatEIG (void) : lambda (), v (), w () { } - FloatEIG (const FloatMatrix& a, bool calc_eigenvectors = true) - : lambda (), v () + FloatEIG (const FloatMatrix& a, bool calc_rev = true, + bool calc_lev = true, bool balance = true) + : lambda (), v (), w () { - init (a, calc_eigenvectors); + init (a, calc_rev, calc_lev, balance); } FloatEIG (const FloatMatrix& a, octave_idx_type& info, - bool calc_eigenvectors = true) - : lambda (), v () + bool calc_rev = true, bool calc_lev = true, bool balance = true) + : lambda (), v (), w () { - info = init (a, calc_eigenvectors); + info = init (a, calc_rev, calc_lev, balance); } FloatEIG (const FloatMatrix& a, const FloatMatrix& b, - bool calc_eigenvectors = true) - : lambda (), v () + bool calc_rev = true, bool calc_lev = true, bool force_qz = false) + : lambda (), v (), w () { - init (a, b, calc_eigenvectors); + init (a, b, calc_rev, calc_lev, force_qz); } FloatEIG (const FloatMatrix& a, const FloatMatrix& b, octave_idx_type& info, - bool calc_eigenvectors = true) - : lambda (), v () + bool calc_rev = true, bool calc_lev = true, bool force_qz = false) + : lambda (), v (), w () { - info = init (a, b, calc_eigenvectors); + info = init (a, b, calc_rev, calc_lev, force_qz); } - FloatEIG (const FloatComplexMatrix& a, bool calc_eigenvectors = true) - : lambda (), v () + FloatEIG (const FloatComplexMatrix& a, bool calc_rev = true, + bool calc_lev = true, bool balance = true) + : lambda (), v (), w () { - init (a, calc_eigenvectors); + init (a, calc_rev, calc_lev, balance); } FloatEIG (const FloatComplexMatrix& a, octave_idx_type& info, - bool calc_eigenvectors = true) - : lambda (), v () + bool calc_rev = true, bool calc_lev = true, bool balance = true) + : lambda (), v (), w () { - info = init (a, calc_eigenvectors); + info = init (a, calc_rev, calc_lev, balance); } FloatEIG (const FloatComplexMatrix& a, const FloatComplexMatrix& b, - bool calc_eigenvectors = true) - : lambda (), v () + bool calc_rev = true, bool calc_lev = true, bool force_qz = false) + : lambda (), v (), w () { - init (a, b, calc_eigenvectors); + init (a, b, calc_rev, calc_lev, force_qz); } FloatEIG (const FloatComplexMatrix& a, const FloatComplexMatrix& b, - octave_idx_type& info, bool calc_eigenvectors = true) - : lambda (), v () + octave_idx_type& info, bool calc_rev = true, bool calc_lev = true, + bool force_qz = false) + : lambda (), v (), w () { - info = init (a, b, calc_eigenvectors); + info = init (a, b, calc_rev, calc_lev, force_qz); } - FloatEIG (const FloatEIG& a) : lambda (a.lambda), v (a.v) { } + FloatEIG (const FloatEIG& a) : lambda (a.lambda), v (a.v), w (a.w) { } FloatEIG& operator = (const FloatEIG& a) { @@ -105,6 +107,7 @@ { lambda = a.lambda; v = a.v; + w = a.w; } return *this; } @@ -112,8 +115,8 @@ ~FloatEIG (void) { } FloatComplexColumnVector eigenvalues (void) const { return lambda; } - - FloatComplexMatrix eigenvectors (void) const { return v; } + FloatComplexMatrix right_eigenvectors (void) const { return v; } + FloatComplexMatrix left_eigenvectors (void) const { return w; } friend std::ostream& operator << (std::ostream& os, const FloatEIG& a); @@ -121,22 +124,34 @@ FloatComplexColumnVector lambda; FloatComplexMatrix v; + FloatComplexMatrix w; - octave_idx_type init (const FloatMatrix& a, bool calc_eigenvectors); + octave_idx_type init (const FloatMatrix& a, bool calc_rev, bool calc_lev, + bool balance); + octave_idx_type init (const FloatMatrix& a, const FloatMatrix& b, - bool calc_eigenvectors); - octave_idx_type init (const FloatComplexMatrix& a, bool calc_eigenvectors); + bool calc_rev, bool calc_lev, bool force_qz); + + octave_idx_type init (const FloatComplexMatrix& a, bool calc_rev, + bool calc_lev, bool balance); + octave_idx_type init (const FloatComplexMatrix& a, - const FloatComplexMatrix& b, bool calc_eigenvectors); + const FloatComplexMatrix& b, + bool calc_rev, bool calc_lev, bool force_qz); - octave_idx_type symmetric_init (const FloatMatrix& a, bool calc_eigenvectors); + octave_idx_type symmetric_init (const FloatMatrix& a, bool calc_rev, + bool calc_lev); + octave_idx_type symmetric_init (const FloatMatrix& a, const FloatMatrix& b, - bool calc_eigenvectors); + bool calc_rev, bool calc_lev); + octave_idx_type hermitian_init (const FloatComplexMatrix& a, - bool calc_eigenvectors); + bool calc_rev, bool calc_lev); + octave_idx_type hermitian_init (const FloatComplexMatrix& a, const FloatComplexMatrix& b, - bool calc_eigenvectors); + bool calc_rev, bool calc_lev); + }; #endif