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

[Jordi Gutiérrez Hermoso]

On 29 April 2013 12:40, Ed Meyer <address@hidden> wrote:
>
>
> On Mon, Apr 29, 2013 at 6:10 AM, Jordi Gutiérrez Hermoso
> <address@hidden> wrote:
>>
>> On 29 April 2013 06:25, Miroslaw Kwasniak <address@hidden>
>> wrote:
>> > it's something wrong whith sparse matrices A(n,n) when n is a multiple
>> > of 65536=2^16.
>> >
>> > Demonstration code ======================================
>> >
>> > for i=1:3;
>> >   for n=i*2^16+(-1:1);
>> >     A=spdiags(ones(n,1),0,n,n);
>> >     t=trace(A);
>> >     printf("n=%8d trace=%8d %s\n",n,t,["ERR";"ok"]((t==n)+1,:));
>> >   endfor;
>> > endfor
>> >
>> > Results ======================================
>> >
>> > n=   65535 trace=   65535 ok
>> > n=   65536 trace=       0 ERR
>> > n=   65537 trace=   65537 ok
>> > n=  131071 trace=  131071 ok
>> > n=  131072 trace=       0 ERR
>> > n=  131073 trace=  131073 ok
>> > n=  196607 trace=  196607 ok
>> > n=  196608 trace=       0 ERR
>> > n=  196609 trace=  196609 ok
>>
>> Confirmed. The problem is that the numel function is limited to
>> returning octave_idx_type, which ordinarily of size 2^32, and
>> certainly is so for Debian. This makes sense, since you can only index
>> that many elements in a matrix. You're hitting the indexing limit. To
>> get 64-bit indexing, you would need to recompile all of Octave's
>> Fortran dependencies with -fdefault-integer-8.
>>
>> I'm not sure exactly what the bug is here. For instance, you can't
>> index your matrix A either, and this is checked for correctly:
>>
>>     A(end)
>>
>> Perhaps the best thing  to do would be to forbid creation of sparse
>> matrices where numel(A) > std::numeric_limits<int>::max().
>>
>> Your matrix is simply too large to be indexed, and this breaks
>> assumptions elsewhere in our code.
>>
>> - Jordi G. H.
>
>
> I'm confused - this is a diagonal sparse matrix so you should be able to
> trace() (or any other op)
> up to n = 2^32, not n^2 = 2^32. The limit on sparse matrices should be
> number of non-zeros < 2^32

All matrices need to be linearly indexable, and of course, this is how
they are actually stored in memory, as a single long array indexed by
a single index. Thus, the total number of indexable elements of a
matrix can't be larger than
std::numeric_limits<octave_idx_type>::max().

There could be some tricks we could do to relax this requirement for
sparse matrices, but it would require some pretty deep surgery of the
current code.

- Jordi G. H.