Mercurial > hg > octave-lyh

diff scripts/sparse/svds.m @ 16414:610e02bf9579

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
use the default ARPACK convergence tolerance (bug #34461) * svds.m: do not divide the input tol by max abs of input matrix, set output flag to return from eigs().
author Ed Meyer <eem2314@gmail.com>
date Wed, 08 Aug 2012 20:51:24 -0700
parents 94d21131fefd
children c8bbab6b9e7a
line wrap: on
line diff
--- a/scripts/sparse/svds.m
+++ b/scripts/sparse/svds.m
@@ -103,7 +103,7 @@
   endif
 
   if (nargin < 4)
-    opts.tol = 1e-10 / root2;
+    opts.tol = 0;   ## use ARPACK default
     opts.disp = 0;
     opts.maxit = 300;
   else
@@ -111,7 +111,7 @@
       error ("svds: OPTS must be a structure");
     endif
     if (!isfield (opts, "tol"))
-      opts.tol = 1e-10 / root2;
+      opts.tol = 0;   ## use ARPACK default
     else
       opts.tol = opts.tol / root2;
     endif
@@ -152,7 +152,6 @@
     b_opts = opts;
     ## Call to eigs is always a symmetric matrix by construction
     b_opts.issym = true;
-    b_opts.tol = opts.tol / max_a;
     b_sigma = sigma;
     if (!ischar (b_sigma))
       b_sigma = b_sigma / max_a;
@@ -235,14 +234,14 @@
     endif
 
     if (nargout > 3)
-      flag = norm (A*v - u*s, 1) > root2 * opts.tol * norm (A, 1);
+      flag = (flag != 0);
     endif
   endif
 
 endfunction
 
 
-%!shared n, k, A, u, s, v, opts, rand_state, randn_state
+%!shared n, k, A, u, s, v, opts, rand_state, randn_state, tol
 %! n = 100;
 %! k = 7;
 %! A = sparse ([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),0.4*n*ones(1,n),ones(1,n-2)]);
@@ -263,13 +262,15 @@
 %! [u2,s2,v2,flag] = svds (A,k);
 %! s2 = diag (s2);
 %! assert (flag, !1);
-%! assert (s2, s(end:-1:end-k+1), 1e-10);
+%! tol = 10 * eps() * norm(s2, 1);
+%! assert (s2, s(end:-1:end-k+1), tol);
 %!
 %!testif HAVE_ARPACK, HAVE_UMFPACK
 %! [u2,s2,v2,flag] = svds (A,k,0,opts);
 %! s2 = diag (s2);
 %! assert (flag, !1);
-%! assert (s2, s(k:-1:1), 1e-10);
+%! tol = 10 * eps() * norm(s2, 1);
+%! assert (s2, s(k:-1:1), tol);
 %!
 %!testif HAVE_ARPACK, HAVE_UMFPACK
 %! idx = floor (n/2);
@@ -278,7 +279,8 @@
 %! [u2,s2,v2,flag] = svds (A,k,sigma,opts);
 %! s2 = diag (s2);
 %! assert (flag, !1);
-%! assert (s2, s((idx+floor(k/2)):-1:(idx-floor(k/2))), 1e-10);
+%! tol = 10 * eps() * norm(s2, 1);
+%! assert (s2, s((idx+floor(k/2)):-1:(idx-floor(k/2))), tol);
 %!
 %!testif HAVE_ARPACK
 %! [u2,s2,v2,flag] = svds (zeros (10), k);