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

[Jordi Gutiérrez Hermoso]

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.