Mercurial > hg > octave-lojdl
changeset 17226:ea9992fd9c89
fix eigs to handle small matrices
* __eigs__.cc: Rename from eigs.cc.
(F__eigs__): Rename from Feigs.
* dldfcn/module-files: Update for file rename.
* eigs.m: New file. Move eigs documentation and tests here from
eigs.cc. Handle small matrices here by calling eig then sorting and
selecting values as needed.
* scripts/sparse/module.mk (sparse_FCN_FILES): Add eigs.m to the list.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Mon, 12 Aug 2013 15:43:49 -0400 |
| parents | 33ce8c381f2c |
| children | a594e0d980eb |
| files | libinterp/dldfcn/__eigs__.cc libinterp/dldfcn/eigs.cc libinterp/dldfcn/module-files scripts/sparse/eigs.m scripts/sparse/module.mk |
| diffstat | 4 files changed, 1132 insertions(+), 926 deletions(-) [+] |
line wrap: on
line diff
rename from libinterp/dldfcn/eigs.cc rename to libinterp/dldfcn/__eigs__.cc --- a/libinterp/dldfcn/eigs.cc +++ b/libinterp/dldfcn/__eigs__.cc @@ -129,179 +129,29 @@ return retval; } -DEFUN_DLD (eigs, args, nargout, +DEFUN_DLD (__eigs__, args, nargout, "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{d} =} eigs (@var{A})\n\ -@deftypefnx {Loadable Function} {@var{d} =} eigs (@var{A}, @var{k})\n\ -@deftypefnx {Loadable Function} {@var{d} =} eigs (@var{A}, @var{k}, @var{sigma})\n\ -@deftypefnx {Loadable Function} {@var{d} =} eigs (@var{A}, @var{k}, @var{sigma}, @var{opts})\n\ -@deftypefnx {Loadable Function} {@var{d} =} eigs (@var{A}, @var{B})\n\ -@deftypefnx {Loadable Function} {@var{d} =} eigs (@var{A}, @var{B}, @var{k})\n\ -@deftypefnx {Loadable Function} {@var{d} =} eigs (@var{A}, @var{B}, @var{k}, @var{sigma})\n\ -@deftypefnx {Loadable Function} {@var{d} =} eigs (@var{A}, @var{B}, @var{k}, @var{sigma}, @var{opts})\n\ -@deftypefnx {Loadable Function} {@var{d} =} eigs (@var{af}, @var{n})\n\ -@deftypefnx {Loadable Function} {@var{d} =} eigs (@var{af}, @var{n}, @var{B})\n\ -@deftypefnx {Loadable Function} {@var{d} =} eigs (@var{af}, @var{n}, @var{k})\n\ -@deftypefnx {Loadable Function} {@var{d} =} eigs (@var{af}, @var{n}, @var{B}, @var{k})\n\ -@deftypefnx {Loadable Function} {@var{d} =} eigs (@var{af}, @var{n}, @var{k}, @var{sigma})\n\ -@deftypefnx {Loadable Function} {@var{d} =} eigs (@var{af}, @var{n}, @var{B}, @var{k}, @var{sigma})\n\ -@deftypefnx {Loadable Function} {@var{d} =} eigs (@var{af}, @var{n}, @var{k}, @var{sigma}, @var{opts})\n\ -@deftypefnx {Loadable Function} {@var{d} =} eigs (@var{af}, @var{n}, @var{B}, @var{k}, @var{sigma}, @var{opts})\n\ -@deftypefnx {Loadable Function} {[@var{V}, @var{d}] =} eigs (@var{A}, @dots{})\n\ -@deftypefnx {Loadable Function} {[@var{V}, @var{d}] =} eigs (@var{af}, @var{n}, @dots{})\n\ -@deftypefnx {Loadable Function} {[@var{V}, @var{d}, @var{flag}] =} eigs (@var{A}, @dots{})\n\ -@deftypefnx {Loadable Function} {[@var{V}, @var{d}, @var{flag}] =} eigs (@var{af}, @var{n}, @dots{})\n\ -Calculate a limited number of eigenvalues and eigenvectors of @var{A},\n\ -based on a selection criteria. The number of eigenvalues and eigenvectors to\n\ -calculate is given by @var{k} and defaults to 6.\n\ -\n\ -By default, @code{eigs} solve the equation\n\ -@tex\n\ -$A \\nu = \\lambda \\nu$,\n\ -@end tex\n\ -@ifinfo\n\ -@code{A * v = lambda * v},\n\ -@end ifinfo\n\ -where\n\ -@tex\n\ -$\\lambda$ is a scalar representing one of the eigenvalues, and $\\nu$\n\ -@end tex\n\ -@ifinfo\n\ -@code{lambda} is a scalar representing one of the eigenvalues, and @code{v}\n\ -@end ifinfo\n\ -is the corresponding eigenvector. If given the positive definite matrix\n\ -@var{B} then @code{eigs} solves the general eigenvalue equation\n\ -@tex\n\ -$A \\nu = \\lambda B \\nu$.\n\ -@end tex\n\ -@ifinfo\n\ -@code{A * v = lambda * B * v}.\n\ -@end ifinfo\n\ -\n\ -The argument @var{sigma} determines which eigenvalues are returned.\n\ -@var{sigma} can be either a scalar or a string. When @var{sigma} is a\n\ -scalar, the @var{k} eigenvalues closest to @var{sigma} are returned. If\n\ -@var{sigma} is a string, it must have one of the following values.\n\ -\n\ -@table @asis\n\ -@item \"lm\"\n\ -Largest Magnitude (default).\n\ -\n\ -@item \"sm\"\n\ -Smallest Magnitude.\n\ -\n\ -@item \"la\"\n\ -Largest Algebraic (valid only for real symmetric problems).\n\ -\n\ -@item \"sa\"\n\ -Smallest Algebraic (valid only for real symmetric problems).\n\ -\n\ -@item \"be\"\n\ -Both Ends, with one more from the high-end if @var{k} is odd (valid only for\n\ -real symmetric problems).\n\ -\n\ -@item \"lr\"\n\ -Largest Real part (valid only for complex or unsymmetric problems).\n\ -\n\ -@item \"sr\"\n\ -Smallest Real part (valid only for complex or unsymmetric problems).\n\ -\n\ -@item \"li\"\n\ -Largest Imaginary part (valid only for complex or unsymmetric problems).\n\ -\n\ -@item \"si\"\n\ -Smallest Imaginary part (valid only for complex or unsymmetric problems).\n\ -@end table\n\ -\n\ -If @var{opts} is given, it is a structure defining possible options that\n\ -@code{eigs} should use. The fields of the @var{opts} structure are:\n\ -\n\ -@table @code\n\ -@item issym\n\ -If @var{af} is given, then flags whether the function @var{af} defines a\n\ -symmetric problem. It is ignored if @var{A} is given. The default is false.\n\ -\n\ -@item isreal\n\ -If @var{af} is given, then flags whether the function @var{af} defines a\n\ -real problem. It is ignored if @var{A} is given. The default is true.\n\ -\n\ -@item tol\n\ -Defines the required convergence tolerance, calculated as\n\ -@code{tol * norm (A)}. The default is @code{eps}.\n\ -\n\ -@item maxit\n\ -The maximum number of iterations. The default is 300.\n\ -\n\ -@item p\n\ -The number of Lanzcos basis vectors to use. More vectors will result in\n\ -faster convergence, but a greater use of memory. The optimal value of\n\ -@code{p} is problem dependent and should be in the range @var{k} to @var{n}.\n\ -The default value is @code{2 * @var{k}}.\n\ -\n\ -@item v0\n\ -The starting vector for the algorithm. An initial vector close to the\n\ -final vector will speed up convergence. The default is for @sc{arpack}\n\ -to randomly generate a starting vector. If specified, @code{v0} must be\n\ -an @var{n}-by-1 vector where @code{@var{n} = rows (@var{A})}\n\ -\n\ -@item disp\n\ -The level of diagnostic printout (0|1|2). If @code{disp} is 0 then\n\ -diagnostics are disabled. The default value is 0.\n\ -\n\ -@item cholB\n\ -Flag if @code{chol (@var{B})} is passed rather than @var{B}. The default is\n\ -false.\n\ -\n\ -@item permB\n\ -The permutation vector of the Cholesky@tie{}factorization of @var{B} if\n\ -@code{cholB} is true. That is @code{chol (@var{B}(permB, permB))}. The\n\ -default is @code{1:@var{n}}.\n\ -\n\ -@end table\n\ -\n\ -It is also possible to represent @var{A} by a function denoted @var{af}.\n\ -@var{af} must be followed by a scalar argument @var{n} defining the length\n\ -of the vector argument accepted by @var{af}. @var{af} can be\n\ -a function handle, an inline function, or a string. When @var{af} is a\n\ -string it holds the name of the function to use.\n\ -\n\ -@var{af} is a function of the form @code{y = af (x)}\n\ -where the required return value of @var{af} is determined by\n\ -the value of @var{sigma}. The four possible forms are\n\ -\n\ -@table @code\n\ -@item A * x\n\ -if @var{sigma} is not given or is a string other than \"sm\".\n\ -\n\ -@item A \\ x\n\ -if @var{sigma} is 0 or \"sm\".\n\ -\n\ -@item (A - sigma * I) \\ x\n\ -for the standard eigenvalue problem, where @code{I} is the identity matrix of\n\ -the same size as @var{A}.\n\ -\n\ -@item (A - sigma * B) \\ x\n\ -for the general eigenvalue problem.\n\ -@end table\n\ -\n\ -The return arguments of @code{eigs} depend on the number of return arguments\n\ -requested. With a single return argument, a vector @var{d} of length @var{k}\n\ -is returned containing the @var{k} eigenvalues that have been found. With\n\ -two return arguments, @var{V} is a @var{n}-by-@var{k} matrix whose columns\n\ -are the @var{k} eigenvectors corresponding to the returned eigenvalues. The\n\ -eigenvalues themselves are returned in @var{d} in the form of a\n\ -@var{n}-by-@var{k} matrix, where the elements on the diagonal are the\n\ -eigenvalues.\n\ -\n\ -Given a third return argument @var{flag}, @code{eigs} returns the status\n\ -of the convergence. If @var{flag} is 0 then all eigenvalues have converged.\n\ -Any other value indicates a failure to converge.\n\ -\n\ -This function is based on the @sc{arpack} package, written by R. Lehoucq,\n\ -K. Maschhoff, D. Sorensen, and C. Yang. For more information see\n\ -@url{http://www.caam.rice.edu/software/ARPACK/}.\n\ -\n\ -@seealso{eig, svds}\n\ +@deftypefn {Loadable Function} {@var{d} =} __eigs__ (@var{A})\n\ +@deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{A}, @var{k})\n\ +@deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{A}, @var{k}, @var{sigma})\n\ +@deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{A}, @var{k}, @var{sigma}, @var{opts})\n\ +@deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{A}, @var{B})\n\ +@deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{A}, @var{B}, @var{k})\n\ +@deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{A}, @var{B}, @var{k}, @var{sigma})\n\ +@deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{A}, @var{B}, @var{k}, @var{sigma}, @var{opts})\n\ +@deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{af}, @var{n})\n\ +@deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{B})\n\ +@deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{k})\n\ +@deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{B}, @var{k})\n\ +@deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{k}, @var{sigma})\n\ +@deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{B}, @var{k}, @var{sigma})\n\ +@deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{k}, @var{sigma}, @var{opts})\n\ +@deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{B}, @var{k}, @var{sigma}, @var{opts})\n\ +@deftypefnx {Loadable Function} {[@var{V}, @var{d}] =} __eigs__ (@var{A}, @dots{})\n\ +@deftypefnx {Loadable Function} {[@var{V}, @var{d}] =} __eigs__ (@var{af}, @var{n}, @dots{})\n\ +@deftypefnx {Loadable Function} {[@var{V}, @var{d}, @var{flag}] =} __eigs__ (@var{A}, @dots{})\n\ +@deftypefnx {Loadable Function} {[@var{V}, @var{d}, @var{flag}] =} __eigs__ (@var{af}, @var{n}, @dots{})\n\ +Undocumented internal function.\n\ @end deftypefn") { octave_value_list retval; @@ -766,756 +616,3 @@ return retval; } - -/* #### SPARSE MATRIX VERSIONS #### */ - -/* -## Real positive definite tests, n must be even -%!shared n, k, A, d0, d2 -%! n = 20; -%! k = 4; -%! A = sparse ([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),4*ones(1,n),ones(1,n-2)]); -%! d0 = eig (A); -%! d2 = sort (d0); -%! [~, idx] = sort (abs (d0)); -%! d0 = d0(idx); -%! rand ("state", 42); # initialize generator to make eigs behavior reproducible -%!testif HAVE_ARPACK -%! d1 = eigs (A, k); -%! assert (d1, d0(end:-1:(end-k+1)), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k+1); -%! assert (d1, d0(end:-1:(end-k)), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "lm"); -%! assert (d1, d0(end:-1:(end-k+1)), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! d1 = eigs (A, k, "sm"); -%! assert (d1, d0(k:-1:1), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "la"); -%! assert (d1, d2(end:-1:(end-k+1)), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "sa"); -%! assert (d1, d2(1:k), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "be"); -%! assert (d1, d2([1:floor(k/2), (end - ceil(k/2) + 1):end]), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k+1, "be"); -%! assert (d1, d2([1:floor((k+1)/2), (end - ceil((k+1)/2) + 1):end]), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! d1 = eigs (A, k, 4.1); -%! [~, idx0] = sort (abs (d0 - 4.1)); -%! [~, idx1] = sort (abs (d1 - 4.1)); -%! assert (d1(idx1), d0(idx0(1:k)), 1e-11); -%!testif HAVE_ARPACK, HAVE_CHOLMOD -%! d1 = eigs (A, speye (n), k, "lm"); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! assert (eigs (A, k, 4.1), eigs (A, speye (n), k, 4.1), 1e-11); -%!testif HAVE_ARPACK -%! opts.cholB = true; -%! d1 = eigs (A, speye (n), k, "lm", opts); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! opts.cholB = true; -%! q = [2:n,1]; -%! opts.permB = q; -%! d1 = eigs (A, speye (n)(q,q), k, "lm", opts); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! opts.cholB = true; -%! d1 = eigs (A, speye (n), k, 4.1, opts); -%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! opts.cholB = true; -%! q = [2:n,1]; -%! opts.permB = q; -%! d1 = eigs (A, speye (n)(q,q), k, 4.1, opts); -%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! assert (eigs (A, k, 4.1), eigs (A, speye (n), k, 4.1), 1e-11); -%!testif HAVE_ARPACK -%! fn = @(x) A * x; -%! opts.issym = 1; opts.isreal = 1; -%! d1 = eigs (fn, n, k, "lm", opts); -%! assert (d1, d0(end:-1:(end-k+1)), 1e-11); -%!testif HAVE_ARPACK -%! fn = @(x) A \ x; -%! opts.issym = 1; opts.isreal = 1; -%! d1 = eigs (fn, n, k, "sm", opts); -%! assert (d1, d0(k:-1:1), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! fn = @(x) (A - 4.1 * eye (n)) \ x; -%! opts.issym = 1; opts.isreal = 1; -%! d1 = eigs (fn, n, k, 4.1, opts); -%! assert (d1, eigs (A, k, 4.1), 1e-11); -%!testif HAVE_ARPACK -%! AA = speye (10); -%! fn = @(x) AA * x; -%! opts.issym = 1; opts.isreal = 1; -%! assert (eigs (fn, 10, AA, 3, "lm", opts), [1; 1; 1], 10*eps); -%!testif HAVE_ARPACK -%! [v1,d1] = eigs (A, k, "lm"); -%! d1 = diag (d1); -%! for i=1:k -%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); -%! endfor -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! [v1,d1] = eigs (A, k, "sm"); -%! d1 = diag (d1); -%! for i=1:k -%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); -%! endfor -%!testif HAVE_ARPACK -%! [v1,d1] = eigs (A, k, "la"); -%! d1 = diag (d1); -%! for i=1:k -%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); -%! endfor -%!testif HAVE_ARPACK -%! [v1,d1] = eigs (A, k, "sa"); -%! d1 = diag (d1); -%! for i=1:k -%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); -%! endfor -%!testif HAVE_ARPACK -%! [v1,d1] = eigs (A, k, "be"); -%! d1 = diag (d1); -%! for i=1:k -%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); -%! endfor -*/ - -/* -## Real unsymmetric tests -%!shared n, k, A, d0 -%! n = 20; -%! k = 4; -%! A = sparse ([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),1:n,-ones(1,n-2)]); -%! d0 = eig (A); -%! [~, idx] = sort (abs (d0)); -%! d0 = d0(idx); -%! rand ("state", 42); % initialize generator to make eigs behavior reproducible -%!testif HAVE_ARPACK -%! d1 = eigs (A, k); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k+1); -%! assert (abs (d1), abs (d0(end:-1:(end-k))),1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "lm"); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! d1 = eigs (A, k, "sm"); -%! assert (abs (d1), abs (d0(1:k)), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "lr"); -%! [~, idx] = sort (real (d0)); -%! d2 = d0(idx); -%! assert (real (d1), real (d2(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "sr"); -%! [~, idx] = sort (real (abs (d0))); -%! d2 = d0(idx); -%! assert (real (d1), real (d2(1:k)), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "li"); -%! [~, idx] = sort (imag (abs (d0))); -%! d2 = d0(idx); -%! assert (sort (imag (d1)), sort (imag (d2(end:-1:(end-k+1)))), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "si"); -%! [~, idx] = sort (imag (abs (d0))); -%! d2 = d0(idx); -%! assert (sort (imag (d1)), sort (imag (d2(1:k))), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! d1 = eigs (A, k, 4.1); -%! [~, idx0] = sort (abs (d0 - 4.1)); -%! [~, idx1] = sort (abs (d1 - 4.1)); -%! assert (abs (d1(idx1)), abs (d0(idx0(1:k))), 1e-11); -%! assert (sort (imag (d1(idx1))), sort (imag (d0(idx0(1:k)))), 1e-11); -%!testif HAVE_ARPACK, HAVE_CHOLMOD -%! d1 = eigs (A, speye (n), k, "lm"); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! opts.cholB = true; -%! d1 = eigs (A, speye (n), k, "lm", opts); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! opts.cholB = true; -%! q = [2:n,1]; -%! opts.permB = q; -%! d1 = eigs (A, speye (n)(q,q), k, "lm", opts); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! opts.cholB = true; -%! d1 = eigs (A, speye (n), k, 4.1, opts); -%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! opts.cholB = true; -%! q = [2:n,1]; -%! opts.permB = q; -%! d1 = eigs (A, speye (n)(q,q), k, 4.1, opts); -%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! assert (abs (eigs (A, k, 4.1)), abs (eigs (A, speye (n), k, 4.1)), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! assert (sort (imag (eigs (A, k, 4.1))), sort (imag (eigs (A, speye (n), k, 4.1))), 1e-11); -%!testif HAVE_ARPACK -%! fn = @(x) A * x; -%! opts.issym = 0; opts.isreal = 1; -%! d1 = eigs (fn, n, k, "lm", opts); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! fn = @(x) A \ x; -%! opts.issym = 0; opts.isreal = 1; -%! d1 = eigs (fn, n, k, "sm", opts); -%! assert (abs (d1), d0(1:k), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! fn = @(x) (A - 4.1 * eye (n)) \ x; -%! opts.issym = 0; opts.isreal = 1; -%! d1 = eigs (fn, n, k, 4.1, opts); -%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); -%!testif HAVE_ARPACK -%! [v1,d1] = eigs (A, k, "lm"); -%! d1 = diag (d1); -%! for i=1:k -%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); -%! endfor -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! [v1,d1] = eigs (A, k, "sm"); -%! d1 = diag (d1); -%! for i=1:k -%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); -%! endfor -%!testif HAVE_ARPACK -%! [v1,d1] = eigs (A, k, "lr"); -%! d1 = diag (d1); -%! for i=1:k -%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); -%! endfor -%!testif HAVE_ARPACK -%! [v1,d1] = eigs (A, k, "sr"); -%! d1 = diag (d1); -%! for i=1:k -%! assert (max (abs ((A - d1(i)*speye (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)*speye (n))*v1(:,i))), 0, 1e-11); -%! endfor -%!testif HAVE_ARPACK -%! [v1,d1] = eigs (A, k, "si"); -%! d1 = diag (d1); -%! for i=1:k -%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); -%! endfor -*/ - -/* -## Complex hermitian tests -%!shared n, k, A, d0 -%! n = 20; -%! k = 4; -%! A = sparse ([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[1i*ones(1,n-2),4*ones(1,n),-1i*ones(1,n-2)]); -%! d0 = eig (A); -%! [~, idx] = sort (abs (d0)); -%! d0 = d0(idx); -%! rand ("state", 42); % initialize generator to make eigs behavior reproducible -%!testif HAVE_ARPACK -%! d1 = eigs (A, k); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k+1); -%! assert (abs (d1), abs (d0(end:-1:(end-k))),1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "lm"); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! d1 = eigs (A, k, "sm"); -%! assert (abs (d1), abs (d0(1:k)), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "lr"); -%! [~, idx] = sort (real (abs (d0))); -%! d2 = d0(idx); -%! assert (real (d1), real (d2(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "sr"); -%! [~, idx] = sort (real (abs (d0))); -%! d2 = d0(idx); -%! assert (real (d1), real (d2(1:k)), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "li"); -%! [~, idx] = sort (imag (abs (d0))); -%! d2 = d0(idx); -%! assert (sort (imag (d1)), sort (imag (d2(end:-1:(end-k+1)))), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "si"); -%! [~, idx] = sort (imag (abs (d0))); -%! d2 = d0(idx); -%! assert (sort (imag (d1)), sort (imag (d2(1:k))), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! d1 = eigs (A, k, 4.1); -%! [~, idx0] = sort (abs (d0 - 4.1)); -%! [~, idx1] = sort (abs (d1 - 4.1)); -%! assert (abs (d1(idx1)), abs (d0(idx0(1:k))), 1e-11); -%! assert (sort (imag (d1(idx1))), sort (imag (d0(idx0(1:k)))), 1e-11); -%!testif HAVE_ARPACK, HAVE_CHOLMOD -%! d1 = eigs (A, speye (n), k, "lm"); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! opts.cholB = true; -%! d1 = eigs (A, speye (n), k, "lm", opts); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! opts.cholB = true; -%! q = [2:n,1]; -%! opts.permB = q; -%! d1 = eigs (A, speye (n)(q,q), k, "lm", opts); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! opts.cholB = true; -%! d1 = eigs (A, speye (n), k, 4.1, opts); -%! assert (abs (abs (d1)), abs (eigs (A, k, 4.1)), 1e-11); -%! assert (sort (imag (abs (d1))), sort (imag (eigs (A, k, 4.1))), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! opts.cholB = true; -%! q = [2:n,1]; -%! opts.permB = q; -%! d1 = eigs (A, speye (n)(q,q), k, 4.1, opts); -%! assert (abs (abs (d1)), abs (eigs (A, k, 4.1)), 1e-11); -%! assert (sort (imag (abs (d1))), sort (imag (eigs (A, k, 4.1))), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! assert (abs (eigs (A, k, 4.1)), abs (eigs (A, speye (n), k, 4.1)), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! assert (sort (imag (eigs (A, k, 4.1))), sort (imag (eigs (A, speye (n), k, 4.1))), 1e-11); -%!testif HAVE_ARPACK -%! fn = @(x) A * x; -%! opts.issym = 0; opts.isreal = 0; -%! d1 = eigs (fn, n, k, "lm", opts); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! fn = @(x) A \ x; -%! opts.issym = 0; opts.isreal = 0; -%! d1 = eigs (fn, n, k, "sm", opts); -%! assert (abs (d1), d0(1:k), 1e-11); -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! fn = @(x) (A - 4.1 * eye (n)) \ x; -%! opts.issym = 0; opts.isreal = 0; -%! d1 = eigs (fn, n, k, 4.1, opts); -%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); -%!testif HAVE_ARPACK -%! [v1,d1] = eigs (A, k, "lm"); -%! d1 = diag (d1); -%! for i=1:k -%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); -%! endfor -%!testif HAVE_ARPACK, HAVE_UMFPACK -%! [v1,d1] = eigs (A, k, "sm"); -%! d1 = diag (d1); -%! for i=1:k -%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); -%! endfor -%!testif HAVE_ARPACK -%! [v1,d1] = eigs (A, k, "lr"); -%! d1 = diag (d1); -%! for i=1:k -%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); -%! endfor -%!testif HAVE_ARPACK -%! [v1,d1] = eigs (A, k, "sr"); -%! d1 = diag (d1); -%! for i=1:k -%! assert (max (abs ((A - d1(i)*speye (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)*speye (n))*v1(:,i))), 0, 1e-11); -%! endfor -%!testif HAVE_ARPACK -%! [v1,d1] = eigs (A, k, "si"); -%! d1 = diag (d1); -%! for i=1:k -%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); -%! endfor -*/ - -/* #### FULL MATRIX VERSIONS #### */ - -/* -## Real positive definite tests, n must be even -%!shared n, k, A, d0, d2 -%! n = 20; -%! k = 4; -%! A = full (sparse ([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),4*ones(1,n),ones(1,n-2)])); -%! d0 = eig (A); -%! d2 = sort (d0); -%! [~, idx] = sort (abs (d0)); -%! d0 = d0(idx); -%! rand ("state", 42); % initialize generator to make eigs behavior reproducible -%!testif HAVE_ARPACK -%! d1 = eigs (A, k); -%! assert (d1, d0(end:-1:(end-k+1)), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k+1); -%! assert (d1, d0(end:-1:(end-k)),1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "lm"); -%! assert (d1, d0(end:-1:(end-k+1)), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "sm"); -%! assert (d1, d0(k:-1:1), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "la"); -%! assert (d1, d2(end:-1:(end-k+1)), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "sa"); -%! assert (d1, d2(1:k), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "be"); -%! assert (d1, d2([1:floor(k/2), (end - ceil(k/2) + 1):end]), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k+1, "be"); -%! assert (d1, d2([1:floor((k+1)/2), (end - ceil((k+1)/2) + 1):end]), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, 4.1); -%! [~, idx0] = sort (abs (d0 - 4.1)); -%! [~, idx1] = sort (abs (d1 - 4.1)); -%! assert (d1(idx1), d0(idx0(1:k)), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, eye (n), k, "lm"); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! assert (eigs (A, k, 4.1), eigs (A, eye (n), k, 4.1), 1e-11); -%!testif HAVE_ARPACK -%! opts.cholB = true; -%! d1 = eigs (A, eye (n), k, "lm", opts); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! opts.cholB = true; -%! q = [2:n,1]; -%! opts.permB = q; -%! d1 = eigs (A, eye (n)(q,q), k, "lm", opts); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! opts.cholB = true; -%! d1 = eigs (A, eye (n), k, 4.1, opts); -%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); -%!testif HAVE_ARPACK -%! opts.cholB = true; -%! q = [2:n,1]; -%! opts.permB = q; -%! d1 = eigs (A, eye (n)(q,q), k, 4.1, opts); -%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); -%!testif HAVE_ARPACK -%! assert (eigs (A, k, 4.1), eigs (A, eye (n), k, 4.1), 1e-11); -%!testif HAVE_ARPACK -%! fn = @(x) A * x; -%! opts.issym = 1; opts.isreal = 1; -%! d1 = eigs (fn, n, k, "lm", opts); -%! assert (d1, d0(end:-1:(end-k+1)), 1e-11); -%!testif HAVE_ARPACK -%! fn = @(x) A \ x; -%! opts.issym = 1; opts.isreal = 1; -%! d1 = eigs (fn, n, k, "sm", opts); -%! assert (d1, d0(k:-1:1), 1e-11); -%!testif HAVE_ARPACK -%! fn = @(x) (A - 4.1 * eye (n)) \ x; -%! opts.issym = 1; opts.isreal = 1; -%! d1 = eigs (fn, n, k, 4.1, opts); -%! assert (d1, eigs (A, k, 4.1), 1e-11); -%!testif HAVE_ARPACK -%! [v1,d1] = eigs (A, k, "lm"); -%! d1 = diag (d1); -%! 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, "sm"); -%! d1 = diag (d1); -%! 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, "la"); -%! d1 = diag (d1); -%! 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, "sa"); -%! d1 = diag (d1); -%! 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, "be"); -%! d1 = diag (d1); -%! for i=1:k -%! assert (max (abs ((A - d1(i)*eye (n))*v1(:,i))), 0, 1e-11); -%! endfor -*/ - -/* -## Real unsymmetric tests -%!shared n, k, A, d0 -%! n = 20; -%! k = 4; -%! A = full (sparse ([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),1:n,-ones(1,n-2)])); -%! d0 = eig (A); -%! [~, idx] = sort (abs (d0)); -%! d0 = d0(idx); -%! rand ("state", 42); % initialize generator to make eigs behavior reproducible -%!testif HAVE_ARPACK -%! d1 = eigs (A, k); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k+1); -%! assert (abs (d1), abs (d0(end:-1:(end-k))),1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "lm"); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "sm"); -%! assert (abs (d1), abs (d0(1:k)), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "lr"); -%! [~, idx] = sort (real (d0)); -%! d2 = d0(idx); -%! assert (real (d1), real (d2(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "sr"); -%! [~, idx] = sort (real (abs (d0))); -%! d2 = d0(idx); -%! assert (real (d1), real (d2(1:k)), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "li"); -%! [~, idx] = sort (imag (abs (d0))); -%! d2 = d0(idx); -%! assert (sort (imag (d1)), sort (imag (d2(end:-1:(end-k+1)))), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "si"); -%! [~, idx] = sort (imag (abs (d0))); -%! d2 = d0(idx); -%! assert (sort (imag (d1)), sort (imag (d2(1:k))), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, 4.1); -%! [~, idx0] = sort (abs (d0 - 4.1)); -%! [~, idx1] = sort (abs (d1 - 4.1)); -%! assert (abs (d1(idx1)), abs (d0(idx0(1:k))), 1e-11); -%! assert (sort (imag (d1(idx1))), sort (imag (d0(idx0(1:k)))), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, eye (n), k, "lm"); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! opts.cholB = true; -%! d1 = eigs (A, eye (n), k, "lm", opts); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! opts.cholB = true; -%! q = [2:n,1]; -%! opts.permB = q; -%! d1 = eigs (A, eye (n)(q,q), k, "lm", opts); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! opts.cholB = true; -%! d1 = eigs (A, eye (n), k, 4.1, opts); -%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); -%!testif HAVE_ARPACK -%! opts.cholB = true; -%! q = [2:n,1]; -%! opts.permB = q; -%! d1 = eigs (A, eye (n)(q,q), k, 4.1, opts); -%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); -%!testif HAVE_ARPACK -%! assert (abs (eigs (A, k, 4.1)), abs (eigs (A, eye (n), k, 4.1)), 1e-11); -%!testif HAVE_ARPACK -%! assert (sort (imag (eigs (A, k, 4.1))), sort (imag (eigs (A, eye (n), k, 4.1))), 1e-11); -%!testif HAVE_ARPACK -%! fn = @(x) A * x; -%! opts.issym = 0; opts.isreal = 1; -%! d1 = eigs (fn, n, k, "lm", opts); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! fn = @(x) A \ x; -%! opts.issym = 0; opts.isreal = 1; -%! d1 = eigs (fn, n, k, "sm", opts); -%! assert (abs (d1), d0(1:k), 1e-11); -%!testif HAVE_ARPACK -%! fn = @(x) (A - 4.1 * eye (n)) \ x; -%! opts.issym = 0; opts.isreal = 1; -%! d1 = eigs (fn, n, k, 4.1, opts); -%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); -%!testif HAVE_ARPACK -%! [v1,d1] = eigs (A, k, "lm"); -%! d1 = diag (d1); -%! 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, "sm"); -%! d1 = diag (d1); -%! 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, "lr"); -%! d1 = diag (d1); -%! 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, "sr"); -%! d1 = diag (d1); -%! 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 -%!testif HAVE_ARPACK -%! [v1,d1] = eigs (A, k, "si"); -%! d1 = diag (d1); -%! for i=1:k -%! assert (max (abs ((A - d1(i)*eye (n))*v1(:,i))), 0, 1e-11); -%! endfor -*/ - -/* -## Complex hermitian tests -%!shared n, k, A, d0 -%! n = 20; -%! k = 4; -%! A = full (sparse ([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[1i*ones(1,n-2),4*ones(1,n),-1i*ones(1,n-2)])); -%! d0 = eig (A); -%! [~, idx] = sort (abs (d0)); -%! d0 = d0(idx); -%! rand ("state", 42); % initialize generator to make eigs behavior reproducible -%!testif HAVE_ARPACK -%! d1 = eigs (A, k); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k+1); -%! assert (abs (d1), abs (d0(end:-1:(end-k))),1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "lm"); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "sm"); -%! assert (abs (d1), abs (d0(1:k)), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "lr"); -%! [~, idx] = sort (real (abs (d0))); -%! d2 = d0(idx); -%! assert (real (d1), real (d2(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "sr"); -%! [~, idx] = sort (real (abs (d0))); -%! d2 = d0(idx); -%! assert (real (d1), real (d2(1:k)), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "li"); -%! [~, idx] = sort (imag (abs (d0))); -%! d2 = d0(idx); -%! assert (sort (imag (d1)), sort (imag (d2(end:-1:(end-k+1)))), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, "si"); -%! [~, idx] = sort (imag (abs (d0))); -%! d2 = d0(idx); -%! assert (sort (imag (d1)), sort (imag (d2(1:k))), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, k, 4.1); -%! [~, idx0] = sort (abs (d0 - 4.1)); -%! [~, idx1] = sort (abs (d1 - 4.1)); -%! assert (abs (d1(idx1)), abs (d0(idx0(1:k))), 1e-11); -%! assert (sort (imag (d1(idx1))), sort (imag (d0(idx0(1:k)))), 1e-11); -%!testif HAVE_ARPACK -%! d1 = eigs (A, eye (n), k, "lm"); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! opts.cholB = true; -%! d1 = eigs (A, eye (n), k, "lm", opts); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! opts.cholB = true; -%! q = [2:n,1]; -%! opts.permB = q; -%! d1 = eigs (A, eye (n)(q,q), k, "lm", opts); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! opts.cholB = true; -%! d1 = eigs (A, eye (n), k, 4.1, opts); -%! assert (abs (abs (d1)), abs (eigs (A, k, 4.1)), 1e-11); -%! assert (sort (imag (abs (d1))), sort (imag (eigs (A, k, 4.1))), 1e-11); -%!testif HAVE_ARPACK -%! opts.cholB = true; -%! q = [2:n,1]; -%! opts.permB = q; -%! d1 = eigs (A, eye (n)(q,q), k, 4.1, opts); -%! assert (abs (abs (d1)), abs (eigs (A, k, 4.1)), 1e-11); -%! assert (sort (imag (abs (d1))), sort (imag (eigs (A, k, 4.1))), 1e-11); -%!testif HAVE_ARPACK -%! assert (abs (eigs (A, k, 4.1)), abs (eigs (A, eye (n), k, 4.1)), 1e-11); -%!testif HAVE_ARPACK -%! assert (sort (imag (eigs (A, k, 4.1))), sort (imag (eigs (A, eye (n), k, 4.1))), 1e-11); -%!testif HAVE_ARPACK -%! fn = @(x) A * x; -%! opts.issym = 0; opts.isreal = 0; -%! d1 = eigs (fn, n, k, "lm", opts); -%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); -%!testif HAVE_ARPACK -%! fn = @(x) A \ x; -%! opts.issym = 0; opts.isreal = 0; -%! d1 = eigs (fn, n, k, "sm", opts); -%! assert (abs (d1), d0(1:k), 1e-11); -%!testif HAVE_ARPACK -%! fn = @(x) (A - 4.1 * eye (n)) \ x; -%! opts.issym = 0; opts.isreal = 0; -%! d1 = eigs (fn, n, k, 4.1, opts); -%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); -%!testif HAVE_ARPACK -%! [v1,d1] = eigs (A, k, "lm"); -%! d1 = diag (d1); -%! 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, "sm"); -%! d1 = diag (d1); -%! 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, "lr"); -%! d1 = diag (d1); -%! 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, "sr"); -%! d1 = diag (d1); -%! 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 -%!testif HAVE_ARPACK -%! [v1,d1] = eigs (A, k, "si"); -%! d1 = diag (d1); -%! for i=1:k -%! assert (max (abs ((A - d1(i)*eye (n))*v1(:,i))), 0, 1e-11); -%! endfor -*/
--- a/libinterp/dldfcn/module-files +++ b/libinterp/dldfcn/module-files @@ -1,6 +1,7 @@ # FILE|CPPFLAGS|LDFLAGS|LIBRARIES __delaunayn__.cc|$(QHULL_CPPFLAGS)|$(QHULL_LDFLAGS)|$(QHULL_LIBS) __dsearchn__.cc +__eigs__.cc|$(ARPACK_CPPFLAGS) $(SPARSE_XCPPFLAGS)|$(ARPACK_LDFLAGS) $(SPARSE_XLDFLAGS)|$(ARPACK_LIBS) $(SPARSE_XLIBS) $(LAPACK_LIBS) $(BLAS_LIBS) __fltk_uigetfile__.cc|$(GRAPHICS_CFLAGS) $(FT2_CPPFLAGS)|$(GRAPHICS_LDFLAGS) $(FT2_LDFLAGS)|$(GRAPHICS_LIBS) $(FT2_LIBS) __glpk__.cc|$(GLPK_CPPFLAGS)|$(GLPK_LDFLAGS)|$(GLPK_LIBS) __init_fltk__.cc|$(GRAPHICS_CFLAGS) $(FT2_CPPFLAGS)|$(GRAPHICS_LDFLAGS) $(FT2_LDFLAGS)|$(GRAPHICS_LIBS) $(FT2_LIBS) $(OPENGL_LIBS) @@ -13,7 +14,6 @@ colamd.cc|$(SPARSE_XCPPFLAGS)|$(SPARSE_XLDFLAGS)|$(SPARSE_XLIBS) convhulln.cc|$(QHULL_CPPFLAGS)|$(QHULL_LDFLAGS)|$(QHULL_LIBS) dmperm.cc|$(SPARSE_XCPPFLAGS)|$(SPARSE_XLDFLAGS)|$(SPARSE_XLIBS) -eigs.cc|$(ARPACK_CPPFLAGS) $(SPARSE_XCPPFLAGS)|$(ARPACK_LDFLAGS) $(SPARSE_XLDFLAGS)|$(ARPACK_LIBS) $(SPARSE_XLIBS) $(LAPACK_LIBS) $(BLAS_LIBS) fftw.cc|$(FFTW_XCPPFLAGS)|$(FFTW_XLDFLAGS)|$(FFTW_XLIBS) qr.cc|$(QRUPDATE_CPPFLAGS) $(SPARSE_XCPPFLAGS)|$(QRUPDATE_LDFLAGS) $(SPARSE_XLDFLAGS)|$(QRUPDATE_LIBS) $(SPARSE_XLIBS) symbfact.cc|$(SPARSE_XCPPFLAGS)|$(SPARSE_XLDFLAGS)|$(SPARSE_XLIBS)
new file mode 100644 --- /dev/null +++ b/scripts/sparse/eigs.m @@ -0,0 +1,1108 @@ +## Copyright (C) 2005-2012 David Bateman +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{d} =} eigs (@var{A}) +## @deftypefnx {Function File} {@var{d} =} eigs (@var{A}, @var{k}) +## @deftypefnx {Function File} {@var{d} =} eigs (@var{A}, @var{k}, @var{sigma}) +## @deftypefnx {Function File} {@var{d} =} eigs (@var{A}, @var{k}, @var{sigma}, @var{opts}) +## @deftypefnx {Function File} {@var{d} =} eigs (@var{A}, @var{B}) +## @deftypefnx {Function File} {@var{d} =} eigs (@var{A}, @var{B}, @var{k}) +## @deftypefnx {Function File} {@var{d} =} eigs (@var{A}, @var{B}, @var{k}, @var{sigma}) +## @deftypefnx {Function File} {@var{d} =} eigs (@var{A}, @var{B}, @var{k}, @var{sigma}, @var{opts}) +## @deftypefnx {Function File} {@var{d} =} eigs (@var{af}, @var{n}) +## @deftypefnx {Function File} {@var{d} =} eigs (@var{af}, @var{n}, @var{B}) +## @deftypefnx {Function File} {@var{d} =} eigs (@var{af}, @var{n}, @var{k}) +## @deftypefnx {Function File} {@var{d} =} eigs (@var{af}, @var{n}, @var{B}, @var{k}) +## @deftypefnx {Function File} {@var{d} =} eigs (@var{af}, @var{n}, @var{k}, @var{sigma}) +## @deftypefnx {Function File} {@var{d} =} eigs (@var{af}, @var{n}, @var{B}, @var{k}, @var{sigma}) +## @deftypefnx {Function File} {@var{d} =} eigs (@var{af}, @var{n}, @var{k}, @var{sigma}, @var{opts}) +## @deftypefnx {Function File} {@var{d} =} eigs (@var{af}, @var{n}, @var{B}, @var{k}, @var{sigma}, @var{opts}) +## @deftypefnx {Function File} {[@var{V}, @var{d}] =} eigs (@var{A}, @dots{}) +## @deftypefnx {Function File} {[@var{V}, @var{d}] =} eigs (@var{af}, @var{n}, @dots{}) +## @deftypefnx {Function File} {[@var{V}, @var{d}, @var{flag}] =} eigs (@var{A}, @dots{}) +## @deftypefnx {Function File} {[@var{V}, @var{d}, @var{flag}] =} eigs (@var{af}, @var{n}, @dots{}) +## Calculate a limited number of eigenvalues and eigenvectors of @var{A}, +## based on a selection criteria. The number of eigenvalues and eigenvectors to +## calculate is given by @var{k} and defaults to 6. +## +## By default, @code{eigs} solve the equation +## @tex +## $A \nu = \lambda \nu$, +## @end tex +## @ifinfo +## @code{A * v = lambda * v}, +## @end ifinfo +## where +## @tex +## $\lambda$ is a scalar representing one of the eigenvalues, and $\nu$ +## @end tex +## @ifinfo +## @code{lambda} is a scalar representing one of the eigenvalues, and @code{v} +## @end ifinfo +## is the corresponding eigenvector. If given the positive definite matrix +## @var{B} then @code{eigs} solves the general eigenvalue equation +## @tex +## $A \nu = \lambda B \nu$. +## @end tex +## @ifinfo +## @code{A * v = lambda * B * v}. +## @end ifinfo +## +## The argument @var{sigma} determines which eigenvalues are returned. +## @var{sigma} can be either a scalar or a string. When @var{sigma} is a +## scalar, the @var{k} eigenvalues closest to @var{sigma} are returned. If +## @var{sigma} is a string, it must have one of the following values. +## +## @table @asis +## @item "lm" +## Largest Magnitude (default). +## +## @item "sm" +## Smallest Magnitude. +## +## @item "la" +## Largest Algebraic (valid only for real symmetric problems). +## +## @item "sa" +## Smallest Algebraic (valid only for real symmetric problems). +## +## @item "be" +## Both Ends, with one more from the high-end if @var{k} is odd (valid only for +## real symmetric problems). +## +## @item "lr" +## Largest Real part (valid only for complex or unsymmetric problems). +## +## @item "sr" +## Smallest Real part (valid only for complex or unsymmetric problems). +## +## @item "li" +## Largest Imaginary part (valid only for complex or unsymmetric problems). +## +## @item "si" +## Smallest Imaginary part (valid only for complex or unsymmetric problems). +## @end table +## +## If @var{opts} is given, it is a structure defining possible options that +## @code{eigs} should use. The fields of the @var{opts} structure are: +## +## @table @code +## @item issym +## If @var{af} is given, then flags whether the function @var{af} defines a +## symmetric problem. It is ignored if @var{A} is given. The default is false. +## +## @item isreal +## If @var{af} is given, then flags whether the function @var{af} defines a +## real problem. It is ignored if @var{A} is given. The default is true. +## +## @item tol +## Defines the required convergence tolerance, calculated as +## @code{tol * norm (A)}. The default is @code{eps}. +## +## @item maxit +## The maximum number of iterations. The default is 300. +## +## @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}. +## The default value is @code{2 * @var{k}}. +## +## @item v0 +## The starting vector for the algorithm. An initial vector close to the +## final vector will speed up convergence. The default is for @sc{arpack} +## to randomly generate a starting vector. If specified, @code{v0} must be +## an @var{n}-by-1 vector where @code{@var{n} = rows (@var{A})} +## +## @item disp +## The level of diagnostic printout (0|1|2). If @code{disp} is 0 then +## diagnostics are disabled. The default value is 0. +## +## @item cholB +## Flag if @code{chol (@var{B})} is passed rather than @var{B}. The default is +## false. +## +## @item permB +## The permutation vector of the Cholesky@tie{}factorization of @var{B} if +## @code{cholB} is true. That is @code{chol (@var{B}(permB, permB))}. The +## default is @code{1:@var{n}}. +## +## @end table +## +## It is also possible to represent @var{A} by a function denoted @var{af}. +## @var{af} must be followed by a scalar argument @var{n} defining the length +## of the vector argument accepted by @var{af}. @var{af} can be +## a function handle, an inline function, or a string. When @var{af} is a +## string it holds the name of the function to use. +## +## @var{af} is a function of the form @code{y = af (x)} +## where the required return value of @var{af} is determined by +## the value of @var{sigma}. The four possible forms are +## +## @table @code +## @item A * x +## if @var{sigma} is not given or is a string other than "sm". +## +## @item A \ x +## if @var{sigma} is 0 or "sm". +## +## @item (A - sigma * I) \ x +## for the standard eigenvalue problem, where @code{I} is the identity matrix of +## the same size as @var{A}. +## +## @item (A - sigma * B) \ x +## for the general eigenvalue problem. +## @end table +## +## The return arguments of @code{eigs} depend on the number of return arguments +## requested. With a single return argument, a vector @var{d} of length @var{k} +## is returned containing the @var{k} eigenvalues that have been found. With +## two return arguments, @var{V} is a @var{n}-by-@var{k} matrix whose columns +## are the @var{k} eigenvectors corresponding to the returned eigenvalues. The +## eigenvalues themselves are returned in @var{d} in the form of a +## @var{n}-by-@var{k} matrix, where the elements on the diagonal are the +## eigenvalues. +## +## Given a third return argument @var{flag}, @code{eigs} returns the status +## of the convergence. If @var{flag} is 0 then all eigenvalues have converged. +## Any other value indicates a failure to converge. +## +## This function is based on the @sc{arpack} package, written by R. Lehoucq, +## K. Maschhoff, D. Sorensen, and C. Yang. For more information see +## @url{http://www.caam.rice.edu/software/ARPACK/}. +## +## @seealso{eig, svds} +## @end deftypefn + +function varargout = eigs (varargin) + + ## For compatibility with Matlab, handle small matrix cases here + ## that ARPACK does not. + + if (nargin == 0) + print_usage (); + endif + + call_eig = false; + offset = 0; + k = 6; + sigma = "lm"; + + if (isnumeric (varargin{1}) && issquare (varargin{1})) + a = varargin{1}; + if (nargin > 1 && isnumeric (varargin{2}) + && issquare (varargin{2}) && size_equal (a, varargin{2})) + b = varargin{2}; + offset = 1; + endif + + if (rows (a) < 9) + call_eig = true; + endif + + if (nargin > 1 + offset) + tmp = varargin{2+offset}; + if (isnumeric (tmp) && isscalar (tmp) && isreal (tmp) + && round (tmp) == tmp) + k = tmp; + + if (rows (a) - k < 3) + call_eig = true; + endif + else + call_eig = false; + endif + + if (nargin > 2 + offset) + tmp = varargin{3+offset}; + if (ischar (tmp) || (isnumeric (tmp) && isscalar (tmp))) + sigma = tmp; + else + call_eig = false; + endif + endif + endif + endif + + if (call_eig) + varargout = cell (1, min (2, max (1, nargout))); + if (offset) + real_valued = isreal (a) && isreal (b); + symmetric = issymmetric (a) && issymmetric (b); + [varargout{:}] = eig (a, b); + else + real_valued = isreal (a); + symmetric = issymmetric (a); + [varargout{:}] = eig (a); + endif + varargout = select (varargout, k, sigma, real_valued, symmetric); + if (nargout == 3) + varargout{3} = 0; + endif + else + varargout = cell (1, max (1, nargout)); + [varargout{:}] = __eigs__ (varargin{:}); + endif + +endfunction + +function out = select (args, k, sigma, real_valued, symmetric) + + if (numel (args) == 1) + d = args{1}; + else + d = diag (args{2}); + endif + + if (ischar (sigma)) + switch (sigma) + case "lm" + [~, idx] = sort (abs (d), "descend"); + + case "sm" + [~, idx] = sort (abs (d), "ascend"); + + case "la" + if (real_valued && symmetric) + [~, idx] = sort (real (d), "descend"); + else + error ("sigma = \"la\" requires real symmetric problem"); + endif + + case "sa" + if (real_valued && symmetric) + [~, idx] = sort (real (d), "ascend"); + else + error ("sigma = \"sa\" requires real symmetric problem"); + endif + + case "be" + if (real_valued && symmetric) + [~, idx] = sort (real (d), "ascend"); + else + error ("sigma = \"be\" requires real symmetric problem"); + endif + + case "lr" + if (! (real_valued || symmetric)) + [~, idx] = sort (real (d), "descend"); + else + error ("sigma = \"lr\" requires complex or unsymmetric problem"); + endif + + case "sr" + if (! (real_valued || symmetric)) + [~, idx] = sort (real (d), "ascend"); + else + error ("sigma = \"sr\" requires complex or unsymmetric problem"); + endif + + case "li" + if (! (real_valued || symmetric)) + [~, idx] = sort (imag (d), "descend"); + else + error ("sigma = \"li\" requires complex or unsymmetric problem"); + endif + + case "si" + if (! (real_valued || symmetric)) + [~, idx] = sort (imag (d), "ascend"); + else + error ("sigma = \"si\" requires complex or unsymmetric problem"); + endif + + otherwise + error ("unrecognized value for sigma: %s", sigma); + endswitch + endif + + d = d(idx); + + n = numel (d); + + k = min (k, n); + + if (strcmp (sigma, 'be')) + tmp = k / 2; + n1 = floor (tmp); + n2 = n - ceil (tmp) + 1; + selection = [1:floor(k/2), n2:n]; + else + selection = 1:k; + endif + + d = d(selection); + + if (numel (args) == 1) + out{1} = d; + else + out{2} = diag (d); + + v = args{1}; + v = v(:,idx); + out{1} = v(selection,:); + endif + +endfunction + +#### SPARSE MATRIX VERSIONS #### + +## Real positive definite tests, n must be even +%!shared n, k, A, d0, d2 +%! n = 20; +%! k = 4; +%! A = sparse ([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),4*ones(1,n),ones(1,n-2)]); +%! d0 = eig (A); +%! d2 = sort (d0); +%! [~, idx] = sort (abs (d0)); +%! d0 = d0(idx); +%! rand ("state", 42); # initialize generator to make eigs behavior reproducible +%!testif HAVE_ARPACK +%! d1 = eigs (A, k); +%! assert (d1, d0(end:-1:(end-k+1)), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k+1); +%! assert (d1, d0(end:-1:(end-k)), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "lm"); +%! assert (d1, d0(end:-1:(end-k+1)), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! d1 = eigs (A, k, "sm"); +%! assert (d1, d0(k:-1:1), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "la"); +%! assert (d1, d2(end:-1:(end-k+1)), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "sa"); +%! assert (d1, d2(1:k), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "be"); +%! assert (d1, d2([1:floor(k/2), (end - ceil(k/2) + 1):end]), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k+1, "be"); +%! assert (d1, d2([1:floor((k+1)/2), (end - ceil((k+1)/2) + 1):end]), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! d1 = eigs (A, k, 4.1); +%! [~, idx0] = sort (abs (d0 - 4.1)); +%! [~, idx1] = sort (abs (d1 - 4.1)); +%! assert (d1(idx1), d0(idx0(1:k)), 1e-11); +%!testif HAVE_ARPACK, HAVE_CHOLMOD +%! d1 = eigs (A, speye (n), k, "lm"); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! assert (eigs (A, k, 4.1), eigs (A, speye (n), k, 4.1), 1e-11); +%!testif HAVE_ARPACK +%! opts.cholB = true; +%! d1 = eigs (A, speye (n), k, "lm", opts); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! opts.cholB = true; +%! q = [2:n,1]; +%! opts.permB = q; +%! d1 = eigs (A, speye (n)(q,q), k, "lm", opts); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! opts.cholB = true; +%! d1 = eigs (A, speye (n), k, 4.1, opts); +%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! opts.cholB = true; +%! q = [2:n,1]; +%! opts.permB = q; +%! d1 = eigs (A, speye (n)(q,q), k, 4.1, opts); +%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! assert (eigs (A, k, 4.1), eigs (A, speye (n), k, 4.1), 1e-11); +%!testif HAVE_ARPACK +%! fn = @(x) A * x; +%! opts.issym = 1; opts.isreal = 1; +%! d1 = eigs (fn, n, k, "lm", opts); +%! assert (d1, d0(end:-1:(end-k+1)), 1e-11); +%!testif HAVE_ARPACK +%! fn = @(x) A \ x; +%! opts.issym = 1; opts.isreal = 1; +%! d1 = eigs (fn, n, k, "sm", opts); +%! assert (d1, d0(k:-1:1), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! fn = @(x) (A - 4.1 * eye (n)) \ x; +%! opts.issym = 1; opts.isreal = 1; +%! d1 = eigs (fn, n, k, 4.1, opts); +%! assert (d1, eigs (A, k, 4.1), 1e-11); +%!testif HAVE_ARPACK +%! AA = speye (10); +%! fn = @(x) AA * x; +%! opts.issym = 1; opts.isreal = 1; +%! assert (eigs (fn, 10, AA, 3, "lm", opts), [1; 1; 1], 10*eps); +%!testif HAVE_ARPACK +%! [v1,d1] = eigs (A, k, "lm"); +%! d1 = diag (d1); +%! for i=1:k +%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); +%! endfor +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! [v1,d1] = eigs (A, k, "sm"); +%! d1 = diag (d1); +%! for i=1:k +%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); +%! endfor +%!testif HAVE_ARPACK +%! [v1,d1] = eigs (A, k, "la"); +%! d1 = diag (d1); +%! for i=1:k +%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); +%! endfor +%!testif HAVE_ARPACK +%! [v1,d1] = eigs (A, k, "sa"); +%! d1 = diag (d1); +%! for i=1:k +%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); +%! endfor +%!testif HAVE_ARPACK +%! [v1,d1] = eigs (A, k, "be"); +%! d1 = diag (d1); +%! for i=1:k +%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); +%! endfor + +## Real unsymmetric tests +%!shared n, k, A, d0 +%! n = 20; +%! k = 4; +%! A = sparse ([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),1:n,-ones(1,n-2)]); +%! d0 = eig (A); +%! [~, idx] = sort (abs (d0)); +%! d0 = d0(idx); +%! rand ("state", 42); % initialize generator to make eigs behavior reproducible +%!testif HAVE_ARPACK +%! d1 = eigs (A, k); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k+1); +%! assert (abs (d1), abs (d0(end:-1:(end-k))),1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "lm"); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! d1 = eigs (A, k, "sm"); +%! assert (abs (d1), abs (d0(1:k)), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "lr"); +%! [~, idx] = sort (real (d0)); +%! d2 = d0(idx); +%! assert (real (d1), real (d2(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "sr"); +%! [~, idx] = sort (real (abs (d0))); +%! d2 = d0(idx); +%! assert (real (d1), real (d2(1:k)), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "li"); +%! [~, idx] = sort (imag (abs (d0))); +%! d2 = d0(idx); +%! assert (sort (imag (d1)), sort (imag (d2(end:-1:(end-k+1)))), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "si"); +%! [~, idx] = sort (imag (abs (d0))); +%! d2 = d0(idx); +%! assert (sort (imag (d1)), sort (imag (d2(1:k))), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! d1 = eigs (A, k, 4.1); +%! [~, idx0] = sort (abs (d0 - 4.1)); +%! [~, idx1] = sort (abs (d1 - 4.1)); +%! assert (abs (d1(idx1)), abs (d0(idx0(1:k))), 1e-11); +%! assert (sort (imag (d1(idx1))), sort (imag (d0(idx0(1:k)))), 1e-11); +%!testif HAVE_ARPACK, HAVE_CHOLMOD +%! d1 = eigs (A, speye (n), k, "lm"); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! opts.cholB = true; +%! d1 = eigs (A, speye (n), k, "lm", opts); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! opts.cholB = true; +%! q = [2:n,1]; +%! opts.permB = q; +%! d1 = eigs (A, speye (n)(q,q), k, "lm", opts); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! opts.cholB = true; +%! d1 = eigs (A, speye (n), k, 4.1, opts); +%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! opts.cholB = true; +%! q = [2:n,1]; +%! opts.permB = q; +%! d1 = eigs (A, speye (n)(q,q), k, 4.1, opts); +%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! assert (abs (eigs (A, k, 4.1)), abs (eigs (A, speye (n), k, 4.1)), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! assert (sort (imag (eigs (A, k, 4.1))), sort (imag (eigs (A, speye (n), k, 4.1))), 1e-11); +%!testif HAVE_ARPACK +%! fn = @(x) A * x; +%! opts.issym = 0; opts.isreal = 1; +%! d1 = eigs (fn, n, k, "lm", opts); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! fn = @(x) A \ x; +%! opts.issym = 0; opts.isreal = 1; +%! d1 = eigs (fn, n, k, "sm", opts); +%! assert (abs (d1), d0(1:k), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! fn = @(x) (A - 4.1 * eye (n)) \ x; +%! opts.issym = 0; opts.isreal = 1; +%! d1 = eigs (fn, n, k, 4.1, opts); +%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); +%!testif HAVE_ARPACK +%! [v1,d1] = eigs (A, k, "lm"); +%! d1 = diag (d1); +%! for i=1:k +%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); +%! endfor +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! [v1,d1] = eigs (A, k, "sm"); +%! d1 = diag (d1); +%! for i=1:k +%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); +%! endfor +%!testif HAVE_ARPACK +%! [v1,d1] = eigs (A, k, "lr"); +%! d1 = diag (d1); +%! for i=1:k +%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); +%! endfor +%!testif HAVE_ARPACK +%! [v1,d1] = eigs (A, k, "sr"); +%! d1 = diag (d1); +%! for i=1:k +%! assert (max (abs ((A - d1(i)*speye (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)*speye (n))*v1(:,i))), 0, 1e-11); +%! endfor +%!testif HAVE_ARPACK +%! [v1,d1] = eigs (A, k, "si"); +%! d1 = diag (d1); +%! for i=1:k +%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); +%! endfor + + +## Complex hermitian tests +%!shared n, k, A, d0 +%! n = 20; +%! k = 4; +%! A = sparse ([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[1i*ones(1,n-2),4*ones(1,n),-1i*ones(1,n-2)]); +%! d0 = eig (A); +%! [~, idx] = sort (abs (d0)); +%! d0 = d0(idx); +%! rand ("state", 42); % initialize generator to make eigs behavior reproducible +%!testif HAVE_ARPACK +%! d1 = eigs (A, k); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k+1); +%! assert (abs (d1), abs (d0(end:-1:(end-k))),1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "lm"); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! d1 = eigs (A, k, "sm"); +%! assert (abs (d1), abs (d0(1:k)), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "lr"); +%! [~, idx] = sort (real (abs (d0))); +%! d2 = d0(idx); +%! assert (real (d1), real (d2(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "sr"); +%! [~, idx] = sort (real (abs (d0))); +%! d2 = d0(idx); +%! assert (real (d1), real (d2(1:k)), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "li"); +%! [~, idx] = sort (imag (abs (d0))); +%! d2 = d0(idx); +%! assert (sort (imag (d1)), sort (imag (d2(end:-1:(end-k+1)))), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "si"); +%! [~, idx] = sort (imag (abs (d0))); +%! d2 = d0(idx); +%! assert (sort (imag (d1)), sort (imag (d2(1:k))), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! d1 = eigs (A, k, 4.1); +%! [~, idx0] = sort (abs (d0 - 4.1)); +%! [~, idx1] = sort (abs (d1 - 4.1)); +%! assert (abs (d1(idx1)), abs (d0(idx0(1:k))), 1e-11); +%! assert (sort (imag (d1(idx1))), sort (imag (d0(idx0(1:k)))), 1e-11); +%!testif HAVE_ARPACK, HAVE_CHOLMOD +%! d1 = eigs (A, speye (n), k, "lm"); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! opts.cholB = true; +%! d1 = eigs (A, speye (n), k, "lm", opts); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! opts.cholB = true; +%! q = [2:n,1]; +%! opts.permB = q; +%! d1 = eigs (A, speye (n)(q,q), k, "lm", opts); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! opts.cholB = true; +%! d1 = eigs (A, speye (n), k, 4.1, opts); +%! assert (abs (abs (d1)), abs (eigs (A, k, 4.1)), 1e-11); +%! assert (sort (imag (abs (d1))), sort (imag (eigs (A, k, 4.1))), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! opts.cholB = true; +%! q = [2:n,1]; +%! opts.permB = q; +%! d1 = eigs (A, speye (n)(q,q), k, 4.1, opts); +%! assert (abs (abs (d1)), abs (eigs (A, k, 4.1)), 1e-11); +%! assert (sort (imag (abs (d1))), sort (imag (eigs (A, k, 4.1))), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! assert (abs (eigs (A, k, 4.1)), abs (eigs (A, speye (n), k, 4.1)), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! assert (sort (imag (eigs (A, k, 4.1))), sort (imag (eigs (A, speye (n), k, 4.1))), 1e-11); +%!testif HAVE_ARPACK +%! fn = @(x) A * x; +%! opts.issym = 0; opts.isreal = 0; +%! d1 = eigs (fn, n, k, "lm", opts); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! fn = @(x) A \ x; +%! opts.issym = 0; opts.isreal = 0; +%! d1 = eigs (fn, n, k, "sm", opts); +%! assert (abs (d1), d0(1:k), 1e-11); +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! fn = @(x) (A - 4.1 * eye (n)) \ x; +%! opts.issym = 0; opts.isreal = 0; +%! d1 = eigs (fn, n, k, 4.1, opts); +%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); +%!testif HAVE_ARPACK +%! [v1,d1] = eigs (A, k, "lm"); +%! d1 = diag (d1); +%! for i=1:k +%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); +%! endfor +%!testif HAVE_ARPACK, HAVE_UMFPACK +%! [v1,d1] = eigs (A, k, "sm"); +%! d1 = diag (d1); +%! for i=1:k +%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); +%! endfor +%!testif HAVE_ARPACK +%! [v1,d1] = eigs (A, k, "lr"); +%! d1 = diag (d1); +%! for i=1:k +%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); +%! endfor +%!testif HAVE_ARPACK +%! [v1,d1] = eigs (A, k, "sr"); +%! d1 = diag (d1); +%! for i=1:k +%! assert (max (abs ((A - d1(i)*speye (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)*speye (n))*v1(:,i))), 0, 1e-11); +%! endfor +%!testif HAVE_ARPACK +%! [v1,d1] = eigs (A, k, "si"); +%! d1 = diag (d1); +%! for i=1:k +%! assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11); +%! endfor + +#### FULL MATRIX VERSIONS #### + +## Real positive definite tests, n must be even +%!shared n, k, A, d0, d2 +%! n = 20; +%! k = 4; +%! A = full (sparse ([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),4*ones(1,n),ones(1,n-2)])); +%! d0 = eig (A); +%! d2 = sort (d0); +%! [~, idx] = sort (abs (d0)); +%! d0 = d0(idx); +%! rand ("state", 42); % initialize generator to make eigs behavior reproducible +%!testif HAVE_ARPACK +%! d1 = eigs (A, k); +%! assert (d1, d0(end:-1:(end-k+1)), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k+1); +%! assert (d1, d0(end:-1:(end-k)),1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "lm"); +%! assert (d1, d0(end:-1:(end-k+1)), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "sm"); +%! assert (d1, d0(k:-1:1), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "la"); +%! assert (d1, d2(end:-1:(end-k+1)), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "sa"); +%! assert (d1, d2(1:k), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "be"); +%! assert (d1, d2([1:floor(k/2), (end - ceil(k/2) + 1):end]), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k+1, "be"); +%! assert (d1, d2([1:floor((k+1)/2), (end - ceil((k+1)/2) + 1):end]), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, 4.1); +%! [~, idx0] = sort (abs (d0 - 4.1)); +%! [~, idx1] = sort (abs (d1 - 4.1)); +%! assert (d1(idx1), d0(idx0(1:k)), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, eye (n), k, "lm"); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! assert (eigs (A, k, 4.1), eigs (A, eye (n), k, 4.1), 1e-11); +%!testif HAVE_ARPACK +%! opts.cholB = true; +%! d1 = eigs (A, eye (n), k, "lm", opts); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! opts.cholB = true; +%! q = [2:n,1]; +%! opts.permB = q; +%! d1 = eigs (A, eye (n)(q,q), k, "lm", opts); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! opts.cholB = true; +%! d1 = eigs (A, eye (n), k, 4.1, opts); +%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); +%!testif HAVE_ARPACK +%! opts.cholB = true; +%! q = [2:n,1]; +%! opts.permB = q; +%! d1 = eigs (A, eye (n)(q,q), k, 4.1, opts); +%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); +%!testif HAVE_ARPACK +%! assert (eigs (A, k, 4.1), eigs (A, eye (n), k, 4.1), 1e-11); +%!testif HAVE_ARPACK +%! fn = @(x) A * x; +%! opts.issym = 1; opts.isreal = 1; +%! d1 = eigs (fn, n, k, "lm", opts); +%! assert (d1, d0(end:-1:(end-k+1)), 1e-11); +%!testif HAVE_ARPACK +%! fn = @(x) A \ x; +%! opts.issym = 1; opts.isreal = 1; +%! d1 = eigs (fn, n, k, "sm", opts); +%! assert (d1, d0(k:-1:1), 1e-11); +%!testif HAVE_ARPACK +%! fn = @(x) (A - 4.1 * eye (n)) \ x; +%! opts.issym = 1; opts.isreal = 1; +%! d1 = eigs (fn, n, k, 4.1, opts); +%! assert (d1, eigs (A, k, 4.1), 1e-11); +%!testif HAVE_ARPACK +%! [v1,d1] = eigs (A, k, "lm"); +%! d1 = diag (d1); +%! 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, "sm"); +%! d1 = diag (d1); +%! 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, "la"); +%! d1 = diag (d1); +%! 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, "sa"); +%! d1 = diag (d1); +%! 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, "be"); +%! d1 = diag (d1); +%! for i=1:k +%! assert (max (abs ((A - d1(i)*eye (n))*v1(:,i))), 0, 1e-11); +%! endfor + +## Real unsymmetric tests +%!shared n, k, A, d0 +%! n = 20; +%! k = 4; +%! A = full (sparse ([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),1:n,-ones(1,n-2)])); +%! d0 = eig (A); +%! [~, idx] = sort (abs (d0)); +%! d0 = d0(idx); +%! rand ("state", 42); % initialize generator to make eigs behavior reproducible +%!testif HAVE_ARPACK +%! d1 = eigs (A, k); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k+1); +%! assert (abs (d1), abs (d0(end:-1:(end-k))),1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "lm"); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "sm"); +%! assert (abs (d1), abs (d0(1:k)), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "lr"); +%! [~, idx] = sort (real (d0)); +%! d2 = d0(idx); +%! assert (real (d1), real (d2(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "sr"); +%! [~, idx] = sort (real (abs (d0))); +%! d2 = d0(idx); +%! assert (real (d1), real (d2(1:k)), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "li"); +%! [~, idx] = sort (imag (abs (d0))); +%! d2 = d0(idx); +%! assert (sort (imag (d1)), sort (imag (d2(end:-1:(end-k+1)))), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "si"); +%! [~, idx] = sort (imag (abs (d0))); +%! d2 = d0(idx); +%! assert (sort (imag (d1)), sort (imag (d2(1:k))), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, 4.1); +%! [~, idx0] = sort (abs (d0 - 4.1)); +%! [~, idx1] = sort (abs (d1 - 4.1)); +%! assert (abs (d1(idx1)), abs (d0(idx0(1:k))), 1e-11); +%! assert (sort (imag (d1(idx1))), sort (imag (d0(idx0(1:k)))), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, eye (n), k, "lm"); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! opts.cholB = true; +%! d1 = eigs (A, eye (n), k, "lm", opts); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! opts.cholB = true; +%! q = [2:n,1]; +%! opts.permB = q; +%! d1 = eigs (A, eye (n)(q,q), k, "lm", opts); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! opts.cholB = true; +%! d1 = eigs (A, eye (n), k, 4.1, opts); +%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); +%!testif HAVE_ARPACK +%! opts.cholB = true; +%! q = [2:n,1]; +%! opts.permB = q; +%! d1 = eigs (A, eye (n)(q,q), k, 4.1, opts); +%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); +%!testif HAVE_ARPACK +%! assert (abs (eigs (A, k, 4.1)), abs (eigs (A, eye (n), k, 4.1)), 1e-11); +%!testif HAVE_ARPACK +%! assert (sort (imag (eigs (A, k, 4.1))), sort (imag (eigs (A, eye (n), k, 4.1))), 1e-11); +%!testif HAVE_ARPACK +%! fn = @(x) A * x; +%! opts.issym = 0; opts.isreal = 1; +%! d1 = eigs (fn, n, k, "lm", opts); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! fn = @(x) A \ x; +%! opts.issym = 0; opts.isreal = 1; +%! d1 = eigs (fn, n, k, "sm", opts); +%! assert (abs (d1), d0(1:k), 1e-11); +%!testif HAVE_ARPACK +%! fn = @(x) (A - 4.1 * eye (n)) \ x; +%! opts.issym = 0; opts.isreal = 1; +%! d1 = eigs (fn, n, k, 4.1, opts); +%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); +%!testif HAVE_ARPACK +%! [v1,d1] = eigs (A, k, "lm"); +%! d1 = diag (d1); +%! 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, "sm"); +%! d1 = diag (d1); +%! 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, "lr"); +%! d1 = diag (d1); +%! 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, "sr"); +%! d1 = diag (d1); +%! 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 +%!testif HAVE_ARPACK +%! [v1,d1] = eigs (A, k, "si"); +%! d1 = diag (d1); +%! for i=1:k +%! assert (max (abs ((A - d1(i)*eye (n))*v1(:,i))), 0, 1e-11); +%! endfor + +## Complex hermitian tests +%!shared n, k, A, d0 +%! n = 20; +%! k = 4; +%! A = full (sparse ([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[1i*ones(1,n-2),4*ones(1,n),-1i*ones(1,n-2)])); +%! d0 = eig (A); +%! [~, idx] = sort (abs (d0)); +%! d0 = d0(idx); +%! rand ("state", 42); % initialize generator to make eigs behavior reproducible +%!testif HAVE_ARPACK +%! d1 = eigs (A, k); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k+1); +%! assert (abs (d1), abs (d0(end:-1:(end-k))),1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "lm"); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "sm"); +%! assert (abs (d1), abs (d0(1:k)), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "lr"); +%! [~, idx] = sort (real (abs (d0))); +%! d2 = d0(idx); +%! assert (real (d1), real (d2(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "sr"); +%! [~, idx] = sort (real (abs (d0))); +%! d2 = d0(idx); +%! assert (real (d1), real (d2(1:k)), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "li"); +%! [~, idx] = sort (imag (abs (d0))); +%! d2 = d0(idx); +%! assert (sort (imag (d1)), sort (imag (d2(end:-1:(end-k+1)))), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, "si"); +%! [~, idx] = sort (imag (abs (d0))); +%! d2 = d0(idx); +%! assert (sort (imag (d1)), sort (imag (d2(1:k))), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, k, 4.1); +%! [~, idx0] = sort (abs (d0 - 4.1)); +%! [~, idx1] = sort (abs (d1 - 4.1)); +%! assert (abs (d1(idx1)), abs (d0(idx0(1:k))), 1e-11); +%! assert (sort (imag (d1(idx1))), sort (imag (d0(idx0(1:k)))), 1e-11); +%!testif HAVE_ARPACK +%! d1 = eigs (A, eye (n), k, "lm"); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! opts.cholB = true; +%! d1 = eigs (A, eye (n), k, "lm", opts); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! opts.cholB = true; +%! q = [2:n,1]; +%! opts.permB = q; +%! d1 = eigs (A, eye (n)(q,q), k, "lm", opts); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! opts.cholB = true; +%! d1 = eigs (A, eye (n), k, 4.1, opts); +%! assert (abs (abs (d1)), abs (eigs (A, k, 4.1)), 1e-11); +%! assert (sort (imag (abs (d1))), sort (imag (eigs (A, k, 4.1))), 1e-11); +%!testif HAVE_ARPACK +%! opts.cholB = true; +%! q = [2:n,1]; +%! opts.permB = q; +%! d1 = eigs (A, eye (n)(q,q), k, 4.1, opts); +%! assert (abs (abs (d1)), abs (eigs (A, k, 4.1)), 1e-11); +%! assert (sort (imag (abs (d1))), sort (imag (eigs (A, k, 4.1))), 1e-11); +%!testif HAVE_ARPACK +%! assert (abs (eigs (A, k, 4.1)), abs (eigs (A, eye (n), k, 4.1)), 1e-11); +%!testif HAVE_ARPACK +%! assert (sort (imag (eigs (A, k, 4.1))), sort (imag (eigs (A, eye (n), k, 4.1))), 1e-11); +%!testif HAVE_ARPACK +%! fn = @(x) A * x; +%! opts.issym = 0; opts.isreal = 0; +%! d1 = eigs (fn, n, k, "lm", opts); +%! assert (abs (d1), abs (d0(end:-1:(end-k+1))), 1e-11); +%!testif HAVE_ARPACK +%! fn = @(x) A \ x; +%! opts.issym = 0; opts.isreal = 0; +%! d1 = eigs (fn, n, k, "sm", opts); +%! assert (abs (d1), d0(1:k), 1e-11); +%!testif HAVE_ARPACK +%! fn = @(x) (A - 4.1 * eye (n)) \ x; +%! opts.issym = 0; opts.isreal = 0; +%! d1 = eigs (fn, n, k, 4.1, opts); +%! assert (abs (d1), eigs (A, k, 4.1), 1e-11); +%!testif HAVE_ARPACK +%! [v1,d1] = eigs (A, k, "lm"); +%! d1 = diag (d1); +%! 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, "sm"); +%! d1 = diag (d1); +%! 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, "lr"); +%! d1 = diag (d1); +%! 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, "sr"); +%! d1 = diag (d1); +%! 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 +%!testif HAVE_ARPACK +%! [v1,d1] = eigs (A, k, "si"); +%! d1 = diag (d1); +%! for i=1:k +%! assert (max (abs ((A - d1(i)*eye (n))*v1(:,i))), 0, 1e-11); +%! endfor + +%!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]);