Mercurial > octave
changeset 24228:1310d8b50ec2
Fix eigs returning an incorrect number or order of eigenvalues (bug #45153, bug #47236).
* liboctave/numberic/eigs-base.cc: Always return the required number of
eigenvalues. Remove ido == 2 for non-generalized problems. Fix the maximum
value for p. Fix the order of eigenvalues. Check the dimension of opts.v0.
* script/sparse/eigs.m: Fix the range of p in the description.
| author | Marco Caliari <marco.caliari@univr.it> |
|---|---|
| date | Fri, 24 Mar 2017 14:41:08 +0100 |
| parents | 0da1bdfbfacb |
| children | 97e628756971 |
| files | liboctave/numeric/eigs-base.cc scripts/sparse/eigs.m |
| diffstat | 2 files changed, 272 insertions(+), 220 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/numeric/eigs-base.cc Sat Nov 11 19:39:23 2017 +0100 +++ b/liboctave/numeric/eigs-base.cc Fri Mar 24 14:41:08 2017 +0100 @@ -626,6 +626,10 @@ resid = ColumnVector (octave_rand::vector (n)); octave_rand::distribution (rand_dist); } + else if (m.cols () != resid.numel ()) + { + (*current_liboctave_error_handler) ("eigs: opts.v0 must be n-by-1"); + } if (n < 3) (*current_liboctave_error_handler) ("eigs: n must be at least 3"); @@ -637,8 +641,8 @@ if (p < 20) p = 20; - if (p > n - 1) - p = n - 1; + if (p > n) + p = n; } if (k < 1 || k > n - 2) @@ -647,9 +651,9 @@ " (must be 0 < k < n-1-1).\n" " Use 'eig (full (A))' instead"); - if (p <= k || p >= n) + if (p <= k || p > n) (*current_liboctave_error_handler) - ("eigs: opts.p must be greater than k and less than n"); + ("eigs: opts.p must be greater than k and less or equal to n"); if (have_b && cholB && ! permB.isempty ()) { @@ -825,7 +829,7 @@ if (info2 == 0) { F77_INT k2 = k / 2; - if (typ != "SM" && typ != "BE") + if (typ != "SM" && typ != "BE" && ! (typ == "SA" && rvec)) { for (F77_INT i = 0; i < k2; i++) { @@ -837,7 +841,7 @@ if (rvec) { - if (typ != "SM" && typ != "BE") + if (typ != "SM" && typ != "BE" && typ != "SA") { OCTAVE_LOCAL_BUFFER (double, dtmp, n); @@ -907,6 +911,10 @@ resid = ColumnVector (octave_rand::vector (n)); octave_rand::distribution (rand_dist); } + else if (m.cols () != resid.numel ()) + { + (*current_liboctave_error_handler) ("eigs: opts.v0 must be n-by-1"); + } if (n < 3) (*current_liboctave_error_handler) ("eigs: n must be at least 3"); @@ -924,13 +932,13 @@ if (p < 20) p = 20; - if (p > n - 1) - p = n - 1; + if (p > n) + p = n; } - if (p <= k || p >= n) + if (p <= k || p > n) (*current_liboctave_error_handler) - ("eigs: opts.p must be greater than k and less than n"); + ("eigs: opts.p must be greater than k and less or equal to n"); if (have_b && cholB && ! permB.isempty ()) { @@ -1066,25 +1074,19 @@ } else { - if (ido == 2) - { - for (F77_INT i = 0; i < n; i++) - workd[iptr(0) + i - 1] = workd[iptr(1) + i - 1]; - } - else - { - double *ip2 = workd+iptr(0)-1; - Matrix tmp(n, 1); - - for (F77_INT i = 0; i < n; i++) - tmp(i,0) = ip2[P[i]]; - - lusolve (L, U, tmp); - - ip2 = workd+iptr(1)-1; - for (F77_INT i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } + // ido cannot be 2 for non-generalized problems + // see dsaupd2 + double *ip2 = workd+iptr(0)-1; + Matrix tmp(n, 1); + + for (F77_INT i = 0; i < n; i++) + tmp(i,0) = ip2[P[i]]; + + lusolve (L, U, tmp); + + ip2 = workd+iptr(1)-1; + for (F77_INT i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); } } else @@ -1191,6 +1193,10 @@ resid = ColumnVector (octave_rand::vector (n)); octave_rand::distribution (rand_dist); } + else if (n != resid.numel ()) + { + (*current_liboctave_error_handler) ("eigs: opts.v0 must be n-by-1"); + } if (n < 3) (*current_liboctave_error_handler) ("eigs: n must be at least 3"); @@ -1202,8 +1208,8 @@ if (p < 20) p = 20; - if (p > n - 1) - p = n - 1; + if (p > n) + p = n; } if (k <= 0 || k >= n - 1) @@ -1212,9 +1218,9 @@ " (must be 0 < k < n-1).\n" " Use 'eig (full (A))' instead"); - if (p <= k || p >= n) + if (p <= k || p > n) (*current_liboctave_error_handler) - ("eigs: opts.p must be greater than k and less than n"); + ("eigs: opts.p must be greater than k and less or equal to n"); if (! have_sigma) { @@ -1444,6 +1450,10 @@ resid = ColumnVector (octave_rand::vector (n)); octave_rand::distribution (rand_dist); } + else if (m.cols () != resid.numel ()) + { + (*current_liboctave_error_handler) ("eigs: opts.v0 must be n-by-1"); + } if (n < 3) (*current_liboctave_error_handler) ("eigs: n must be at least 3"); @@ -1455,8 +1465,8 @@ if (p < 20) p = 20; - if (p > n - 1) - p = n - 1; + if (p > n) + p = n; } if (k <= 0 || k >= n - 1) @@ -1465,9 +1475,9 @@ " (must be 0 < k < n-1).\n" " Use 'eig (full (A))' instead"); - if (p <= k || p >= n) + if (p <= k || p > n) (*current_liboctave_error_handler) - ("eigs: opts.p must be greater than k and less than n"); + ("eigs: opts.p must be greater than k and less or equal to n"); if (have_b && cholB && ! permB.isempty ()) { @@ -1550,12 +1560,15 @@ { F77_INT tmp_info = octave::to_f77_int (info); + // on exit, ip(4) <= k + 1 is the number of converged eigenvalues + // see dnaupd2 F77_FUNC (dnaupd, DNAUPD) (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, tmp_info F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + // k is not changed info = tmp_info; @@ -1639,65 +1652,45 @@ for (F77_INT i = 0; i < k+1; i++) dr[i] = di[i] = 0.; + F77_INT k0 = k; // original number of eigenvalues required F77_FUNC (dneupd, DNEUPD) (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar, sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n, F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + // on exit, if (and only if) rvec == true, k may have been increased by one + // and be equal to ip(4), see dngets if (f77_exception_encountered) (*current_liboctave_error_handler) ("eigs: unrecoverable exception encountered in dneupd"); - eig_val.resize (k+1); + if (! rvec && ip(4) > k) + k = ip(4); + + eig_val.resize (k); Complex *d = eig_val.fortran_vec (); if (info2 == 0) { - F77_INT jj = 0; - for (F77_INT i = 0; i < k+1; i++) + for (F77_INT i = 0; i < k; i++) + d[i] = Complex (dr[i], di[i]); + + if (! rvec) { - if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0) - jj++; - else - d[i-jj] = Complex (dr[i], di[i]); - } - if (jj == 0 && ! rvec) - for (F77_INT i = 0; i < k; i++) - d[i] = d[i+1]; - - F77_INT k2 = k / 2; - for (F77_INT i = 0; i < k2; i++) - { - Complex dtmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = dtmp; - } - eig_val.resize (k); - - if (rvec) - { - OCTAVE_LOCAL_BUFFER (double, dtmp, n); - + // ARPACK seems to give the eigenvalues in reversed order + F77_INT k2 = k / 2; for (F77_INT i = 0; i < k2; i++) { - F77_INT off1 = i * n; - F77_INT off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (F77_INT j = 0; j < n; j++) - dtmp[j] = z[off1 + j]; - - for (F77_INT j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (F77_INT j = 0; j < n; j++) - z[off2 + j] = dtmp[j]; + Complex dtmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = dtmp; } - + } + else + { + // When eigenvectors required, ARPACK seems to give the right order eig_vec.resize (n, k); F77_INT i = 0; while (i < k) @@ -1728,6 +1721,11 @@ if (note3) eig_vec = utsolve (bt, permB, eig_vec); } + if (k0 < k) + { + eig_val.resize (k0); + eig_vec.resize (n, k0); + } } else (*current_liboctave_error_handler) ("eigs: error %d in dneupd", info2); @@ -1774,6 +1772,10 @@ resid = ColumnVector (octave_rand::vector (n)); octave_rand::distribution (rand_dist); } + else if (m.cols () != resid.numel ()) + { + (*current_liboctave_error_handler) ("eigs: opts.v0 must be n-by-1"); + } if (n < 3) (*current_liboctave_error_handler) ("eigs: n must be at least 3"); @@ -1785,8 +1787,8 @@ if (p < 20) p = 20; - if (p > n - 1) - p = n - 1; + if (p > n) + p = n; } if (k <= 0 || k >= n - 1) @@ -1795,9 +1797,9 @@ " (must be 0 < k < n-1).\n" " Use 'eig (full (A))' instead"); - if (p <= k || p >= n) + if (p <= k || p > n) (*current_liboctave_error_handler) - ("eigs: opts.p must be greater than k and less than n"); + ("eigs: opts.p must be greater than k and less or equal to n"); if (have_b && cholB && ! permB.isempty ()) { @@ -1860,12 +1862,15 @@ { F77_INT tmp_info = octave::to_f77_int (info); + // on exit, ip(4) <= k + 1 is the number of converged eigenvalues + // see dnaupd2 F77_FUNC (dnaupd, DNAUPD) (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, tmp_info F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + // k is not changed info = tmp_info; @@ -1933,32 +1938,26 @@ } else { - if (ido == 2) - { - for (F77_INT i = 0; i < n; i++) - workd[iptr(0) + i - 1] = workd[iptr(1) + i - 1]; - } - else - { - double *ip2 = workd+iptr(0)-1; - Matrix tmp(n, 1); - - for (F77_INT i = 0; i < n; i++) - tmp(i,0) = ip2[P[i]]; - - lusolve (L, U, tmp); - - ip2 = workd+iptr(1)-1; - for (F77_INT i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } + // ido cannot be 2 for non-generalized problems + // see dnaupd2 + double *ip2 = workd+iptr(0)-1; + Matrix tmp(n, 1); + + for (F77_INT i = 0; i < n; i++) + tmp(i,0) = ip2[P[i]]; + + lusolve (L, U, tmp); + + ip2 = workd+iptr(1)-1; + for (F77_INT i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); } } else { if (info < 0) (*current_liboctave_error_handler) - ("eigs: error %d in dsaupd", info); + ("eigs: error %d in dnaupd", info); break; } @@ -1993,65 +1992,45 @@ for (F77_INT i = 0; i < k+1; i++) dr[i] = di[i] = 0.; + F77_INT k0 = k; // original number of eigenvalues required F77_FUNC (dneupd, DNEUPD) (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar, sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n, F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + // on exit, if (and only if) rvec == true, k may have been increased by one + // and be equal to ip(4), see dngets if (f77_exception_encountered) (*current_liboctave_error_handler) ("eigs: unrecoverable exception encountered in dneupd"); - eig_val.resize (k+1); + if (! rvec && ip(4) > k) + k = ip(4); + + eig_val.resize (k); Complex *d = eig_val.fortran_vec (); if (info2 == 0) { - F77_INT jj = 0; - for (F77_INT i = 0; i < k+1; i++) + for (F77_INT i = 0; i < k; i++) + d[i] = Complex (dr[i], di[i]); + + if (! rvec) { - if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0) - jj++; - else - d[i-jj] = Complex (dr[i], di[i]); - } - if (jj == 0 && ! rvec) - for (F77_INT i = 0; i < k; i++) - d[i] = d[i+1]; - - F77_INT k2 = k / 2; - for (F77_INT i = 0; i < k2; i++) - { - Complex dtmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = dtmp; - } - eig_val.resize (k); - - if (rvec) - { - OCTAVE_LOCAL_BUFFER (double, dtmp, n); - + // ARPACK seems to give the eigenvalues in reversed order + F77_INT k2 = k / 2; for (F77_INT i = 0; i < k2; i++) { - F77_INT off1 = i * n; - F77_INT off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (F77_INT j = 0; j < n; j++) - dtmp[j] = z[off1 + j]; - - for (F77_INT j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (F77_INT j = 0; j < n; j++) - z[off2 + j] = dtmp[j]; + Complex dtmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = dtmp; } - + } + else + { + // When eigenvectors required, ARPACK seems to give the right order eig_vec.resize (n, k); F77_INT i = 0; while (i < k) @@ -2079,6 +2058,11 @@ } } } + if (k0 < k) + { + eig_val.resize (k0); + eig_vec.resize (n, k0); + } } else (*current_liboctave_error_handler) ("eigs: error %d in dneupd", info2); @@ -2112,6 +2096,10 @@ resid = ColumnVector (octave_rand::vector (n)); octave_rand::distribution (rand_dist); } + else if (n != resid.numel ()) + { + (*current_liboctave_error_handler) ("eigs: opts.v0 must be n-by-1"); + } if (n < 3) (*current_liboctave_error_handler) ("eigs: n must be at least 3"); @@ -2123,8 +2111,8 @@ if (p < 20) p = 20; - if (p > n - 1) - p = n - 1; + if (p > n) + p = n; } if (k <= 0 || k >= n - 1) @@ -2133,9 +2121,9 @@ " (must be 0 < k < n-1).\n" " Use 'eig (full (A))' instead"); - if (p <= k || p >= n) + if (p <= k || p > n) (*current_liboctave_error_handler) - ("eigs: opts.p must be greater than k and less than n"); + ("eigs: opts.p must be greater than k and less or equal to n"); if (! have_sigma) { @@ -2195,12 +2183,15 @@ { F77_INT tmp_info = octave::to_f77_int (info); + // on exit, ip(4) <= k + 1 is the number of converged eigenvalues + // see dnaupd2 F77_FUNC (dnaupd, DNAUPD) - (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, tmp_info F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + // k is not changed info = tmp_info; @@ -2284,65 +2275,45 @@ for (F77_INT i = 0; i < k+1; i++) dr[i] = di[i] = 0.; + F77_INT k0 = k; // original number of eigenvalues required F77_FUNC (dneupd, DNEUPD) (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar, sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n, F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + // on exit, if (and only if) rvec == true, k may have been increased by one + // and be equal to ip(4), see dngets if (f77_exception_encountered) (*current_liboctave_error_handler) ("eigs: unrecoverable exception encountered in dneupd"); - eig_val.resize (k+1); + if (! rvec && ip(4) > k) + k = ip(4); + + eig_val.resize (k); Complex *d = eig_val.fortran_vec (); if (info2 == 0) { - F77_INT jj = 0; - for (F77_INT i = 0; i < k+1; i++) + for (F77_INT i = 0; i < k; i++) + d[i] = Complex (dr[i], di[i]); + + if (! rvec) { - if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0) - jj++; - else - d[i-jj] = Complex (dr[i], di[i]); - } - if (jj == 0 && ! rvec) - for (F77_INT i = 0; i < k; i++) - d[i] = d[i+1]; - - F77_INT k2 = k / 2; - for (F77_INT i = 0; i < k2; i++) - { - Complex dtmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = dtmp; - } - eig_val.resize (k); - - if (rvec) - { - OCTAVE_LOCAL_BUFFER (double, dtmp, n); - + // ARPACK seems to give the eigenvalues in reversed order + octave_idx_type k2 = k / 2; for (F77_INT i = 0; i < k2; i++) { - F77_INT off1 = i * n; - F77_INT off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (F77_INT j = 0; j < n; j++) - dtmp[j] = z[off1 + j]; - - for (F77_INT j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (F77_INT j = 0; j < n; j++) - z[off2 + j] = dtmp[j]; + Complex dtmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = dtmp; } - + } + else + { + // ARPACK seems to give the eigenvalues in reversed order eig_vec.resize (n, k); F77_INT i = 0; while (i < k) @@ -2370,6 +2341,11 @@ } } } + if (k0 < k) + { + eig_val.resize (k0); + eig_vec.resize (n, k0); + } } else (*current_liboctave_error_handler) ("eigs: error %d in dneupd", info2); @@ -2416,6 +2392,10 @@ cresid(i) = Complex (rr(i),ri(i)); octave_rand::distribution (rand_dist); } + else if (m.cols () != cresid.numel ()) + { + (*current_liboctave_error_handler) ("eigs: opts.v0 must be n-by-1"); + } if (n < 3) (*current_liboctave_error_handler) ("eigs: n must be at least 3"); @@ -2427,8 +2407,8 @@ if (p < 20) p = 20; - if (p > n - 1) - p = n - 1; + if (p > n) + p = n; } if (k <= 0 || k >= n - 1) @@ -2437,9 +2417,9 @@ " (must be 0 < k < n-1).\n" " Use 'eig (full (A))' instead"); - if (p <= k || p >= n) + if (p <= k || p > n) (*current_liboctave_error_handler) - ("eigs: opts.p must be greater than k and less than n"); + ("eigs: opts.p must be greater than k and less or equal to n"); if (have_b && cholB && ! permB.isempty ()) { @@ -2705,6 +2685,10 @@ cresid(i) = Complex (rr(i),ri(i)); octave_rand::distribution (rand_dist); } + else if (m.cols () != cresid.numel ()) + { + (*current_liboctave_error_handler) ("eigs: opts.v0 must be n-by-1"); + } if (n < 3) (*current_liboctave_error_handler) ("eigs: n must be at least 3"); @@ -2716,8 +2700,8 @@ if (p < 20) p = 20; - if (p > n - 1) - p = n - 1; + if (p > n) + p = n; } if (k <= 0 || k >= n - 1) @@ -2726,9 +2710,9 @@ " (must be 0 < k < n-1).\n" " Use 'eig (full (A))' instead"); - if (p <= k || p >= n) + if (p <= k || p > n) (*current_liboctave_error_handler) - ("eigs: opts.p must be greater than k and less than n"); + ("eigs: opts.p must be greater than k and less or equal to n"); if (have_b && cholB && ! permB.isempty ()) { @@ -2866,33 +2850,26 @@ } else { - if (ido == 2) - { - for (F77_INT i = 0; i < n; i++) - workd[iptr(0) + i - 1] = - workd[iptr(1) + i - 1]; - } - else - { - Complex *ip2 = workd+iptr(0)-1; - ComplexMatrix tmp(n, 1); - - for (F77_INT i = 0; i < n; i++) - tmp(i,0) = ip2[P[i]]; - - lusolve (L, U, tmp); - - ip2 = workd+iptr(1)-1; - for (F77_INT i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } + // ido cannot be 2 for non-generalized problems + // see znaup2 + Complex *ip2 = workd+iptr(0)-1; + ComplexMatrix tmp(n, 1); + + for (F77_INT i = 0; i < n; i++) + tmp(i,0) = ip2[P[i]]; + + lusolve (L, U, tmp); + + ip2 = workd+iptr(1)-1; + for (F77_INT i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); } } else { if (info < 0) (*current_liboctave_error_handler) - ("eigs: error %d in dsaupd", info); + ("eigs: error %d in znaupd", info); break; } @@ -3005,6 +2982,10 @@ cresid(i) = Complex (rr(i),ri(i)); octave_rand::distribution (rand_dist); } + else if (n != cresid.numel ()) + { + (*current_liboctave_error_handler) ("eigs: opts.v0 must be n-by-1"); + } if (n < 3) (*current_liboctave_error_handler) @@ -3017,8 +2998,8 @@ if (p < 20) p = 20; - if (p > n - 1) - p = n - 1; + if (p > n) + p = n; } if (k <= 0 || k >= n - 1) @@ -3027,9 +3008,9 @@ " (must be 0 < k < n-1).\n" " Use 'eig (full (A))' instead"); - if (p <= k || p >= n) + if (p <= k || p > n) (*current_liboctave_error_handler) - ("eigs: opts.p must be greater than k and less than n"); + ("eigs: opts.p must be greater than k and less or equal to n"); if (! have_sigma) {
--- a/scripts/sparse/eigs.m Sat Nov 11 19:39:23 2017 +0100 +++ b/scripts/sparse/eigs.m Fri Mar 24 14:41:08 2017 +0100 @@ -124,7 +124,8 @@ ## @item p ## The number of Lanzcos basis vectors to use. More vectors will result in ## faster convergence, but a greater use of memory. The optimal value of -## @code{p} is problem dependent and should be in the range @var{k} to @var{n}. +## @code{p} is problem dependent and should be in the range +## @code{@var{k} + 1} to @var{n}. ## The default value is @code{2 * @var{k}}. ## ## @item v0 @@ -1229,6 +1230,12 @@ %! for i=1:k %! assert (max (abs ((A - d1(i)*eye (n))*v1(:,i))), 0, 1e-11); %! endfor +%!testif HAVE_ARPACK +%! [v1,d1] = eigs (A, k, "li"); +%! d1 = diag (d1); +%! for i=1:k +%! assert (max (abs ((A - d1(i)*eye (n))*v1(:,i))), 0, 1e-11); +%! endfor %!test %! A = 2 * diag (ones (10, 1)) - diag (ones (9, 1), 1) - diag (ones (9, 1), -1); @@ -1236,6 +1243,66 @@ %! reseig = eig (A, B); %! [~, idx] = sort (abs (reseig), "ascend"); %! assert (eigs (A, B, 10, 0), reseig (idx)); +%!testif HAVE_ARPACK +%! A = eye (9); +%! A(1, 1) = 0; +%! A(1, 9) = 1; +%! [V, L] = eigs (A, 4, -1); +%! assert (!any (isnan (diag (L)))) +%! assert (any (abs (diag (L)) <= 2 * eps)) +%!testif HAVE_ARPACK +%! A = diag (ones (9, 1), 1); +%! A(10,:) = [-1, zeros(1, 8), -1]; +%! opts.v0 = (1:10)'; +%! typ = "lr"; +%! [v, m] = eigs (A, 5, typ, opts); +%! assert (sort (real (diag (m))), ... +%! [-0.081751; 0.514038; 0.514038; 0.880290; 0.880290], 1e-4) +%! m = eigs (A, 5, typ, opts); +%! assert (sort (real (m)), ... +%! [-0.081751; 0.514038; 0.514038; 0.880290; 0.880290], 1e-4) +%! typ = "li"; +%! [v, m] = eigs (A, 5, typ, opts); +%! assert (sort (abs (imag (diag (m)))), ... +%! [0.75447; 0.78972; 0.78972; 0.96518; 0.96518], 1e-4) +%! m = eigs (A, 5, typ, opts); +%! assert (sort (abs (imag (m))), ... +%! [0.75447; 0.78972; 0.78972; 0.96518; 0.96518], 1e-4) +%! typ = "sr"; +%! [v, m] = eigs (A, 5, typ, opts); +%! assert (sort (real (diag (m))), ... +%! [-1.12180; -1.12180; -0.69077; -0.08175; -0.08175], 1e-4) +%! m = eigs (A, 5, typ, opts); +%! assert (sort (real (m)), ... +%! [-1.12180; -1.12180; -0.69077; -0.69077; -0.08175], 1e-4) +%! typ = "si"; +%! [v, m] = eigs (A, 5, typ, opts); +%! assert (sort (abs (imag (diag (m)))), ... +%! [0.25552; 0.25552; 0.30282; 0.30282; 0.75447], 1e-4) +%! m = eigs (A, 5, typ, opts); +%! assert (sort (abs (imag (m))), ... +%! [0.25552; 0.25552; 0.30282; 0.30282; 0.75447], 1e-4) +%! typ = "lm"; +%! [v, m] = eigs (A, 5, typ, opts); +%! assert (sort (abs (diag (m))), ... +%! [0.96863; 0.96863; 1.02294; 1.15054; 1.15054], 1e-4) +%! m = eigs (A, 5, typ, opts); +%! assert (sort (abs (m)), ... +%! [0.96863; 1.02294; 1.02294; 1.15054; 1.15054], 1e-4) +%! typ = "sm"; +%! [v, m] = eigs (A, 5, typ, opts); +%! assert (sort (abs (diag (m))), ... +%! [0.93092; 0.93092; 0.94228; 0.94228; 0.96863], 1e-4) +%! m = eigs (A, 5, typ, opts); +%! assert (sort (abs (m)), ... +%! [0.93092; 0.93092; 0.94228; 0.94228; 0.96863], 1e-4) +%!testif HAVE_ARPACK +%! A = toeplitz (sparse ([2, 1, zeros(1,8)])); +%! opts.v0 = (1:10)'; +%! [v, m] = eigs (A, 3, "sa", opts); +%! assert (diag (m), [0.081014; 0.317493; 0.690279], 1e-4); +%! m = eigs (A, 3, "sa", opts); +%! assert (m, [0.081014; 0.317493; 0.690279], 1e-4); %!test %! X = [70 47 42 39 50 73 79 23; @@ -1259,3 +1326,7 @@ %!assert (eigs (diag (1:5), 5, "sa"), [1;2;3;4;5]) %!assert (eigs (diag (1:5), 5, "la"), [5;4;3;2;1]) %!assert (eigs (diag (1:5), 3, "be"), [1;4;5]) +%!error +%! A = rand (10); +%! opts.v0 = ones (8, 1); +%! eigs (A, 4, "sm", opts);