# HG changeset patch
# User David Bateman <dbateman@free.fr>
# Date 1371687445 -7200
# Node ID 8fce0ed4894ad8b8980e00f6a3e0d3f91bf37075
# Parent  e205f5ea826ac9b59078f2346d4f26516a2ddfd7
Specialize is_empty and numel methods for sparse matrices (debian bug #706376)

* ov-base.h (virtual bool is_empty (void) const) : Make method virtual
* ov-base-sparse.h (bool is_empty (void) const)) : Declare new method
(octave_idx_type numel (void) const): New method.
* ov-base-sparse.cc (template <class T> bool octave_base_sparse<T>:is_empty
(void) const)) : Define new method
* ov-bool-sparse.h (idx_vector index_vector (void) const): Start specialization
for sparse matrices
* lo-array-gripes.cc (void gripe_index_value (void)): Clarify error message
* data.cc : Add test to isempty
* idx-vector.cc (idx_vector::idx_vector_rep::idx_vector_rep (const
Sparse<bool>&) : Fix sparse indexing bug
* __sprand_impl__.m : Don't use randparm when there are more than sizemax() elements.
* sprand.m, sprandn.m : Add a test for large, low density matrices

diff --git a/libinterp/interpfcn/data.cc b/libinterp/interpfcn/data.cc
--- a/libinterp/interpfcn/data.cc
+++ b/libinterp/interpfcn/data.cc
@@ -3342,6 +3342,11 @@
   return retval;
 }
 
+/*
+%% Debian bug #706376
+%!assert (isempty (speye(2^16)), false)
+*/
+
 DEFUN (isnumeric, args, ,
   "-*- texinfo -*-\n\
 @deftypefn {Built-in Function} {} isnumeric (@var{x})\n\
diff --git a/libinterp/octave-value/ov-base-sparse.cc b/libinterp/octave-value/ov-base-sparse.cc
--- a/libinterp/octave-value/ov-base-sparse.cc
+++ b/libinterp/octave-value/ov-base-sparse.cc
@@ -278,6 +278,15 @@
 
 template <class T>
 bool
+octave_base_sparse<T>::is_empty (void) const
+{
+  dim_vector dv = dims ();
+
+  return (dv.any_zero ());
+}
+
+template <class T>
+bool
 octave_base_sparse<T>::print_as_scalar (void) const
 {
   dim_vector dv = dims ();
diff --git a/libinterp/octave-value/ov-base-sparse.h b/libinterp/octave-value/ov-base-sparse.h
--- a/libinterp/octave-value/ov-base-sparse.h
+++ b/libinterp/octave-value/ov-base-sparse.h
@@ -72,6 +72,8 @@
 
   ~octave_base_sparse (void) { }
 
+  octave_idx_type numel (void) const { return dims ().safe_numel (); }
+
   octave_idx_type nnz (void) const { return matrix.nnz (); }
 
   octave_idx_type nzmax (void) const { return matrix.nzmax (); }
@@ -137,6 +139,8 @@
 
   bool is_defined (void) const { return true; }
 
+  bool is_empty (void) const;
+
   bool is_constant (void) const { return true; }
 
   bool is_true (void) const;
diff --git a/libinterp/octave-value/ov-base.h b/libinterp/octave-value/ov-base.h
--- a/libinterp/octave-value/ov-base.h
+++ b/libinterp/octave-value/ov-base.h
@@ -331,7 +331,7 @@
 
   virtual bool is_defined (void) const { return false; }
 
-  bool is_empty (void) const { return numel () == 0; }
+  virtual bool is_empty (void) const { return numel () == 0; }
 
   virtual bool is_cell (void) const { return false; }
 
diff --git a/libinterp/octave-value/ov-bool-sparse.h b/libinterp/octave-value/ov-bool-sparse.h
--- a/libinterp/octave-value/ov-bool-sparse.h
+++ b/libinterp/octave-value/ov-bool-sparse.h
@@ -83,9 +83,9 @@
 
   octave_base_value *try_narrowing_conversion (void);
 
-  // FIXME Adapt idx_vector to allow sparse logical indexing!!
+  // FIXME Adapt idx_vector to allow sparse logical indexing without overflow!!
   idx_vector index_vector (void) const
-    { return idx_vector (bool_array_value ()); }
+    { return idx_vector (matrix); }
 
   builtin_type_t builtin_type (void) const { return btyp_bool; }
 
diff --git a/liboctave/array/idx-vector.cc b/liboctave/array/idx-vector.cc
--- a/liboctave/array/idx-vector.cc
+++ b/liboctave/array/idx-vector.cc
@@ -423,30 +423,27 @@
 }
 
 idx_vector::idx_vector_rep::idx_vector_rep (const Sparse<bool>& bnda)
-  : data (0), len (0), ext (0), aowner (0), orig_dims ()
+  : data (0), len (bnda.nnz ()), ext (0), aowner (0), orig_dims ()
 {
-  for (octave_idx_type i = 0, l = bnda.nnz (); i < l; i++)
-    if (bnda.data (i)) len++;
+  const dim_vector dv = bnda.dims ();
 
-  dim_vector dv = bnda.dims ();
-
-  orig_dims = ((dv.length () == 2 && dv(0) == 1)
-               ? dim_vector (1, len) : orig_dims = dim_vector (len, 1));
+  if (! dv.all_zero ())
+    orig_dims = ((dv.length () == 2 && dv(0) == 1)
+                 ? dim_vector (1, len) : dim_vector (len, 1));
 
   if (len != 0)
     {
       octave_idx_type *d = new octave_idx_type [len];
 
-      octave_idx_type nnz = bnda.nnz ();
+      octave_idx_type k = 0;
+      octave_idx_type nc = bnda.cols ();
+      octave_idx_type nr = bnda.rows ();
 
-      octave_idx_type k = 0;
-      // FIXME: I hope this is OK, i.e. the element iterated this way are correctly ordered.
-      for (octave_idx_type i = 0; i < nnz; i++)
-        {
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = bnda.cidx(j); i < bnda.cidx(j+1); i++)
           if (bnda.data (i))
-            d[k++] = bnda.cidx (i) + bnda.rows () * bnda.ridx (i);
-        }
-
+            d[k++] = j * nr + bnda.ridx (i);
+ 
       data = d;
 
       ext = d[k-1] + 1;
diff --git a/liboctave/util/lo-array-gripes.cc b/liboctave/util/lo-array-gripes.cc
--- a/liboctave/util/lo-array-gripes.cc
+++ b/liboctave/util/lo-array-gripes.cc
@@ -130,7 +130,11 @@
   const char *err_id = error_id_invalid_index;
 
   (*current_liboctave_error_with_id_handler)
-    (err_id, "subscript indices must be either positive integers or logicals");
+#ifdef USE_64_BIT_IDX_T
+    (err_id, "subscript indices must be either positive integers less than 2^63 or logicals");
+#else
+    (err_id, "subscript indices must be either positive integers less than 2^31 or logicals");
+#endif
 }
 
 // FIXME -- the following is a common error message to resize,
diff --git a/scripts/sparse/private/__sprand_impl__.m b/scripts/sparse/private/__sprand_impl__.m
--- a/scripts/sparse/private/__sprand_impl__.m
+++ b/scripts/sparse/private/__sprand_impl__.m
@@ -55,9 +55,27 @@
 
   mn = m*n;
   k = round (d*mn);
-  idx = randperm (mn, k);
+  if (mn > sizemax ())
+    ## randperm will overflow, so use alternative methods
+
+    idx = unique (fix (rand (min (k*1.01, k+10), 1) * mn)) + 1;
 
-  [i, j] = ind2sub ([m, n], idx);
+    ## idx contains random numbers in [1,mn]
+    ## generate 1% or 10 more random values than necessary in order to
+    ## reduce the probability that there are less than k distinct
+    ## values; maybe a better strategy could be used but I don't think
+    ## it's worth the price
+    
+    ## actual number of entries in S
+    k = min (length (idx), k);
+    j = floor ((idx(1:k) - 1) / m);
+    i = idx(1:k) - j * m;
+  else
+    idx = randperm (mn, k);
+    [i, j] = ind2sub ([m, n], idx);
+  endif
+
+
   S = sparse (i, j, randfun (k, 1), m, n);
 
 endfunction
diff --git a/scripts/sparse/sprand.m b/scripts/sparse/sprand.m
--- a/scripts/sparse/sprand.m
+++ b/scripts/sparse/sprand.m
@@ -81,3 +81,7 @@
 %!error sprand (3, 3, -1)
 %!error sprand (3, 3, 2)
 
+%% Test very large, very low density matrix doesn't fail 
+%!test
+%! s = sprand(1e6,1e6,1e-7);
+
diff --git a/scripts/sparse/sprandn.m b/scripts/sparse/sprandn.m
--- a/scripts/sparse/sprandn.m
+++ b/scripts/sparse/sprandn.m
@@ -72,3 +72,6 @@
 %!error sprandn (3, 3, -1)
 %!error sprandn (3, 3, 2)
 
+%% Test very large, very low density matrix doesn't fail 
+%!test
+%! s = sprandn(1e6,1e6,1e-7);