[Pkg-octave-devel] Bug#706376: Bug#706376: Bug#706376: Bug#706376: Bug#706376: octave: sparse matrix n*2^16

[David Bateman]

On 06/16/2013 03:59 AM, Jordi Gutiérrez Hermoso wrote:
> I'm saying the real problem is that we assume linear indexing works
> for all matrix types, including sparse matrices. I claim that this is
> the real problem. 
Who is assuming linear indexing works for all matrix types ? Where
exactly is that stated ? If its in then manual that its a bug in the
manual as linear indexing can never work correctly in this case because
the underlying 32 integer type won't allowed it and in fact octave as it
stand throws an error !!!

> We can patch around this problem by avoiding linear
> indexing,

The bug report was in "trace" that called "isempty" on a sparse matrix.
Neither function needs "numel" or "linear indexing". We aren't patching
around anything, we are fixing a bug

>  but this is just treating the symptoms, not the disease.

So you advocate everyone moving to 64-bit indexing ? In that case what
happens when we get the same bug report for s(2^32, 2^32) ? Is that a
"disease" too.  Certainly there is a limitation on linear indexing and
numel for sparse matrices. We should document it and make as many things
as possible work with sparse matrices and be done with it.

> While I don't deny that we can make some progress masking the
> symptoms, the disease itself should also be treated somehow.

There is no disease, and unless you want to artificially limit the size
of sparse matrices that can be treated such that numel is less than 2^31
for 32 bit indexing and 2^63 for 64 bit indexing. Why do this which
makes sparse matrices much less useful, so there is no solution for what
you call a disease

>> So essentially you're saying that sparse matrices with
>> 32-bit indexing and numel larger than 2^31 are useless!!
> I'm saying that they will fail in other unexpected ways,

Isn't that the definition of a bug.

>  and we shouln't mask symptoms. 

We never tried to. Look at the code in dim-vector.h

<quote>
  // The following function will throw a std::bad_alloc ()
  // exception if the requested size is larger than can be indexed by
  // octave_idx_type. This may be smaller than the actual amount of
  // memory that can be safely allocated on a system.  However, if we
  // don't fail here, we can end up with a mysterious crash inside a
  // function that is iterating over an array using octave_idx_type
  // indices.

  octave_idx_type safe_numel (void) const;
</quote>

The numel method of Sparse<T> calls this method that is supposed to
throw an error. However as the builtin Fnumel is calling args(1).numel()
which is calling dims().numel() the sparse safe version of numel isn't
being called. The solution is to add a numel method to
octave_base_sparse that calls dims().safe_numel() instead. So this is a
bug as well.

As for linear indexing, if you look in idx-vector.cc you'll see that in
the convert_index functions an error is returned like

octave_idx_type idx = static_cast<octave_idx_type>(d)
bool err = static_cast<double> (idx) != d;

So as expected

s = speye (2^17);
s (2^32)

throws an error

error: subscript indices must be either positive integers or logicals

You might think this is a little cryptic but I wouldn't call it "masking
a symptom". I propose modifying this error to read

error: subscript indices must be either positive integers less than 2^31
or logicals

for 32 bit indexing and

error: subscript indices must be either positive integers less than 2^63
or logicals

for 64 bit indexing. See the attached changeset

David



-------------- next part --------------
# HG changeset patch
# User David Bateman <dbateman at free.fr>
# Date 1371664645 -7200
# Node ID 95d4850cddd551c2747855aea8be357ab6f84bac
# Parent  3542e106c496f04c339946ff6b9d67c3c5305e7f
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
* lo-array-gripes.cc (vois gripe_index_value (void)): Clarify error message

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/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,