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

[Jordi Gutiérrez Hermoso]

On 19 June 2013 14:06, David Bateman <david at bateman.eu> wrote:
> 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 ?

We are assuming it in our code. In numel for one. And in places like
whatever processes A(idx) for a logical index idx. We're not making
special cases in these places, "if (sparse_type)
dont_linearly_index();"

Each of these places that assume that linear indexing works needs to
be patched to check for sparse types.

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

We are fixing one symptom of a larger bug, a bug that is present in
many different locations.

>>  but this is just treating the symptoms, not the disease.
>
> So you advocate everyone moving to 64-bit indexing ?

That would delay the problem for a nontrivial amount of time. It would
be nice, but it wouldn't fix the problem.

There are other things we could do.

(1) Avoid linear indexing for sparse matrices as much as possible,
i.e. check everywhere we can think of for sparse matrix types. You've
mentioned a few more places where this should be checked.

(2) Warn when creating sparse matrices with large indices that some
operations may not work, or clearly error out when those operations
are attempted.

(3) Abstract octave_idx_type so that it doesn't actually use 32-bit
ints for sparse matrices.

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

Well, numel needs "if(sparse_type) {weep();}" or whatever.

>>> 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.

Yes, the real bug: that we have a tacit assumption in our code that
linear indexing works.

>>  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.

Yep, I suppose that's my dont_linearly_index() function suggested
above.

> 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.

I think you'll still need to do something about a more realistic
situation of linear indexing with a logical matrix, which also ends up
translating into linear indexing thanks to our underlying bug: the
assumption that linear indexing works. In this case, there shouldn't
be an error at all, like Ed suggested, since we have enough
information in a logical matrix to avoid linear indexing.

> See the attached changeset

It looks good. Can you push it? Also, note that we have an actual
Octave bug for it:

    https://savannah.gnu.org/bugs/?38936

Contrary to apperances, I don't mean to be unhelpful, so please let me
know if you can't push the fix yourself.

- Jordi G. H.