Mercurial > octave
changeset 25696:3b6691ff0f59
Fix eigs for generalized problems with user defined function (bug #54167).
* eigs-base.cc (EigsRealSymmetricFunc, EigsRealNonSymmetricFunc,
EigsComplexNonSymmetricFunc): Add B as input argument and manage mode = 1 and
mode = 3.
* eigs-base.h: Declare EigsRealSymmetricFunc, EigsRealNonSymmetricFunc, and
EigsComplexNonSymmetricFunc template functions.
* __eigs__.cc: Call EigsRealSymmetricFunc, EigsRealNonSymmetricFunc, and
EigsComplexNonSymmetricFunc with the new input argument B.
* scripts/sparse/eigs.m: Add BIST tests.
| author | Marco Caliari <marco.caliari@univr.it> |
|---|---|
| date | Wed, 27 Jun 2018 08:36:55 +0200 |
| parents | 038fb01854a0 |
| children | 91f0416c2ba7 |
| files | libinterp/dldfcn/__eigs__.cc liboctave/numeric/eigs-base.cc liboctave/numeric/eigs-base.h scripts/sparse/eigs.m |
| diffstat | 4 files changed, 630 insertions(+), 87 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/dldfcn/__eigs__.cc Thu Jul 19 22:43:17 2018 +0200 +++ b/libinterp/dldfcn/__eigs__.cc Wed Jun 27 08:36:55 2018 +0200 @@ -182,6 +182,7 @@ bool sym_tested = false; bool cholB = false; bool a_is_sparse = false; + bool b_is_sparse = false; ColumnVector permB; int arg_offset = 0; double tol = std::numeric_limits<double>::epsilon (); @@ -265,6 +266,13 @@ if (args(1+arg_offset).iscomplex ()) { b_arg = 1+arg_offset; + if (args(b_arg).issparse ()) + { + bscm = (args(b_arg).sparse_complex_matrix_value ()); + b_is_sparse = true; + } + else + bcm = (args(b_arg).complex_matrix_value ()); have_b = true; b_is_complex = true; arg_offset++; @@ -272,6 +280,13 @@ else { b_arg = 1+arg_offset; + if (args(b_arg).issparse ()) + { + bsmm = (args(b_arg).sparse_matrix_value ()); + b_is_sparse = true; + } + else + bmm = (args(b_arg).matrix_value ()); have_b = true; arg_offset++; } @@ -374,14 +389,14 @@ { if (a_is_complex || b_is_complex) { - if (a_is_sparse) + if (b_is_sparse) bscm = args(b_arg).sparse_complex_matrix_value (); else bcm = args(b_arg).complex_matrix_value (); } else { - if (a_is_sparse) + if (b_is_sparse) bsmm = args(b_arg).sparse_matrix_value (); else bmm = args(b_arg).matrix_value (); @@ -400,10 +415,18 @@ ComplexColumnVector eig_val; if (have_a_fun) - nconv = EigsComplexNonSymmetricFunc - (eigs_complex_func, n, typ, sigma, k, p, info, eig_vec, - eig_val, cresid, octave_stdout, tol, (nargout > 1), cholB, - disp, maxit); + { + if (b_is_sparse) + nconv = EigsComplexNonSymmetricFunc + (eigs_complex_func, n, typ, sigma, k, p, info, eig_vec, + eig_val, bscm, permB, cresid, octave_stdout, tol, + (nargout > 1), cholB, disp, maxit); + else + nconv = EigsComplexNonSymmetricFunc + (eigs_complex_func, n, typ, sigma, k, p, info, eig_vec, + eig_val, bcm, permB, cresid, octave_stdout, tol, + (nargout > 1), cholB, disp, maxit); + } else if (have_sigma) { if (a_is_sparse) @@ -443,10 +466,18 @@ ComplexColumnVector eig_val; if (have_a_fun) - nconv = EigsComplexNonSymmetricFunc - (eigs_complex_func, n, typ, sigma, k, p, info, eig_vec, - eig_val, cresid, octave_stdout, tol, (nargout > 1), cholB, - disp, maxit); + { + if (b_is_sparse) + nconv = EigsComplexNonSymmetricFunc + (eigs_complex_func, n, typ, sigma, k, p, info, eig_vec, + eig_val, bscm, permB, cresid, octave_stdout, tol, + (nargout > 1), cholB, disp, maxit); + else + nconv = EigsComplexNonSymmetricFunc + (eigs_complex_func, n, typ, sigma, k, p, info, eig_vec, + eig_val, bcm, permB, cresid, octave_stdout, tol, + (nargout > 1), cholB, disp, maxit); + } else { if (a_is_sparse) @@ -474,10 +505,18 @@ ColumnVector eig_val; if (have_a_fun) - nconv = EigsRealSymmetricFunc - (eigs_func, n, typ, sigmar, k, p, info, eig_vec, - eig_val, resid, octave_stdout, tol, (nargout > 1), - cholB, disp, maxit); + { + if (b_is_sparse) + nconv = EigsRealSymmetricFunc + (eigs_func, n, typ, sigmar, k, p, info, eig_vec, + eig_val, bsmm, permB, resid, octave_stdout, tol, + (nargout > 1), cholB, disp, maxit); + else + nconv = EigsRealSymmetricFunc + (eigs_func, n, typ, sigmar, k, p, info, eig_vec, + eig_val, bmm, permB, resid, octave_stdout, tol, + (nargout > 1), cholB, disp, maxit); + } else if (have_sigma) { if (a_is_sparse) @@ -516,10 +555,18 @@ ComplexColumnVector eig_val; if (have_a_fun) - nconv = EigsRealNonSymmetricFunc - (eigs_func, n, typ, sigmar, k, p, info, eig_vec, - eig_val, resid, octave_stdout, tol, (nargout > 1), - cholB, disp, maxit); + { + if (b_is_sparse) + nconv = EigsRealNonSymmetricFunc + (eigs_func, n, typ, sigmar, k, p, info, eig_vec, + eig_val, bsmm, permB, resid, octave_stdout, tol, + (nargout > 1), cholB, disp, maxit); + else + nconv = EigsRealNonSymmetricFunc + (eigs_func, n, typ, sigmar, k, p, info, eig_vec, + eig_val, bmm, permB, resid, octave_stdout, tol, + (nargout > 1), cholB, disp, maxit); + } else if (have_sigma) { if (a_is_sparse)
--- a/liboctave/numeric/eigs-base.cc Thu Jul 19 22:43:17 2018 +0200 +++ b/liboctave/numeric/eigs-base.cc Wed Jun 27 08:36:55 2018 +0200 @@ -1202,23 +1202,29 @@ return ip(4); } +template <typename M> octave_idx_type EigsRealSymmetricFunc (EigsFunc fun, octave_idx_type n_arg, const std::string& _typ, double sigma, octave_idx_type k_arg, octave_idx_type p_arg, octave_idx_type& info, Matrix& eig_vec, - ColumnVector& eig_val, ColumnVector& resid, + ColumnVector& eig_val, const M& _b, + ColumnVector& permB, ColumnVector& resid, std::ostream& os, double tol, bool rvec, - bool /* cholB */, int disp, int maxit) + bool cholB, int disp, int maxit) { F77_INT n = octave::to_f77_int (n_arg); F77_INT k = octave::to_f77_int (k_arg); F77_INT p = octave::to_f77_int (p_arg); + M b(_b); std::string typ (_typ); bool have_sigma = (sigma ? true : false); + bool have_b = ! b.isempty (); + bool note3 = false; char bmat = 'I'; F77_INT mode = 1; int err = 0; + M bt; if (resid.isempty ()) { @@ -1254,6 +1260,23 @@ (*current_liboctave_error_handler) ("eigs: opts.p must be greater than k and less than or equal to n"); + if (have_b && cholB && ! permB.isempty ()) + { + // Check the we really have a permutation vector + if (permB.numel () != n) + (*current_liboctave_error_handler) ("eigs: permB vector invalid"); + + Array<bool> checked (dim_vector (n, 1), false); + for (F77_INT i = 0; i < n; i++) + { + octave_idx_type bidx = static_cast<octave_idx_type> (permB(i)); + + if (checked(bidx) || bidx < 0 || bidx >= n + || octave::math::x_nint (bidx) != bidx) + (*current_liboctave_error_handler) ("eigs: permB vector invalid"); + } + } + if (! have_sigma) { if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" @@ -1265,19 +1288,53 @@ (*current_liboctave_error_handler) ("eigs: invalid sigma value for real symmetric problem"); + if (typ != "SM" && have_b) + note3 = true; + if (typ == "SM") { typ = "LM"; sigma = 0.; mode = 3; + if (have_b) + bmat = 'G'; } } else if (! std::abs (sigma)) - typ = "SM"; + { + typ = "SM"; + if (have_b) + bmat = 'G'; + } else { typ = "LM"; mode = 3; + if (have_b) + bmat = 'G'; + } + + if (mode == 1 && have_b) + { + // See Note 3 dsaupd + note3 = true; + if (cholB) + { + bt = b; + b = b.transpose (); + if (permB.isempty ()) + { + permB = ColumnVector (n); + for (F77_INT i = 0; i < n; i++) + permB(i) = i; + } + } + else + { + if (! make_cholb (b, bt, permB)) + (*current_liboctave_error_handler) + ("eigs: The matrix B is not positive definite"); + } } Array<F77_INT> ip (dim_vector (11, 1)); @@ -1352,20 +1409,85 @@ if (ido == -1 || ido == 1 || ido == 2) { - double *ip2 = workd + iptr(0) - 1; - ColumnVector x(n); - - for (F77_INT i = 0; i < n; i++) - x(i) = *ip2++; - - ColumnVector y = fun (x, err); - - if (err) - return false; - - ip2 = workd + iptr(1) - 1; - for (F77_INT i = 0; i < n; i++) - *ip2++ = y(i); + if (have_b) + { + if (mode == 1) // regular mode with factorized B + { + Matrix mtmp (n,1); + for (F77_INT i = 0; i < n; i++) + mtmp(i,0) = workd[i + iptr(0) - 1]; + + mtmp = utsolve (bt, permB, mtmp); + ColumnVector y = fun (mtmp, err); + + if (err) + return false; + + mtmp = ltsolve (b, permB, y); + + for (F77_INT i = 0; i < n; i++) + workd[i+iptr(1)-1] = mtmp(i,0); + } + else // shift-invert mode + { + if (ido == -1) + { + OCTAVE_LOCAL_BUFFER (double, dtmp, n); + + vector_product (b, workd+iptr(0)-1, dtmp); + + ColumnVector x(n); + + for (F77_INT i = 0; i < n; i++) + x(i) = dtmp[i]; + + ColumnVector y = fun (x, err); + + if (err) + return false; + + double *ip2 = workd + iptr(1) - 1; + for (F77_INT i = 0; i < n; i++) + ip2[i] = y(i); + } + else if (ido == 2) + vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1); + else + { + double *ip2 = workd+iptr(2)-1; + ColumnVector x(n); + + for (F77_INT i = 0; i < n; i++) + x(i) = *ip2++; + + ColumnVector y = fun (x, err); + + if (err) + return false; + + ip2 = workd + iptr(1) - 1; + for (F77_INT i = 0; i < n; i++) + *ip2++ = y(i); + } + } + } + else + { + double *ip2 = workd + iptr(0) - 1; + ColumnVector x(n); + + for (F77_INT i = 0; i < n; i++) + x(i) = *ip2++; + + ColumnVector y = fun (x, err); + + if (err) + return false; + + ip2 = workd + iptr(1) - 1; + for (F77_INT i = 0; i < n; i++) + *ip2++ = y(i); + } } else { @@ -1408,7 +1530,7 @@ for (F77_INT i = ip(4); i < k; i++) d[i] = octave::numeric_limits<double>::NaN (); F77_INT k2 = ip(4) / 2; - if (typ != "SM" && typ != "BE") + if (mode == 3 || (mode == 1 && typ != "SM" && typ != "BE")) { for (F77_INT i = 0; i < k2; i++) { @@ -1426,7 +1548,7 @@ for (F77_INT j = 0; j < n; j++) z[off1 + j] = octave::numeric_limits<double>::NaN (); } - if (typ != "SM" && typ != "BE") + if (mode == 3 || (mode == 1 && typ != "SM" && typ != "BE")) { OCTAVE_LOCAL_BUFFER (double, dtmp, n); @@ -1448,7 +1570,9 @@ z[off2 + j] = dtmp[j]; } } - } + if (note3) + eig_vec = utsolve (bt, permB, eig_vec); + } } else (*current_liboctave_error_handler) ("eigs: error %d in dseupd", info2); @@ -2179,24 +2303,30 @@ return ip(4); } +template <typename M> octave_idx_type EigsRealNonSymmetricFunc (EigsFunc fun, octave_idx_type n_arg, const std::string& _typ, double sigmar, octave_idx_type k_arg, octave_idx_type p_arg, octave_idx_type& info, ComplexMatrix& eig_vec, - ComplexColumnVector& eig_val, ColumnVector& resid, + ComplexColumnVector& eig_val, const M& _b, + ColumnVector& permB, ColumnVector& resid, std::ostream& os, double tol, bool rvec, - bool /* cholB */, int disp, int maxit) + bool cholB, int disp, int maxit) { F77_INT n = octave::to_f77_int (n_arg); F77_INT k = octave::to_f77_int (k_arg); F77_INT p = octave::to_f77_int (p_arg); + M b(_b); std::string typ (_typ); bool have_sigma = (sigmar ? true : false); - char bmat = 'I'; double sigmai = 0.; F77_INT mode = 1; + bool have_b = ! b.isempty (); + bool note3 = false; + char bmat = 'I'; int err = 0; + M bt; if (resid.isempty ()) { @@ -2232,6 +2362,23 @@ (*current_liboctave_error_handler) ("eigs: opts.p must be greater than k and less than or equal to n"); + if (have_b && cholB && ! permB.isempty ()) + { + // Check the we really have a permutation vector + if (permB.numel () != n) + (*current_liboctave_error_handler) ("eigs: permB vector invalid"); + + Array<bool> checked (dim_vector (n, 1), false); + for (F77_INT i = 0; i < n; i++) + { + octave_idx_type bidx = static_cast<octave_idx_type> (permB(i)); + + if (checked(bidx) || bidx < 0 || bidx >= n + || octave::math::x_nint (bidx) != bidx) + (*current_liboctave_error_handler) ("eigs: permB vector invalid"); + } + } + if (! have_sigma) { if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" @@ -2243,19 +2390,53 @@ (*current_liboctave_error_handler) ("eigs: invalid sigma value for unsymmetric problem"); + if (typ != "SM" && have_b) + note3 = true; + if (typ == "SM") { typ = "LM"; sigmar = 0.; mode = 3; + if (have_b) + bmat = 'G'; } + } + else if (! std::abs (sigmar)) + { + typ = "SM"; + if (have_b) + bmat = 'G'; } - else if (! std::abs (sigmar)) - typ = "SM"; else { typ = "LM"; mode = 3; + if (have_b) + bmat = 'G'; + } + + if (mode == 1 && have_b) + { + // See Note 3 dsaupd + note3 = true; + if (cholB) + { + bt = b; + b = b.transpose (); + if (permB.isempty ()) + { + permB = ColumnVector (n); + for (F77_INT i = 0; i < n; i++) + permB(i) = i; + } + } + else + { + if (! make_cholb (b, bt, permB)) + (*current_liboctave_error_handler) + ("eigs: The matrix B is not positive definite"); + } } Array<F77_INT> ip (dim_vector (11, 1)); @@ -2334,20 +2515,85 @@ if (ido == -1 || ido == 1 || ido == 2) { - double *ip2 = workd + iptr(0) - 1; - ColumnVector x(n); - - for (F77_INT i = 0; i < n; i++) - x(i) = *ip2++; - - ColumnVector y = fun (x, err); - - if (err) - return false; - - ip2 = workd + iptr(1) - 1; - for (F77_INT i = 0; i < n; i++) - *ip2++ = y(i); + if (have_b) + { + if (mode == 1) // regular mode with factorized B + { + Matrix mtmp (n,1); + for (F77_INT i = 0; i < n; i++) + mtmp(i,0) = workd[i + iptr(0) - 1]; + + mtmp = utsolve (bt, permB, mtmp); + ColumnVector y = fun (mtmp, err); + + if (err) + return false; + + mtmp = ltsolve (b, permB, y); + + for (F77_INT i = 0; i < n; i++) + workd[i+iptr(1)-1] = mtmp(i,0); + } + else // shift-invert mode + { + if (ido == -1) + { + OCTAVE_LOCAL_BUFFER (double, dtmp, n); + + vector_product (b, workd+iptr(0)-1, dtmp); + + ColumnVector x(n); + + for (F77_INT i = 0; i < n; i++) + x(i) = dtmp[i]; + + ColumnVector y = fun (x, err); + + if (err) + return false; + + double *ip2 = workd + iptr(1) - 1; + for (F77_INT i = 0; i < n; i++) + ip2[i] = y(i); + } + else if (ido == 2) + vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1); + else + { + double *ip2 = workd+iptr(2)-1; + ColumnVector x(n); + + for (F77_INT i = 0; i < n; i++) + x(i) = *ip2++; + + ColumnVector y = fun (x, err); + + if (err) + return false; + + ip2 = workd + iptr(1) - 1; + for (F77_INT i = 0; i < n; i++) + *ip2++ = y(i); + } + } + } + else + { + double *ip2 = workd + iptr(0) - 1; + ColumnVector x(n); + + for (F77_INT i = 0; i < n; i++) + x(i) = *ip2++; + + ColumnVector y = fun (x, err); + + if (err) + return false; + + ip2 = workd + iptr(1) - 1; + for (F77_INT i = 0; i < n; i++) + *ip2++ = y(i); + } } else { @@ -2485,6 +2731,8 @@ eig_vec(jj,ii) = Complex (octave::numeric_limits<double>::NaN (), 0.); } + if (note3) + eig_vec = utsolve (bt, permB, eig_vec); } if (k0 < k) { @@ -3119,24 +3367,29 @@ return ip(4); } +template <typename M> octave_idx_type EigsComplexNonSymmetricFunc (EigsComplexFunc fun, octave_idx_type n_arg, const std::string& _typ, Complex sigma, octave_idx_type k_arg, octave_idx_type p_arg, octave_idx_type& info, ComplexMatrix& eig_vec, - ComplexColumnVector& eig_val, - ComplexColumnVector& cresid, std::ostream& os, - double tol, bool rvec, bool /* cholB */, - int disp, int maxit) + ComplexColumnVector& eig_val, const M& _b, + ColumnVector& permB, ComplexColumnVector& cresid, + std::ostream& os, double tol, bool rvec, + bool cholB, int disp, int maxit) { F77_INT n = octave::to_f77_int (n_arg); F77_INT k = octave::to_f77_int (k_arg); F77_INT p = octave::to_f77_int (p_arg); + M b(_b); std::string typ (_typ); bool have_sigma = (std::abs (sigma) ? true : false); + F77_INT mode = 1; + bool have_b = ! b.isempty (); + bool note3 = false; char bmat = 'I'; - F77_INT mode = 1; int err = 0; + M bt; if (cresid.isempty ()) { @@ -3176,6 +3429,23 @@ (*current_liboctave_error_handler) ("eigs: opts.p must be greater than k and less than or equal to n"); + if (have_b && cholB && ! permB.isempty ()) + { + // Check the we really have a permutation vector + if (permB.numel () != n) + (*current_liboctave_error_handler) ("eigs: permB vector invalid"); + + Array<bool> checked (dim_vector (n, 1), false); + for (F77_INT i = 0; i < n; i++) + { + octave_idx_type bidx = static_cast<octave_idx_type> (permB(i)); + + if (checked(bidx) || bidx < 0 || bidx >= n + || octave::math::x_nint (bidx) != bidx) + (*current_liboctave_error_handler) ("eigs: permB vector invalid"); + } + } + if (! have_sigma) { if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" @@ -3187,19 +3457,53 @@ (*current_liboctave_error_handler) ("eigs: invalid sigma value for complex problem"); + if (typ != "SM" && have_b) + note3 = true; + if (typ == "SM") { typ = "LM"; sigma = 0.; mode = 3; + if (have_b) + bmat ='G'; } } else if (! std::abs (sigma)) - typ = "SM"; + { + typ = "SM"; + if (have_b) + bmat = 'G'; + } else { typ = "LM"; mode = 3; + if (have_b) + bmat = 'G'; + } + + if (mode == 1 && have_b) + { + // See Note 3 dsaupd + note3 = true; + if (cholB) + { + bt = b; + b = b.hermitian (); + if (permB.isempty ()) + { + permB = ColumnVector (n); + for (F77_INT i = 0; i < n; i++) + permB(i) = i; + } + } + else + { + if (! make_cholb (b, bt, permB)) + (*current_liboctave_error_handler) + ("eigs: The matrix B is not positive definite"); + } } Array<F77_INT> ip (dim_vector (11, 1)); @@ -3276,20 +3580,85 @@ if (ido == -1 || ido == 1 || ido == 2) { - Complex *ip2 = workd + iptr(0) - 1; - ComplexColumnVector x(n); - - for (F77_INT i = 0; i < n; i++) - x(i) = *ip2++; - - ComplexColumnVector y = fun (x, err); - - if (err) - return false; - - ip2 = workd + iptr(1) - 1; - for (F77_INT i = 0; i < n; i++) - *ip2++ = y(i); + if (have_b) + { + if (mode == 1) // regular mode with factorized B + { + ComplexMatrix mtmp (n,1); + for (F77_INT i = 0; i < n; i++) + mtmp(i,0) = workd[i + iptr(0) - 1]; + + mtmp = utsolve (bt, permB, mtmp); + ComplexColumnVector y = fun (mtmp, err); + + if (err) + return false; + + mtmp = ltsolve (b, permB, y); + + for (F77_INT i = 0; i < n; i++) + workd[i+iptr(1)-1] = mtmp(i,0); + } + else // shift-invert mode + { + if (ido == -1) + { + OCTAVE_LOCAL_BUFFER (Complex, ctmp, n); + + vector_product (b, workd+iptr(0)-1, ctmp); + + ComplexColumnVector x(n); + + for (F77_INT i = 0; i < n; i++) + x(i) = ctmp[i]; + + ComplexColumnVector y = fun (x, err); + + if (err) + return false; + + Complex *ip2 = workd+iptr(1)-1; + for (F77_INT i = 0; i < n; i++) + ip2[i] = y(i); + } + else if (ido == 2) + vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1); + else + { + Complex *ip2 = workd+iptr(2)-1; + ComplexColumnVector x(n); + + for (F77_INT i = 0; i < n; i++) + x(i) = *ip2++; + + ComplexColumnVector y = fun (x, err); + + if (err) + return false; + + ip2 = workd + iptr(1) - 1; + for (F77_INT i = 0; i < n; i++) + *ip2++ = y(i); + } + } + } + else + { + Complex *ip2 = workd + iptr(0) - 1; + ComplexColumnVector x(n); + + for (F77_INT i = 0; i < n; i++) + x(i) = *ip2++; + + ComplexColumnVector y = fun (x, err); + + if (err) + return false; + + ip2 = workd + iptr(1) - 1; + for (F77_INT i = 0; i < n; i++) + *ip2++ = y(i); + } } else { @@ -3377,6 +3746,8 @@ for (F77_INT j = 0; j < n; j++) z[off2 + j] = ctmp[j]; } + if (note3) + eig_vec = utsolve (bt, permB, eig_vec); } } else @@ -3409,6 +3780,15 @@ template octave_idx_type +EigsRealSymmetricFunc<Matrix> +(EigsFunc fun, octave_idx_type n, const std::string& _typ, double sigma, + octave_idx_type k, octave_idx_type p, octave_idx_type& info, + Matrix& eig_vec, ColumnVector& eig_val, const Matrix& _b, + ColumnVector& permB, ColumnVector& resid, std::ostream& os, double tol, + bool rvec, bool cholB, int disp, int maxit); + +template +octave_idx_type EigsRealNonSymmetricMatrix<Matrix> (const Matrix& m, const std::string typ, octave_idx_type k, octave_idx_type p, octave_idx_type& info, ComplexMatrix& eig_vec, @@ -3425,6 +3805,15 @@ ColumnVector& resid, std::ostream& os, double tol, bool rvec, bool cholB, int disp, int maxit); +template +octave_idx_type +EigsRealNonSymmetricFunc<Matrix> +(EigsFunc fun, octave_idx_type n, const std::string& _typ, double sigmar, + octave_idx_type k, octave_idx_type p, octave_idx_type& info, + ComplexMatrix& eig_vec, ComplexColumnVector& eig_val, const Matrix& _b, + ColumnVector& permB, ColumnVector& resid, std::ostream& os, double tol, + bool rvec, bool cholB, int disp, int maxit); + // SparseMatrix template @@ -3447,6 +3836,15 @@ template octave_idx_type +EigsRealSymmetricFunc<SparseMatrix> +(EigsFunc fun, octave_idx_type n, const std::string& _typ, double sigma, + octave_idx_type k, octave_idx_type p, octave_idx_type& info, + Matrix& eig_vec, ColumnVector& eig_val, const SparseMatrix& _b, + ColumnVector& permB, ColumnVector& resid, std::ostream& os, double tol, + bool rvec, bool cholB, int disp, int maxit); + +template +octave_idx_type EigsRealNonSymmetricMatrix<SparseMatrix> (const SparseMatrix& m, const std::string typ, octave_idx_type k, octave_idx_type p, octave_idx_type& info, ComplexMatrix& eig_vec, @@ -3463,6 +3861,15 @@ ColumnVector& resid, std::ostream& os, double tol, bool rvec, bool cholB, int disp, int maxit); +template +octave_idx_type +EigsRealNonSymmetricFunc<SparseMatrix> +(EigsFunc fun, octave_idx_type n, const std::string& _typ, double sigmar, + octave_idx_type k, octave_idx_type p, octave_idx_type& info, + ComplexMatrix& eig_vec, ComplexColumnVector& eig_val, + const SparseMatrix& _b, ColumnVector& permB, ColumnVector& resid, + std::ostream& os, double tol, bool rvec, bool cholB, int disp, int maxit); + // ComplexMatrix template @@ -3483,6 +3890,15 @@ ComplexColumnVector& cresid, std::ostream& os, double tol, bool rvec, bool cholB, int disp, int maxit); +template +octave_idx_type +EigsComplexNonSymmetricFunc<ComplexMatrix> +(EigsComplexFunc fun, octave_idx_type n, const std::string& _typ, Complex sigma, + octave_idx_type k, octave_idx_type p, octave_idx_type& info, + ComplexMatrix& eig_vec, ComplexColumnVector& eig_val, + const ComplexMatrix& _b, ColumnVector& permB, ComplexColumnVector& cresid, + std::ostream& os, double tol, bool rvec, bool cholB, int disp, int maxit); + // SparseComplexMatrix template @@ -3503,4 +3919,13 @@ ColumnVector& permB, ComplexColumnVector& cresid, std::ostream& os, double tol, bool rvec, bool cholB, int disp, int maxit); +template +octave_idx_type +EigsComplexNonSymmetricFunc<SparseComplexMatrix> +(EigsComplexFunc fun, octave_idx_type n, const std::string& _typ, Complex sigma, + octave_idx_type k, octave_idx_type p, octave_idx_type& info, + ComplexMatrix& eig_vec, ComplexColumnVector& eig_val, + const SparseComplexMatrix& _b, ColumnVector& permB, + ComplexColumnVector& cresid, std::ostream& os, double tol, bool rvec, + bool cholB, int disp, int maxit); #endif
--- a/liboctave/numeric/eigs-base.h Thu Jul 19 22:43:17 2018 +0200 +++ b/liboctave/numeric/eigs-base.h Wed Jun 27 08:36:55 2018 +0200 @@ -60,14 +60,16 @@ std::ostream& os, double tol, bool rvec, bool cholB, int disp, int maxit); +template <typename M> extern OCTAVE_API octave_idx_type EigsRealSymmetricFunc (EigsFunc fun, octave_idx_type n, const std::string& _typ, double sigma, octave_idx_type k, octave_idx_type p, octave_idx_type& info, Matrix& eig_vec, - ColumnVector& eig_val, ColumnVector& resid, + ColumnVector& eig_val, const M& _b, + ColumnVector& permB, ColumnVector& resid, std::ostream& os, double tol, bool rvec, - bool /* cholB */, int disp, int maxit); + bool cholB, int disp, int maxit); template <typename M> octave_idx_type @@ -90,14 +92,16 @@ std::ostream& os, double tol, bool rvec, bool cholB, int disp, int maxit); +template <typename M> extern OCTAVE_API octave_idx_type EigsRealNonSymmetricFunc (EigsFunc fun, octave_idx_type n, const std::string& _typ, double sigmar, octave_idx_type k, octave_idx_type p, octave_idx_type& info, ComplexMatrix& eig_vec, - ComplexColumnVector& eig_val, ColumnVector& resid, + ComplexColumnVector& eig_val, const M& _b, + ColumnVector& permB, ColumnVector& resid, std::ostream& os, double tol, bool rvec, - bool /* cholB */, int disp, int maxit); + bool cholB, int disp, int maxit); template <typename M> octave_idx_type @@ -122,14 +126,15 @@ std::ostream& os, double tol, bool rvec, bool cholB, int disp, int maxit); +template <typename M> extern OCTAVE_API octave_idx_type EigsComplexNonSymmetricFunc (EigsComplexFunc fun, octave_idx_type n, const std::string& _typ, Complex sigma, octave_idx_type k, octave_idx_type p, octave_idx_type& info, ComplexMatrix& eig_vec, - ComplexColumnVector& eig_val, - ComplexColumnVector& cresid, std::ostream& os, - double tol, bool rvec, bool /* cholB */, - int disp, int maxit); + ComplexColumnVector& eig_val, const M& _b, + ColumnVector& permB, ComplexColumnVector& cresid, + std::ostream& os, double tol, bool rvec, + bool cholB, int disp, int maxit); #endif
--- a/scripts/sparse/eigs.m Thu Jul 19 22:43:17 2018 +0200 +++ b/scripts/sparse/eigs.m Wed Jun 27 08:36:55 2018 +0200 @@ -1480,11 +1480,77 @@ %! j_B = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]; %! v_B = [3, 10i, 1, 8i, 7, 6i, 5, 4i, 9, 7i]; %! B = sparse(i_B, j_B, v_B); # not SPD -%! [Evectors Evalues] = eigs(A, B, 5, 'SM'); # call_eig is true +%! [Evectors, Evalues] = eigs(A, B, 5, "SM"); # call_eig is true %! ResidualVectors = A * Evectors - B * Evectors * Evalues; %! RelativeErrors = norm (ResidualVectors, "columns") ./ ... %! norm (A * Evectors, "columns"); -%! assert (RelativeErrors, zeros (1, 5)) +%! assert (RelativeErrors, zeros (1, 5)); %!testif HAVE_ARPACK %! A = rand (8); %! eigs (A, 6, "lr"); # this failed on 4.2.x +%!testif HAVE_ARPACK +%! A = rand (10); +%! B = rand (10); +%! B = B * B'; +%! opts.v0 = (1:10)'; +%! [Evector, Evalues] = eigs (A, B, 4, "LM", opts); +%! Afun = @(x) A * x; +%! [Evector_f Evalues_f] = eigs (Afun, 10, B, 4, "LM", opts); +%! assert (Evector, Evector_f); +%! assert (Evalues, Evalues_f); +%!testif HAVE_ARPACK +%! A = rand (10); +%! B = rand (10); +%! B = B * B'; +%! opts.v0 = (1:10)'; +%! [Evector, Evalues] = eigs (A, B, 4, "SM", opts); +%! Afun = @(x) A \ x; +%! [Evector_f Evalues_f] = eigs (Afun, 10, B, 4, "SM", opts); +%! assert (Evector, Evector_f); +%! assert (Evalues, Evalues_f); +%!testif HAVE_ARPACK +%! A = rand (10); +%! A = A * A'; +%! B = rand (10); +%! B = B * B'; +%! opts.v0 = (1:10)'; +%! [Evector, Evalues] = eigs (A, B, 4, "LM", opts); +%! Afun = @(x) A * x; +%! opts.issym = true; +%! [Evector_f Evalues_f] = eigs (Afun, 10, B, 4, "LM", opts); +%! assert (Evector, Evector_f); +%! assert (Evalues, Evalues_f); +%!testif HAVE_ARPACK +%! A = rand (10); +%! A = A * A'; +%! B = rand (10); +%! B = B * B'; +%! opts.v0 = (1:10)'; +%! [Evector, Evalues] = eigs (A, B, 4, "SM", opts); +%! Afun = @(x) A \ x; +%! opts.issym = true; +%! [Evector_f Evalues_f] = eigs (Afun, 10, B, 4, "SM", opts); +%! assert (Evector, Evector_f, 100*eps); +%! assert (Evalues, Evalues_f, 100*eps); +%!testif HAVE_ARPACK +%! A = rand (10) + 1i * rand (10); +%! B = rand (10) + 1i * rand (10); +%! B = B * B'; +%! opts.v0 = (1:10)'; +%! [Evector, Evalues] = eigs (A, B, 4, "LM", opts); +%! Afun = @(x) A * x; +%! opts.isreal = false; +%! [Evector_f Evalues_f] = eigs (Afun, 10, B, 4, "LM", opts); +%! assert (Evector, Evector_f); +%! assert (Evalues, Evalues_f); +%!testif HAVE_ARPACK +%! A = rand (10) + 1i * rand (10); +%! B = rand (10) + 1i * rand (10); +%! B = B * B'; +%! opts.v0 = (1:10)'; +%! [Evector, Evalues] = eigs (A, B, 4, "SM", opts); +%! Afun = @(x) A \ x; +%! opts.isreal = false; +%! [Evector_f, Evalues_f] = eigs (Afun, 10, B, 4, "SM", opts); +%! assert (Evector, Evector_f); +%! assert (Evalues, Evalues_f);