Mercurial > octave

diff liboctave/util/oct-binmap.h @ 18789:dccbc8bff5cb stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix binmap for sparse-scalar or scalar-sparse operations (bug #40813). * oct-binmap.h (binmap (Sparse, Scalar), binmap (Scalar, Sparse)): Check that function is sparsity preserving before using sparse algorithm. Initialize retval row, column indices from original sparse array. Call maybe_compress (true) to remove zero elements that function may have produced. * data.cc (Fatan2): Only preserve sparsity if *first* argument is sparse. Add %!tests to verify sparse operations. * data.cc (do_hypot): Call binmap with sparse inputs whenever at least one of the inputs is sparse. * data.cc (Fhypot): Add %!tests to verify sparse operations. * data.cc (Frem): Call binmap with sparse inputs whenever at least one of the inputs is sparse. Add %!tests to verify sparse operations.
author Rik <rik@octave.org>
date Sun, 01 Jun 2014 21:41:58 -0700
parents 658d23da2c46
children 59975c3cea6b
line wrap: on
line diff
--- a/liboctave/util/oct-binmap.h	Sun Jun 01 15:42:17 2014 +0200
+++ b/liboctave/util/oct-binmap.h	Sun Jun 01 21:41:58 2014 -0700
@@ -218,17 +218,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
@@ -236,17 +249,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)
@@ -295,8 +321,8 @@
                   jx++;
                   jx_lt_max = jx < jx_max;
                 }
-              else if (! jx_lt_max ||
-                  (jy_lt_max && (ys.ridx (jy) < xs.ridx (jx))))
+              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));