Mercurial > octave

changeset 25302:94b4c19c0fd9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Merge stable to default.
author John W. Eaton <jwe@octave.org>
date Tue, 24 Apr 2018 11:02:55 -0400
parents 5ca8abb27644 (current diff) 1a632692a58e (diff)
children d73dfd34230b
files
diffstat 4 files changed, 128 insertions(+), 80 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/numeric/eigs-base.cc	Mon Apr 23 17:41:21 2018 -0700
+++ b/liboctave/numeric/eigs-base.cc	Tue Apr 24 11:02:55 2018 -0400
@@ -210,7 +210,7 @@
   octave_idx_type info;
   octave::math::sparse_chol<SparseMatrix> fact (b, info, false);
 
-  if (fact.P () != 0)
+  if (info != 0)
     return false;
   else
     {
@@ -248,7 +248,7 @@
   octave_idx_type info;
   octave::math::sparse_chol<SparseComplexMatrix> fact (b, info, false);
 
-  if (fact.P () != 0)
+  if (info != 0)
     return false;
   else
     {
--- a/m4/acinclude.m4	Mon Apr 23 17:41:21 2018 -0700
+++ b/m4/acinclude.m4	Tue Apr 24 11:02:55 2018 -0400
@@ -838,7 +838,11 @@
           warn_$1=
           AC_DEFINE([HAVE_]m4_toupper([$1]), 1,
             [Define to 1 if $2 is available.])], [$8])
+      else
+        m4_toupper([$1])_LIBS=
       fi
+    else
+      m4_toupper([$1])_LIBS=
     fi
     m4_ifnblank([$6], [AC_LANG_POP($6)])
     CPPFLAGS="$ac_octave_save_CPPFLAGS"
@@ -855,7 +859,6 @@
   AC_SUBST(m4_toupper([$1])_LIBS)
   if test -n "$warn_$1"; then
     OCTAVE_CONFIGURE_WARNING([warn_$1])
-    m4_toupper([$1])_LIBS=
   fi
 ])
 dnl
--- a/scripts/sparse/eigs.m	Mon Apr 23 17:41:21 2018 -0700
+++ b/scripts/sparse/eigs.m	Tue Apr 24 11:02:55 2018 -0400
@@ -215,7 +215,7 @@
       offset = 1;
     endif
 
-    if (rows (a) < 9)
+    if (rows (a) < 13)
       call_eig = true;
     endif
 
@@ -224,10 +224,12 @@
       if (isnumeric (tmp) && isscalar (tmp) && isreal (tmp)
           && round (tmp) == tmp)
         k = tmp;
-
-        if (rows (a) - k < 3)
-          call_eig = true;
-        endif
+        p = 2 * k;
+      elseif (isfield (tmp, "p"))
+          p = tmp.p;
+      endif
+      if (p >= rows (a))
+        call_eig = true;
       else
         call_eig = false;
       endif
@@ -238,6 +240,24 @@
           sigma = tolower (tmp);
         elseif (isnumeric (tmp) && isscalar (tmp))
           sigma = tmp;
+        elseif (isfield (tmp, "p"))
+          p = tmp.p;
+        endif
+        if (p >= rows (a))
+          call_eig = true;
+        else
+          call_eig = false;
+        endif
+      endif
+      if (nargin > 3 + offset)
+        tmp = varargin{4+offset};
+        if (isfield (tmp, "p"))
+          p = tmp.p;
+        else
+          p = 2 * k;
+        endif
+        if (p >= rows (a))
+          call_eig = true;
         else
           call_eig = false;
         endif
@@ -305,28 +325,28 @@
         endif
 
       case "lr"
-        if (! (real_valued || symmetric))
+        if (! (real_valued && symmetric))
           [~, idx] = sort (real (d), "descend");
         else
           error ('eigs: SIGMA = "lr" requires complex or unsymmetric problem');
         endif
 
       case "sr"
-        if (! (real_valued || symmetric))
+        if (! (real_valued && symmetric))
           [~, idx] = sort (real (d), "ascend");
         else
           error ('eigs: SIGMA = "sr" requires complex or unsymmetric problem');
         endif
 
       case "li"
-        if (! (real_valued || symmetric))
+        if (! (real_valued && symmetric))
           [~, idx] = sort (imag (d), "descend");
         else
           error ('eigs: SIGMA = "li" requires complex or unsymmetric problem');
         endif
 
       case "si"
-        if (! (real_valued || symmetric))
+        if (! (real_valued && symmetric))
           [~, idx] = sort (imag (d), "ascend");
         else
           error ('eigs: SIGMA = "si" requires complex or unsymmetric problem');
@@ -752,15 +772,15 @@
 %!testif HAVE_ARPACK, HAVE_UMFPACK
 %! A = toeplitz (sparse (1:10));
 %! B = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, 10));
-%! [v, d] = eigs (A, B, 5, "lm");
-%! for i = 1:5
+%! [v, d] = eigs (A, B, 4, "lm");
+%! for i = 1:4
 %!   assert (A * v(:,i), d(i, i) * B * v(:,i), 1e-12)
 %! endfor
 %! ddiag = diag (d);
 %! [ddiag, idx] = sort (ddiag);
 %! v = v(:, idx);
 %! R = chol (B);
-%! [v1, d1] = eigs (R' \ A / R, 5, "lm");
+%! [v1, d1] = eigs (R' \ A / R, 4, "lm");
 %! d1diag = diag (d1);
 %! [d1diag, idx] = sort (d1diag);
 %! v1 = v1(:, idx);
@@ -768,15 +788,15 @@
 %!testif HAVE_ARPACK, HAVE_UMFPACK
 %! A = spdiags ([[1./(2:11)]',[-5:-2:-23]',[1:10]'],-1:1,10,10);
 %! B = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, 10));
-%! [v, d] = eigs (A, B, 5, "lm");
-%! for i = 1:5
+%! [v, d] = eigs (A, B, 4, "lm");
+%! for i = 1:4
 %!   assert (A * v(:,i), d(i, i) * B * v(:,i), 1e-12);
 %! endfor
 %! ddiag = diag (d);
 %! [ddiag, idx] = sort (ddiag);
 %! v = v(:, idx);
 %! R = chol (B);
-%! [v1, d1] = eigs (R' \ A / R, 5, "lm");
+%! [v1, d1] = eigs (R' \ A / R, 4, "lm");
 %! d1diag = diag (d1);
 %! [d1diag, idx] = sort (d1diag);
 %! v1 = v1(:, idx);
@@ -785,15 +805,15 @@
 %! A = spdiags ([[1./(2:11)]',[-5:-2:-23]',[1:10]'],-1:1,10,10) -...
 %! 1i * spdiags ([[1./(2:11)]',[-5:-2:-23]',[1:10]'],-1:1,10,10);
 %! B = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, 10));
-%! [v, d] = eigs (A, B, 5, "lm");
-%! for i = 1:5
+%! [v, d] = eigs (A, B, 4, "lm");
+%! for i = 1:4
 %!   assert (A * v(:,i), d(i, i) * B * v(:,i), 1e-12);
 %! endfor
 %! ddiag = diag (d);
 %! [ddiag, idx] = sort (ddiag);
 %! v = v(:, idx);
 %! R = chol (B);
-%! [v1, d1] = eigs (R' \ A / R, 5, "lm");
+%! [v1, d1] = eigs (R' \ A / R, 4, "lm");
 %! d1diag = diag (d1);
 %! [d1diag, idx] = sort (d1diag);
 %! v1 = v1(:, idx);
@@ -801,15 +821,15 @@
 %!testif HAVE_ARPACK, HAVE_UMFPACK
 %! A = toeplitz (sparse (1:10));
 %! B = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, 10));
-%! [v, d] = eigs (A, B, 5, 1);
-%! for i = 1:5
+%! [v, d] = eigs (A, B, 4, 1);
+%! for i = 1:4
 %!   assert (A * v(:,i), d(i, i) * B * v(:,i), 1e-12);
 %! endfor
 %! ddiag = diag (d);
 %! [ddiag, idx] = sort (ddiag);
 %! v = v(:, idx);
 %! R = chol (B);
-%! [v1, d1] = eigs (R' \ A / R, 5, 1);
+%! [v1, d1] = eigs (R' \ A / R, 4, 1);
 %! d1diag = diag (d1);
 %! [d1diag, idx] = sort (d1diag);
 %! v1 = v1(:, idx);
@@ -817,15 +837,15 @@
 %!testif HAVE_ARPACK, HAVE_UMFPACK
 %! A = spdiags ([[1./(2:11)]',[-5:-2:-23]',[1:10]'],-1:1,10,10);
 %! B = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, 10));
-%! [v, d] = eigs (A, B, 5, 1);
-%! for i = 1:5
+%! [v, d] = eigs (A, B, 4, 1);
+%! for i = 1:4
 %!   assert (A * v(:,i), d(i, i) * B * v(:,i), 1e-12);
 %! endfor
 %! ddiag = diag (d);
 %! [ddiag, idx] = sort (ddiag);
 %! v = v(:, idx);
 %! R = chol (B);
-%! [v1, d1] = eigs (R' \ A / R, 5, 1);
+%! [v1, d1] = eigs (R' \ A / R, 4, 1);
 %! d1diag = diag (d1);
 %! [d1diag, idx] = sort (d1diag);
 %! v1 = v1(:, idx);
@@ -834,15 +854,15 @@
 %! A = spdiags ([[1./(2:11)]',[-5:-2:-23]',[1:10]'],-1:1,10,10) -...
 %! 1i * spdiags ([[1./(2:11)]',[-5:-2:-23]',[1:10]'],-1:1,10,10);
 %! B = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, 10));
-%! [v, d] = eigs (A, B, 5, 1);
-%! for i = 1:5
+%! [v, d] = eigs (A, B, 4, 1);
+%! for i = 1:4
 %!   assert (A * v(:,i), d(i, i) * B * v(:,i), 1e-12);
 %! endfor
 %! ddiag = diag (d);
 %! [ddiag, idx] = sort (ddiag);
 %! v = v(:, idx);
 %! R = chol (B);
-%! [v1, d1] = eigs (R' \ A / R, 5, 1);
+%! [v1, d1] = eigs (R' \ A / R, 4, 1);
 %! d1diag = diag (d1);
 %! [d1diag, idx] = sort (d1diag);
 %! v1 = v1(:, idx);
@@ -852,8 +872,8 @@
 %! B = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, 10));
 %! R = chol (B);
 %! opts.cholB = R;
-%! [v, d] = eigs (A, R, 5, "lm", opts);
-%! for i = 1:5
+%! [v, d] = eigs (A, R, 4, "lm", opts);
+%! for i = 1:4
 %!   assert (A * v(:,i), d(i, i) * B * v(:,i), 1e-12)
 %! endfor
 %!testif HAVE_ARPACK, HAVE_UMFPACK
@@ -862,8 +882,8 @@
 %! [R, ~, permB] = chol (B, "vector");
 %! opts.cholB = R;
 %! opts.permB = permB;
-%! [v, d] = eigs (A, R, 5, "lm", opts);
-%! for i = 1:5
+%! [v, d] = eigs (A, R, 4, "lm", opts);
+%! for i = 1:4
 %!   assert (A * v(:,i), d(i, i) * B * v(:,i), 1e-12)
 %! endfor
 
@@ -1250,10 +1270,10 @@
 
 %!test
 %! A = 2 * diag (ones (10, 1)) - diag (ones (9, 1), 1) - diag (ones (9, 1), -1);
-%! B = diag (ones (10, 1));
+%! B = eye (10);
 %! reseig = eig (A, B);
 %! [~, idx] = sort (abs (reseig), "ascend");
-%! assert (eigs (A, B, 10, 0), reseig (idx));
+%! assert (eigs (A, B, 4, 0), reseig (idx(4:-1:1)), 8 * eps);
 %!testif HAVE_ARPACK
 %! A = eye (9);
 %! A(1, 1) = 0;
@@ -1266,47 +1286,47 @@
 %! A(10,:) = [-1, zeros(1, 8), -1];
 %! opts.v0 = (1:10)';
 %! typ = "lr";
-%! [v, m] = eigs (A, 5, typ, opts);
+%! [v, m] = eigs (A, 4, 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);
+%!         [0.514038; 0.514038; 0.880290; 0.880290], 1e-4);
+%! m = eigs (A, 4, typ, opts);
 %! assert (sort (real (m)), ...
-%!         [-0.081751; 0.514038; 0.514038; 0.880290; 0.880290], 1e-4);
+%!         [0.514038; 0.514038; 0.880290; 0.880290], 1e-4);
 %! typ = "li";
-%! [v, m] = eigs (A, 5, typ, opts);
+%! [v, m] = eigs (A, 4, 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);
+%!         [0.78972; 0.78972; 0.96518; 0.96518], 1e-4);
+%! m = eigs (A, 4, typ, opts);
 %! assert (sort (abs (imag (m))), ...
-%!         [0.75447; 0.78972; 0.78972; 0.96518; 0.96518], 1e-4);
+%!         [0.78972; 0.78972; 0.96518; 0.96518], 1e-4);
 %! typ = "sr";
-%! [v, m] = eigs (A, 5, typ, opts);
+%! [v, m] = eigs (A, 4, 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);
+%!         [-1.12180; -1.12180; -0.69077; -0.69077], 1e-4);
+%! m = eigs (A, 4, typ, opts);
 %! assert (sort (real (m)), ...
-%!         [-1.12180; -1.12180; -0.69077; -0.69077; -0.08175], 1e-4);
+%!         [-1.12180; -1.12180; -0.69077; -0.69077], 1e-4);
 %! typ = "si";
-%! [v, m] = eigs (A, 5, typ, opts);
+%! [v, m] = eigs (A, 4, 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);
+%!         [0.25552; 0.25552; 0.30282; 0.30282], 1e-4);
+%! m = eigs (A, 4, typ, opts);
 %! assert (sort (abs (imag (m))), ...
-%!         [0.25552; 0.25552; 0.30282; 0.30282; 0.75447], 1e-4);
+%!         [0.25552; 0.25552; 0.30282; 0.30282], 1e-4);
 %! typ = "lm";
-%! [v, m] = eigs (A, 5, typ, opts);
+%! [v, m] = eigs (A, 4, 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);
+%!         [1.02294;  1.02294; 1.15054; 1.15054], 1e-4);
+%! m = eigs (A, 4, typ, opts);
 %! assert (sort (abs (m)), ...
-%!         [0.96863; 1.02294; 1.02294; 1.15054; 1.15054], 1e-4);
+%!         [1.02294; 1.02294; 1.15054; 1.15054], 1e-4);
 %! typ = "sm";
-%! [v, m] = eigs (A, 5, typ, opts);
+%! [v, m] = eigs (A, 4, 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);
+%!         [0.93092; 0.93092; 0.94228; 0.94228], 1e-4);
+%! m = eigs (A, 4, typ, opts);
 %! assert (sort (abs (m)), ...
-%!         [0.93092; 0.93092; 0.94228; 0.94228; 0.96863], 1e-4);
+%!         [0.93092; 0.93092; 0.94228; 0.94228], 1e-4);
 %!testif HAVE_ARPACK
 %! A = toeplitz (sparse ([2, 1, zeros(1,8)]));
 %! opts.v0 = (1:10)';
@@ -1329,14 +1349,14 @@
 %! Z = X * X';
 %! r = rank (Z);
 %! assert (r, 8);
-%! [V, D] = eigs (Z, r, "lm");
+%! [V, D] = eigs (Z, r, "lm"); # call_eig is true
 %! ZZ = V * D * V';
 %! tmp = abs (Z - ZZ);
 %! assert (max (tmp(:)) < 5e-11);
 
-%!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])
+%!assert (eigs (diag (1:5), 5, "sa"), [1;2;3;4;5]) # call_eig is true
+%!assert (eigs (diag (1:5), 5, "la"), [5;4;3;2;1]) # call_eig is true
+%!assert (eigs (diag (1:5), 3, "be"), [1;4;5]) # call_eig is true
 %!error
 %! A = rand (10);
 %! opts.v0 = ones (8, 1);
@@ -1443,3 +1463,24 @@
 %! warning ("off", "Octave:eigs:UnconvergedEigenvalues", "local");
 %! d = eigs (Afun, 100, 6, "lm", opts);
 %! assert (d(6), NaN+1i*NaN);
+%!testif HAVE_ARPACK
+%! A = sparse (magic (10));
+%! B = sparse (magic (10)); # not HPD
+%! fail ("eigs (A, B, 4)", "eigs: The matrix B is not positive definite")
+%!testif HAVE_ARPACK
+%! i_A = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10];
+%! j_A = [1, 2, 3, 4, 5, 6, 7,  8, 9, 10];
+%! v_A = [1, 2i, 3, 4i, 5, 6i, 7, 8, 9, 10i];
+%! A = sparse(i_A, j_A, v_A);
+%! i_B = [1,2, 3, 4, 5, 6, 7, 8, 9, 10];
+%! j_B = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10];
+%! v_B = [3, 10i, 1, 8i, 7, 6i, 5, 4i, 9, 7i];
+%! B = sparse(i_B, j_B, v_B); # not SPD
+%! [Evectors Evalues] = eigs(A, B, 5, 'SM'); # call_eig is true
+%! ResidualVectors = A * Evectors - B * Evectors * Evalues;
+%! RelativeErrors = norm (ResidualVectors, "columns") ./ ...
+%! norm (A * Evectors, "columns");
+%! assert (RelativeErrors, zeros (1, 5))
+%!testif HAVE_ARPACK
+%! A = rand (8);
+%! eigs (A, 6, "lr"); # this failed on 4.2.x
--- a/scripts/testfun/speed.m	Mon Apr 23 17:41:21 2018 -0700
+++ b/scripts/testfun/speed.m	Tue Apr 24 11:02:55 2018 -0400
@@ -142,9 +142,10 @@
 ##
 ## Type @kbd{example ("speed")} to see some real examples or
 ## @kbd{demo ("speed")} to run them.
+##
 ## @end deftypefn
 
-## Programming Note: All variables for speed() must use the internal prefix "__".
+## Programming Note: All variables for speed must use the internal prefix "__".
 ## Shared variables are eval'ed into the current workspace and therefore might
 ## collide with the names used in the speed.m function itself.
 
@@ -219,19 +220,19 @@
       fflush (stdout);
     endif
 
-    eval (["__t = time();" __f1 "__v1=ans; __t = time()-__t;"]);
+    eval (["__tid = tic();" __f1 "__v1=ans; __t = toc(__tid);"]);
     if (__t < 0.25)
-      eval (["__t2 = time();" __f1 "__t2 = time()-__t2;"]);
-      eval (["__t3 = time();" __f1 "__t3 = time()-__t3;"]);
+      eval (["__tid = tic();" __f1 "__t2 = toc(__tid);"]);
+      eval (["__tid = tic();" __f1 "__t3 = toc(__tid);"]);
       __t = min ([__t, __t2, __t3]);
     endif
     __tnew(k) = __t;
 
     if (! isempty (__f2))
-      eval (["__t = time();" __f2 "__v2=ans; __t = time()-__t;"]);
+      eval (["__tid = tic();" __f2 "__v2=ans; __t = toc(__tid);"]);
       if (__t < 0.25)
-        eval (["__t2 = time();" __f2 "__t2 = time()-__t2;"]);
-        eval (["__t3 = time();" __f2 "__t3 = time()-__t3;"]);
+        eval (["__tid = tic();" __f2 "__t2 = toc(__tid);"]);
+        eval (["__tid = tic();" __f2 "__t3 = toc(__tid);"]);
         __t = min ([__t, __t2, __t3]);
       endif
       __torig(k) = __t;
@@ -287,19 +288,21 @@
   elseif (do_display)
 
     subplot (1, 2, 1);
-    semilogx (__test_n, __torig./__tnew,
-             ["-*r;" strrep(__f1, ";", ".") " / " strrep(__f2, ";", ".") ";"],
-              __test_n, __tnew./__torig,
-             ["-*g;", strrep(__f2, ";", ".") " / " strrep(__f1, ";", ".") ";"]);
+    semilogx (__test_n, __tnew./__torig, "-*g", 
+              __test_n, __torig./__tnew, "-*r");
+    legend ({[strrep(__f1, ";", ".") " / " strrep(__f2, ";", ".")],
+             [strrep(__f2, ";", ".") " / " strrep(__f1, ";", ".")]},
+            "location", "northwest");
     title ("Speedup Ratio");
     xlabel ("test length");
     ylabel ("speedup ratio");
 
     subplot (1, 2, 2);
-    loglog (__test_n, __tnew*1000,
-            ["*-g;" strrep(__f1,";",".") ";"],
-            __test_n, __torig*1000,
-            ["*-r;" strrep(__f2,";",".") ";"]);
+    loglog (__test_n, __tnew*1000, "*-g",
+            __test_n, __torig*1000, "*-r");
+    legend ({strrep(__f1,";","."),
+             strrep(__f2,";",".")},
+            "location", "northwest");
     title ({"Execution Times", ["init: " __init]});
     xlabel ("test length");
     ylabel ("best execution time (ms)");
@@ -430,7 +433,8 @@
 %! assert (length (T_f2) > 10);
 
 %!test
-%! [order, n, T_f1, T_f2] = speed ("sum (x)", "", [100, 1000], "v = 0; for i = 1:length (x), v += x(i); endfor");
+%! [order, n, T_f1, T_f2] = speed ("sum (x)", "", [100, 1000], ...
+%!                            "v = 0; for i = 1:length (x), v += x(i); endfor");
 %! assert (isstruct (order));
 %! assert (size (order), [1, 1]);
 %! assert (fieldnames (order), {"p"; "a"});