changeset 18982:03c2671493f9

maint: Periodic merge of gui-release to default.
author Rik <rik@octave.org>
date Mon, 02 Jun 2014 11:17:59 -0700
parents 96527fbdb8d4 (current diff) 332450e56698 (diff)
children 1514f5337781
files NEWS libinterp/corefcn/data.cc libinterp/corefcn/find.cc
diffstat 6 files changed, 224 insertions(+), 111 deletions(-) [+]
line wrap: on
line diff
--- a/NEWS
+++ b/NEWS
@@ -123,6 +123,11 @@
 
       allow_noninteger_range_as_index
       do_braindead_shortcircuit_evaluation
+      
+    The internal function atan2 of the sparse matrix class has been deprecated
+    in Octave 4.0 and will be removed from Octave 4.4 (or whatever version is
+    the second major release after 4.0).  Use the Fatan2 function with sparse
+    inputs as a replacement.
 
 ---------------------------------------------------------
 
--- a/libinterp/corefcn/data.cc
+++ b/libinterp/corefcn/data.cc
@@ -244,12 +244,9 @@
         }
       else
         {
-          bool a0_scalar = args(0).is_scalar_type ();
-          bool a1_scalar = args(1).is_scalar_type ();
-          if (a0_scalar && a1_scalar)
+          if (args(0).is_scalar_type () && args(1).is_scalar_type ())
             retval = atan2 (args(0).scalar_value (), args(1).scalar_value ());
-          else if ((a0_scalar || args(0).is_sparse_type ())
-                   && (a1_scalar || args(1).is_sparse_type ()))
+          else if (args(0).is_sparse_type ())
             {
               SparseMatrix m0 = args(0).sparse_matrix_value ();
               SparseMatrix m1 = args(1).sparse_matrix_value ();
@@ -292,6 +289,40 @@
 %! x = single ([1, 3, 1, 1, 1, 1, 3, 1]);
 %! assert (atan2 (y, x), v, sqrt (eps ("single")));
 
+## Test sparse implementations
+%!shared xs
+%! xs = sparse (0:3);
+%!test
+%! y = atan2 (1, xs);
+%! assert (issparse (y), false);
+%! assert (nnz (y), 4);
+%! assert (y, atan2 (1, 0:3));
+%!test
+%! y = atan2 (0, xs);
+%! assert (issparse (y), false);
+%! assert (nnz (y), 0);
+%! assert (y, zeros (1,4));
+%!test
+%! y = atan2 (xs, 1);
+%! assert (issparse (y));
+%! assert (nnz (y), 3);
+%! assert (y, sparse (atan2 (0:3, 1)));
+%!test
+%! y = atan2 (xs, 0);
+%! assert (issparse (y));
+%! assert (nnz (y), 3);
+%! assert (y, sparse (atan2 (0:3, 0)));
+%!test
+%! y = atan2 (xs, sparse (ones (1, 4)));
+%! assert (issparse (y));
+%! assert (nnz (y), 3);
+%! assert (y, sparse (atan2 (0:3, ones (1,4))));
+%!test
+%! y = atan2 (xs, sparse (zeros (1,4)));
+%! assert (issparse (y));
+%! assert (nnz (y), 3);
+%! assert (y, sparse (atan2 (0:3, zeros (1,4))));
+
 %!error atan2 ()
 %!error atan2 (1, 2, 3)
 */
@@ -328,12 +359,9 @@
         }
       else
         {
-          bool a0_scalar = arg0.is_scalar_type ();
-          bool a1_scalar = arg1.is_scalar_type ();
-          if (a0_scalar && a1_scalar)
+          if (arg0.is_scalar_type () && arg1.is_scalar_type ())
             retval = hypot (arg0.scalar_value (), arg1.scalar_value ());
-          else if ((a0_scalar || arg0.is_sparse_type ())
-                   && (a1_scalar || arg1.is_sparse_type ()))
+          else if (arg0.is_sparse_type () || arg1.is_sparse_type ())
             {
               SparseMatrix m0 = arg0.sparse_matrix_value ();
               SparseMatrix m1 = arg1.sparse_matrix_value ();
@@ -398,6 +426,35 @@
 %!assert (size (hypot (1, 2)), [1, 1])
 %!assert (hypot (1:10, 1:10), sqrt (2) * [1:10], 16*eps)
 %!assert (hypot (single (1:10), single (1:10)), single (sqrt (2) * [1:10]))
+
+## Test sparse implementations
+%!shared xs
+%! xs = sparse (0:3);
+%!test
+%! y = hypot (1, xs);
+%! assert (nnz (y), 4);
+%! assert (y, sparse (hypot (1, 0:3)));
+%!test
+%! y = hypot (0, xs);
+%! assert (nnz (y), 3);
+%! assert (y, xs);
+%!test
+%! y = hypot (xs, 1);
+%! assert (nnz (y), 4);
+%! assert (y, sparse (hypot (0:3, 1)));
+%!test
+%! y = hypot (xs, 0);
+%! assert (nnz (y), 3);
+%! assert (y, xs);
+%!test
+%! y = hypot (sparse ([0 0]), sparse ([0 1]));
+%! assert (nnz (y), 1);
+%! assert (y, sparse ([0 1]));
+%!test
+%! y = hypot (sparse ([0 1]), sparse ([0 0]));
+%! assert (nnz (y), 1);
+%! assert (y, sparse ([0 1]));
+
 */
 
 template<typename T, typename ET>
@@ -594,12 +651,9 @@
         }
       else
         {
-          bool a0_scalar = args(0).is_scalar_type ();
-          bool a1_scalar = args(1).is_scalar_type ();
-          if (a0_scalar && a1_scalar)
+          if (args(0).is_scalar_type () && args(1).is_scalar_type ())
             retval = xrem (args(0).scalar_value (), args(1).scalar_value ());
-          else if ((a0_scalar || args(0).is_sparse_type ())
-                   && (a1_scalar || args(1).is_sparse_type ()))
+          else if (args(0).is_sparse_type () || args(1).is_sparse_type ())
             {
               SparseMatrix m0 = args(0).sparse_matrix_value ();
               SparseMatrix m1 = args(1).sparse_matrix_value ();
@@ -620,26 +674,52 @@
 }
 
 /*
+%!assert (size (fmod (zeros (0, 2), zeros (0, 2))), [0, 2])
+%!assert (size (fmod (rand (2, 3, 4), zeros (2, 3, 4))), [2, 3, 4])
+%!assert (size (fmod (rand (2, 3, 4), 1)), [2, 3, 4])
+%!assert (size (fmod (1, rand (2, 3, 4))), [2, 3, 4])
+%!assert (size (fmod (1, 2)), [1, 1])
+
 %!assert (rem ([1, 2, 3; -1, -2, -3], 2), [1, 0, 1; -1, 0, -1])
 %!assert (rem ([1, 2, 3; -1, -2, -3], 2 * ones (2, 3)),[1, 0, 1; -1, 0, -1])
 %!assert (rem (uint8 ([1, 2, 3; -1, -2, -3]), uint8 (2)), uint8 ([1, 0, 1; -1, 0, -1]))
 %!assert (uint8 (rem ([1, 2, 3; -1, -2, -3], 2 * ones (2, 3))),uint8 ([1, 0, 1; -1, 0, -1]))
 
+## Test sparse implementations
+%!shared xs
+%! xs = sparse (0:3);
+%!test
+%! y = rem (11, xs);
+%! assert (nnz (y), 3);
+%! assert (y, sparse (rem (11, 0:3)));
+%!test
+%! y = rem (0, xs);
+%! assert (nnz (y), 0);
+%! assert (y, sparse (zeros (1,4)));
+%!test
+%! y = rem (xs, 2);
+%! assert (nnz (y), 2);
+%! assert (y, sparse (rem (0:3, 2)));
+%!test
+%! y = rem (xs, 1);
+%! assert (nnz (y), 0);
+%! assert (y, sparse (rem (0:3, 1)));
+%!test
+%! y = rem (sparse ([11 11 11 11]), xs);
+%! assert (nnz (y), 3);
+%! assert (y, sparse (rem (11, 0:3)));
+%!test
+%! y = rem (sparse ([0 0 0 0]), xs);
+%! assert (nnz (y), 0);
+%! assert (y, sparse (zeros (1,4)));
+
 %!error rem (uint (8), int8 (5))
 %!error rem (uint8 ([1, 2]), uint8 ([3, 4, 5]))
 %!error rem ()
 %!error rem (1, 2, 3)
 %!error rem ([1, 2], [3, 4, 5])
 %!error rem (i, 1)
-*/
-
-/*
-
-%!assert (size (fmod (zeros (0, 2), zeros (0, 2))), [0, 2])
-%!assert (size (fmod (rand (2, 3, 4), zeros (2, 3, 4))), [2, 3, 4])
-%!assert (size (fmod (rand (2, 3, 4), 1)), [2, 3, 4])
-%!assert (size (fmod (1, rand (2, 3, 4))), [2, 3, 4])
-%!assert (size (fmod (1, 2)), [1, 1])
+
 */
 
 DEFALIAS (fmod, rem)
@@ -728,12 +808,9 @@
         }
       else
         {
-          bool a0_scalar = args(0).is_scalar_type ();
-          bool a1_scalar = args(1).is_scalar_type ();
-          if (a0_scalar && a1_scalar)
+          if (args(0).is_scalar_type () && args(1).is_scalar_type ())
             retval = xmod (args(0).scalar_value (), args(1).scalar_value ());
-          else if ((a0_scalar || args(0).is_sparse_type ())
-                   && (a1_scalar || args(1).is_sparse_type ()))
+          else if (args(0).is_sparse_type () || args(1).is_sparse_type ())
             {
               SparseMatrix m0 = args(0).sparse_matrix_value ();
               SparseMatrix m1 = args(1).sparse_matrix_value ();
--- a/libinterp/corefcn/find.cc
+++ b/libinterp/corefcn/find.cc
@@ -184,7 +184,7 @@
     {
       // No items found.  Fixup return dimensions for Matlab compatibility.
       // The behavior to match is documented in Array.cc (Array<T>::find).
-      if ((nr == 0 && nc == 0) || nr == 1 & nc == 1)
+      if ((nr == 0 && nc == 0) || (nr == 1 && nc == 1))
         {
           idx.resize (0, 0);
 
@@ -294,12 +294,13 @@
     }
   else
     {
-      // FIXME: Is this case even possible?  A scalar permutation matrix seems to devolve
-      //        to a scalar full matrix, at least from the Octave command line.  Perhaps
-      //        this function could be called internally from C++ with such a matrix.
+      // FIXME: Is this case even possible?  A scalar permutation matrix seems
+      // to devolve to a scalar full matrix, at least from the Octave command
+      // line.  Perhaps this function could be called internally from C++ with
+      // such a matrix.
       // No items found.  Fixup return dimensions for Matlab compatibility.
       // The behavior to match is documented in Array.cc (Array<T>::find).
-      if ((nr == 0 && nc == 0) || nr == 1 & nc == 1)
+      if ((nr == 0 && nc == 0) || (nr == 1 && nc == 1))
         {
           idx.resize (0, 0);
 
--- a/liboctave/array/dSparse.h
+++ b/liboctave/array/dSparse.h
@@ -120,10 +120,13 @@
   friend OCTAVE_API SparseMatrix real (const SparseComplexMatrix& a);
   friend OCTAVE_API SparseMatrix imag (const SparseComplexMatrix& a);
 
-  friend OCTAVE_API SparseMatrix atan2 (const double& x, const SparseMatrix& y);
-  friend OCTAVE_API SparseMatrix atan2 (const SparseMatrix& x, const double& y);
+  friend OCTAVE_API SparseMatrix atan2 (const double& x, const SparseMatrix& y)
+                                        GCC_ATTR_DEPRECATED ;
+  friend OCTAVE_API SparseMatrix atan2 (const SparseMatrix& x, const double& y)
+                                        GCC_ATTR_DEPRECATED ; 
   friend OCTAVE_API SparseMatrix atan2 (const SparseMatrix& x,
-                                        const SparseMatrix& y);
+                                        const SparseMatrix& y)
+                                        GCC_ATTR_DEPRECATED ;
 
   SparseMatrix transpose (void) const
   {
--- a/liboctave/util/oct-binmap.h
+++ b/liboctave/util/oct-binmap.h
@@ -219,17 +219,30 @@
 Sparse<U>
 binmap (const T& x, const Sparse<R>& ys, F fcn)
 {
-  octave_idx_type nz = ys.nnz ();
-  Sparse<U> retval (ys.rows (), ys.cols (), nz);
-  for (octave_idx_type i = 0; i < nz; i++)
+  R yzero = R ();
+  U fz = fcn (x, yzero);
+
+  if (fz == U ())  // Sparsity preserving fcn
     {
+      octave_idx_type nz = ys.nnz ();
+      Sparse<U> retval (ys.rows (), ys.cols (), nz);
+      copy_or_memcpy (nz, ys.ridx (), retval.ridx ());
+      copy_or_memcpy (ys.cols () + 1, ys.cidx (), retval.cidx ());
+
+      for (octave_idx_type i = 0; i < nz; i++)
+        {
+          octave_quit ();
+          // FIXME: Could keep track of whether fcn call results in a 0.
+          //        If no zeroes are created could skip maybe_compress()
+          retval.xdata (i) = fcn (x, ys.data (i));
+        }
+
       octave_quit ();
-      retval.xdata (i) = fcn (x, ys.data (i));
+      retval.maybe_compress (true);
+      return retval;
     }
-
-  octave_quit ();
-  retval.maybe_compress ();
-  return retval;
+  else
+    return Sparse<U> (binmap<U, T, R, F> (x, ys.array_value (), fcn));
 }
 
 // Sparse-scalar
@@ -237,17 +250,30 @@
 Sparse<U>
 binmap (const Sparse<T>& xs, const R& y, F fcn)
 {
-  octave_idx_type nz = xs.nnz ();
-  Sparse<U> retval (xs.rows (), xs.cols (), nz);
-  for (octave_idx_type i = 0; i < nz; i++)
+  T xzero = T ();
+  U fz = fcn (xzero, y);
+
+  if (fz == U ())  // Sparsity preserving fcn
     {
+      octave_idx_type nz = xs.nnz ();
+      Sparse<U> retval (xs.rows (), xs.cols (), nz);
+      copy_or_memcpy (nz, xs.ridx (), retval.ridx ());
+      copy_or_memcpy (xs.cols () + 1, xs.cidx (), retval.cidx ());
+
+      for (octave_idx_type i = 0; i < nz; i++)
+        {
+          octave_quit ();
+          // FIXME: Could keep track of whether fcn call results in a 0.
+          //        If no zeroes are created could skip maybe_compress()
+          retval.xdata (i) = fcn (xs.data (i), y);
+        }
+
       octave_quit ();
-      retval.xdata (i) = fcn (xs.data (i), y);
+      retval.maybe_compress (true);
+      return retval;
     }
-
-  octave_quit ();
-  retval.maybe_compress ();
-  return retval;
+  else
+    return Sparse<U> (binmap<U, T, R, F> (xs.array_value (), y, fcn));
 }
 
 // Sparse-Sparse (treats singletons as scalars)
@@ -264,77 +290,61 @@
 
   T xzero = T ();
   R yzero = R ();
+  U fz = fcn (xzero, yzero);
 
-  U fz = fcn (xzero, yzero);
   if (fz == U ())
     {
-      // Sparsity-preserving function. Do it efficiently.
+      // Sparsity-preserving function.  Do it efficiently.
       octave_idx_type nr = xs.rows ();
       octave_idx_type nc = xs.cols ();
-      Sparse<T> retval (nr, nc);
+      Sparse<T> retval (nr, nc, xs.nnz () + ys.nnz ());
 
       octave_idx_type nz = 0;
-      // Count nonzeros.
       for (octave_idx_type j = 0; j < nc; j++)
         {
           octave_quit ();
-          octave_idx_type ix = xs.cidx (j);
-          octave_idx_type iy = ys.cidx (j);
-          octave_idx_type ux = xs.cidx (j+1);
-          octave_idx_type uy = ys.cidx (j+1);
-          while (ix != ux || iy != uy)
+
+          octave_idx_type jx = xs.cidx (j);
+          octave_idx_type jx_max = xs.cidx (j+1);
+          bool jx_lt_max = jx < jx_max;
+
+          octave_idx_type jy = ys.cidx (j);
+          octave_idx_type jy_max = ys.cidx (j+1);
+          bool jy_lt_max = jy < jy_max;
+
+          while (jx_lt_max || jy_lt_max)
             {
-              octave_idx_type rx = xs.ridx (ix);
-              octave_idx_type ry = ys.ridx (ix);
-              ix += rx <= ry;
-              iy += ry <= rx;
+              if (! jy_lt_max
+                  || (jx_lt_max && (xs.ridx (jx) < ys.ridx (jy))))
+                {
+                  retval.xridx (nz) = xs.ridx (jx);
+                  retval.xdata (nz) = fcn (xs.data (jx), yzero);
+                  jx++;
+                  jx_lt_max = jx < jx_max;
+                }
+              else if (! jx_lt_max
+                       || (jy_lt_max && (ys.ridx (jy) < xs.ridx (jx))))
+                {
+                  retval.xridx (nz) = ys.ridx (jy);
+                  retval.xdata (nz) = fcn (xzero, ys.data (jy));
+                  jy++;
+                  jy_lt_max = jy < jy_max;
+                }
+              else
+                {
+                  retval.xridx (nz) = xs.ridx (jx);
+                  retval.xdata (nz) = fcn (xs.data (jx), ys.data (jy));
+                  jx++;
+                  jx_lt_max = jx < jx_max;
+                  jy++;
+                  jy_lt_max = jy < jy_max;
+                }
               nz++;
             }
-
           retval.xcidx (j+1) = nz;
         }
 
-      // Allocate space.
-      retval.change_capacity (retval.xcidx (nc));
-
-      // Fill.
-      nz = 0;
-      for (octave_idx_type j = 0; j < nc; j++)
-        {
-          octave_quit ();
-          octave_idx_type ix = xs.cidx (j);
-          octave_idx_type iy = ys.cidx (j);
-          octave_idx_type ux = xs.cidx (j+1);
-          octave_idx_type uy = ys.cidx (j+1);
-          while (ix != ux || iy != uy)
-            {
-              octave_idx_type rx = xs.ridx (ix);
-              octave_idx_type ry = ys.ridx (ix);
-              if (rx == ry)
-                {
-                  retval.xridx (nz) = rx;
-                  retval.xdata (nz) = fcn (xs.data (ix), ys.data (iy));
-                  ix++;
-                  iy++;
-                }
-              else if (rx < ry)
-                {
-                  retval.xridx (nz) = rx;
-                  retval.xdata (nz) = fcn (xs.data (ix), yzero);
-                  ix++;
-                }
-              else if (ry < rx)
-                {
-                  retval.xridx (nz) = ry;
-                  retval.xdata (nz) = fcn (xzero, ys.data (iy));
-                  iy++;
-                }
-
-              nz++;
-            }
-        }
-
-      retval.maybe_compress ();
+      retval.maybe_compress (true);
       return retval;
     }
   else
--- a/scripts/sparse/eigs.m
+++ b/scripts/sparse/eigs.m
@@ -359,7 +359,7 @@
 
     v = args{1};
     v = v(:,idx);
-    out{1} = v(selection,:);
+    out{1} = v(:,selection);
   endif
 
 endfunction
@@ -1114,8 +1114,25 @@
 %! [~, idx] = sort (abs (reseig), "ascend");
 %! assert (eigs (A, B, 10, 0), reseig (idx))
 
+%!test
+%! X = [70 47 42 39 50 73 79 23;
+%!      19 52 61 80 36 76 63 68;
+%!      14 34 66 65 29  4 72  9;
+%!      24  8 78 49 58 54 43 33;
+%!      62 69 32 31 40 46 22 28;
+%!      48 12 45 59 10 17 15 25;
+%!      64 67 77 56 13 55 41 74;
+%!      37 38 18 21 11  3 71  7;
+%!       5 35 16  1 51 27 26 44;
+%!      30 57 60 75  2 53 20  6];
+%! Z = X * X';
+%! r = rank (Z);
+%! assert (r, 8);
+%! [V, D] = eigs (Z, r, "lm");
+%! 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]);
-
-