changeset 27661:2213e82bb9d3

svds should *always* return real singular values (bug #57185). * svds.m: Use real() on output of eigs() to guarantee real-valued results. Add BIST test for bug #57185.
author Rik <rik@octave.org>
date Sat, 09 Nov 2019 11:33:55 -0800
parents 30e84a3d58e5
children 3b078b750181
files scripts/sparse/svds.m
diffstat 1 files changed, 13 insertions(+), 6 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/sparse/svds.m	Sat Nov 09 12:45:07 2019 +0100
+++ b/scripts/sparse/svds.m	Sat Nov 09 11:33:55 2019 -0800
@@ -179,12 +179,14 @@
       b_k = k;  # Normal case, find just the k largest eigenvalues
     endif
 
+    ## FIXME: singular values are always real.  Until eigs() guarantees
+    ## real results for symmetric inputs (bug #57196), use real() here.
     if (nargout > 1)
       [V, s, flag] = eigs ([sparse(m,m), b; b', sparse(n,n)],
                            b_k, b_sigma, b_opts);
-      s = diag (s);
+      s = real (diag (s));
     else
-      s = eigs ([sparse(m,m), b; b', sparse(n,n)], b_k, b_sigma, b_opts);
+      s = real (eigs ([sparse(m,m), b; b', sparse(n,n)], b_k, b_sigma, b_opts));
     endif
 
     if (ischar (sigma))
@@ -271,21 +273,21 @@
 %! rand ("state", 42);
 %! opts.v0 = rand (2*n,1);  # Initialize eigs ARPACK starting vector
 %!                          # to guarantee reproducible results
-%!
+
 %!testif HAVE_ARPACK
 %! [u2,s2,v2,flag] = svds (A,k);
 %! s2 = diag (s2);
 %! assert (flag, ! 1);
 %! 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);
 %! tol = 10 * eps() * norm(s2, 1);
 %! assert (s2, s(k:-1:1), tol);
-%!
+
 %!testif HAVE_ARPACK, HAVE_UMFPACK
 %! idx = floor (n/2);
 %! % Don't put sigma right on a singular value or there are convergence issues
@@ -295,7 +297,7 @@
 %! assert (flag, ! 1);
 %! 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);
 %! assert (u2, eye (10, k));
@@ -306,6 +308,11 @@
 %! s = svds (speye (10));
 %! assert (s, ones (6, 1), 8*eps);
 
+%!testif HAVE_ARPACK <57185>
+%! z = complex (ones (10), ones (10));
+%! s = svds (z);
+%! assert (isreal (s));
+
 %!test
 %! ## Restore random number generator seed at end of tests
 %! rand ("state", rand_state);