# HG changeset patch # User Rik # Date 1401684118 25200 # Node ID dccbc8bff5cb55c599df63c1f37db0860a8e72b4 # Parent 658d23da2c46622eb4c31822b480c5706f64734a 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. diff -r 658d23da2c46 -r dccbc8bff5cb libinterp/corefcn/data.cc --- a/libinterp/corefcn/data.cc Sun Jun 01 15:42:17 2014 +0200 +++ b/libinterp/corefcn/data.cc Sun Jun 01 21:41:58 2014 -0700 @@ -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) */ @@ -327,12 +358,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 (); @@ -397,6 +425,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 @@ -593,12 +650,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 (); @@ -619,26 +673,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) @@ -727,12 +807,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 (); diff -r 658d23da2c46 -r dccbc8bff5cb liboctave/util/oct-binmap.h --- 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 binmap (const T& x, const Sparse& ys, F fcn) { - octave_idx_type nz = ys.nnz (); - Sparse 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 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 (binmap (x, ys.array_value (), fcn)); } // Sparse-scalar @@ -236,17 +249,30 @@ Sparse binmap (const Sparse& xs, const R& y, F fcn) { - octave_idx_type nz = xs.nnz (); - Sparse 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 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 (binmap (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));