Bug#706376: [Pkg-octave-devel] Bug#706376: Bug#706376: Bug#706376: Bug#706376: Bug#706376: Bug#706376: Bug#706376: octave: sparse matrix n*2^16
On 06/20/2013 04:50 PM, Jordi Gutiérrez Hermoso wrote: I prefer the principle of consistency and as few special cases as possible. More special cases increase code complexity and makes it more difficult to maintain. - Jordi G. H. Ok so a partially undid this changeset and used the sparse specialize version of the is_empty method to replace the version in ov-base.h D. -- To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: Bug#706376: Bug#706376: Bug#706376: Bug#706376: Bug#706376: octave: sparse matrix n*2^16
On 06/19/2013 08:32 PM, David Bateman wrote: On 06/20/2013 01:10 AM, David Bateman wrote: I'd like to add some tests first and see if any other bugs have turned up after this change. For example the changes use made to sprand and sprandn 2 years ago to call randperm also overflows. At the moment I'm getting 791 failed tests with make check is that normal ? David Ok, it seems Jaroslav's code for idx_vector(Sparsebool hasn't been used much in the last 5 years as it was completely wrong and when I started using it, it caused 791 failures in make check. I've fixed his code as it makes sense to use it and pushed my changeset at http://hg.savannah.gnu.org/hgweb/octave/rev/8fce0ed4894a This change seems OK to me, but is there some reason to not use dim_vector dv = dims (); return (dv.any_zero ()); as the default definition for is_empty? jwe -- To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: Bug#706376: Bug#706376: Bug#706376: Bug#706376: Bug#706376: octave: sparse matrix n*2^16
Le 20 juin 2013 à 15:01, John W. Eaton j...@octave.org a écrit : On 06/19/2013 08:32 PM, David Bateman wrote: On 06/20/2013 01:10 AM, David Bateman wrote: I'd like to add some tests first and see if any other bugs have turned up after this change. For example the changes use made to sprand and sprandn 2 years ago to call randperm also overflows. At the moment I'm getting 791 failed tests with make check is that normal ? David Ok, it seems Jaroslav's code for idx_vector(Sparsebool hasn't been used much in the last 5 years as it was completely wrong and when I started using it, it caused 791 failures in make check. I've fixed his code as it makes sense to use it and pushed my changeset at http://hg.savannah.gnu.org/hgweb/octave/rev/8fce0ed4894a This change seems OK to me, but is there some reason to not use dim_vector dv = dims (); return (dv.any_zero ()); as the default definition for is_empty? jwe Just the principal of minimum impact on existing code. Making this change everywhere will make this slightly faster too, though I doubt the calculation of numel is the limiting factor in most codes D. -- To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: Bug#706376: Bug#706376: Bug#706376: Bug#706376: Bug#706376: Bug#706376: octave: sparse matrix n*2^16
On 20 June 2013 10:41, David Bateman adb...@gmail.com wrote: Le 20 juin 2013 à 15:01, John W. Eaton j...@octave.org a écrit : On 06/19/2013 08:32 PM, David Bateman wrote: Ok, it seems Jaroslav's code for idx_vector(Sparsebool hasn't been used much in the last 5 years as it was completely wrong and when I started using it, it caused 791 failures in make check. I've fixed his code as it makes sense to use it and pushed my changeset at http://hg.savannah.gnu.org/hgweb/octave/rev/8fce0ed4894a This change seems OK to me, but is there some reason to not use dim_vector dv = dims (); return (dv.any_zero ()); as the default definition for is_empty? Just the principal of minimum impact on existing code. I prefer the principle of consistency and as few special cases as possible. More special cases increase code complexity and makes it more difficult to maintain. - Jordi G. H. -- To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: Bug#706376: Bug#706376: Bug#706376: octave: sparse matrix n*2^16
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 SparseT 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_castoctave_idx_type(d)
bool err = static_castdouble (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
# HG changeset patch
# User David Bateman dbate...@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_sparseT: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_sparseT::is_empty (void) const
+{
+ dim_vector dv = dims ();
+
+ return (dv.any_zero ());
+}
+
+template class T
+bool
octave_base_sparseT::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 (); }
@@
Bug#706376: [Pkg-octave-devel] Bug#706376: Bug#706376: Bug#706376: Bug#706376: Bug#706376: octave: sparse matrix n*2^16
On 19 June 2013 14:06, David Bateman da...@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 SparseT 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_castoctave_idx_type(d)
bool err = static_castdouble (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.
--
To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: Bug#706376: Bug#706376: Bug#706376: Bug#706376: octave: sparse matrix n*2^16
On 06/19/2013 09:53 PM, Jordi Gutiérrez Hermoso wrote:
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.
Huh ? How can a logical matrix linear index overflow ? The matrix is
limited to 2^31-1 elements in any case and so can never overflow! I
presume you means sparse logical matrix linear indexing. This is one
of the FIXMEs the the sparse code that has existed for a long while and
is documented in ov-bool-sparse.h
quote
// FIXME Adapt idx_vector to allow sparse logical indexing!!
idx_vector index_vector (void) const
{ return idx_vector (bool_array_value ()); }
/quote
and in the projects page (as copied from the old PROJECTS file)
quote
Sparse logical indexing in idx_vector class so that something like
'a=sprandn(1e6,1e6,1e-6); a(a1) = 0' won't cause a memory overflow.
/quote
At the time I wrote the above line, it seemed a major effort to do this
correct. It seems that Yaroslav added some code to treat sparse logical
index vectors but didn't connect it to the idx_vector method of
octave_sparse_bool_matrix. It should be used as it'll improve the speed
of the creation of idx_vector, but it doesn't help as the underlying
data type of idx_vector in this case is Arrayoctave_idx_type and so a
single index (even if its a sparse logical matrix) can result in an
overflow. This will take major reworking of the idx_vector class and the
uses of single index vectors in SparseT. In any case this throws and
error at this point and so won't silently fail
It looks good. Can you push it?
I'd like to add some tests first and see if any other bugs have turned
up after this change. For example the changes use made to sprand and
sprandn 2 years ago to call randperm also overflows. At the moment I'm
getting 791 failed tests with make check is that normal ?
David
--
To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org
with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: Bug#706376: Bug#706376: Bug#706376: Bug#706376: octave: sparse matrix n*2^16
On 06/20/2013 01:10 AM, David Bateman wrote: I'd like to add some tests first and see if any other bugs have turned up after this change. For example the changes use made to sprand and sprandn 2 years ago to call randperm also overflows. At the moment I'm getting 791 failed tests with make check is that normal ? David Ok, it seems Jaroslav's code for idx_vector(Sparsebool hasn't been used much in the last 5 years as it was completely wrong and when I started using it, it caused 791 failures in make check. I've fixed his code as it makes sense to use it and pushed my changeset at http://hg.savannah.gnu.org/hgweb/octave/rev/8fce0ed4894a David -- To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: Bug#706376: Bug#706376: octave: sparse matrix n*2^16
On 05/06/2013 04:24 PM, Jordi Gutiérrez Hermoso wrote:
On 4 May 2013 20:50, Ed Meyer eem2...@gmail.com wrote:
I think I see why numel() is getting called from trace(). isempty()
is called which calls is_empty() which should be virtual. Were it to
call the sparse version of is_empty() there would be no problem
because it tests the row column dimensions instead of calling
numel(). So shouldn't is_empty() be virtual?
Sure, this patches this hole. What about linear indexing, how will you
patch that one?
- Jordi G. H.
That a bit of a specious argument. Because we can't solve problem B we
shouldn't solve problem A. Taking this argument to the absurd this
shouldn't work either
n = 2^16;
s = sparse (1:n,1:n,1);
t = s * s;
as t will have more elements than can be indexed with octave_idx_type.
Clear it works. So essentially you're saying that sparse matrices with
32-bit indexing and numel larger than 2^31 are useless!!
A lot of attention was made in the sparse implementation to not rely
either on linear indexing or the value of numel to avoid this issue. I'd
think that any function or operator that relies on either for sparse
matrices, or at least when it doesn't have to, is buggy.
Ed is right the is_empty method of the octave_base_value class should be
specialized for sparse matrices as in the attached changeset as it has
no need to rely on numel. Then any function that relies on isempty
should now work for sparse matrices.
David
# HG changeset patch
# User David Bateman dbate...@free.fr
# Date 1371308306 -7200
# Node ID d5c61e2aefa428e9b51d5671ff1e72e7de90a500
# Parent 1c8b6ab2c8ae71f1fdef51016472396a6665f483
Specialize is_empty method for sparse matrices
* ov-base.h (virtual bool is_empty (void) const) : Make method virtual
* ov-base-sparse.h (bool is_empty (void) const)) : Declare new method
* ov-base-sparse.cc (template class T bool octave_base_sparseT:is_empty
(void) const)) : Define new method
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_sparseT::is_empty (void) const
+{
+ dim_vector dv = dims ();
+
+ return (dv.any_zero ());
+}
+
+template class T
+bool
octave_base_sparseT::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
@@ -137,6 +137,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; }
Bug#706376: [Pkg-octave-devel] Bug#706376: Bug#706376: Bug#706376: Bug#706376: octave: sparse matrix n*2^16
On 15 June 2013 11:00, David Bateman da...@bateman.eu wrote: Ed is right the is_empty method of the octave_base_value class should be specialized for sparse matrices as in the attached changeset as it has no need to rely on numel. Then any function that relies on isempty should now work for sparse matrices. Now what should happen if numel is actually called for sparse matrices? And how do you propose to fix linear indexing? I'm not saying you shouldn't do it, just that it's difficult and boring. Also, can't you push your own patches? I thought you were still in the access list. - Jordi G. H. -- To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: Bug#706376: Bug#706376: Bug#706376: octave: sparse matrix n*2^16
On 15 June 2013 11:00, David Bateman da...@bateman.eu wrote: That a bit of a specious argument. Because we can't solve problem B we shouldn't solve problem A. Taking this argument to the absurd this shouldn't work either n = 2^16; s = sparse (1:n,1:n,1); t = s * s; You're strawmanning me. 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. We can patch around this problem by avoiding linear indexing, but this is just treating the symptoms, not the disease. While I don't deny that we can make some progress masking the symptoms, the disease itself should also be treated somehow. 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, and we shouln't mask symptoms. By all means, avoiding calling numel whenever possible, but don't mask symptoms. Should numel or similar actually get called for such sparse matrices, at the very least there should be a warning. I would prefer if the warning occurred at the time that the sparse matrix is created, until we can fix the actual underlying problem. A lot of attention was made in the sparse implementation to not rely either on linear indexing or the value of numel to avoid this issue. It's great you avoided that, but the users of the sparse matrix type might not immediately think of avoiding it. - Jordi G. H. -- To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: Bug#706376: Bug#706376: octave: sparse matrix n*2^16
On 4 May 2013 20:50, Ed Meyer eem2...@gmail.com wrote: I think I see why numel() is getting called from trace(). isempty() is called which calls is_empty() which should be virtual. Were it to call the sparse version of is_empty() there would be no problem because it tests the row column dimensions instead of calling numel(). So shouldn't is_empty() be virtual? Sure, this patches this hole. What about linear indexing, how will you patch that one? - Jordi G. H. -- To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: Bug#706376: octave: sparse matrix n*2^16
On Tue, Apr 30, 2013 at 6:37 AM, Jordi Gutiérrez Hermoso jord...@octave.org wrote: On 29 April 2013 17:51, Ed Meyer eem2...@gmail.com wrote: I'm not saying numel() is to blame or should be changed, only that I see no reason to ever use numel when handling sparse matrices unless you are converting it to full in which case the current behavior is ok. Well, look at the current implementation of trace (), in which numel is a perfectly reasonable function to call. If you don't want to ever call numel () for sparse matrices, then you would need to rewrite this function to check for sparsity as well as any other function that calls numel (). I think I see why numel() is getting called from trace(). isempty() is called which calls is_empty() which should be virtual. Were it to call the sparse version of is_empty() there would be no problem because it tests the row column dimensions instead of calling numel(). So shouldn't is_empty() be virtual? -- Ed Meyer
Bug#706376: [Pkg-octave-devel] Bug#706376: Bug#706376: octave: sparse matrix n*2^16
On 29 April 2013 17:51, Ed Meyer eem2...@gmail.com wrote: I'm not saying numel() is to blame or should be changed, only that I see no reason to ever use numel when handling sparse matrices unless you are converting it to full in which case the current behavior is ok. Well, look at the current implementation of trace (), in which numel is a perfectly reasonable function to call. If you don't want to ever call numel () for sparse matrices, then you would need to rewrite this function to check for sparsity as well as any other function that calls numel (). Here again, why would you ever want A(idx) for a sparse matrix? Because this is the only way to do arbitrary indexing, say, indexing with a logical index: n = 2^16; A = sprandsym (n, 1e-7); idx = A 0.5; A(idx) *= -1; The alternative to arbitrary indexing is looping. I agree that a special index type is the wrong approach; what I'm saying is that with sparse matrices you should never run into this problem in the first place if you don't try to treat them the same as full. Otherwise why have sparse matrices at all? It is desirable to have sparse matrices behave like ordinary matrices under most circumstances, such as indexing and when writing the trace function. Otherwise, you would have to write special code all over the place to make sure that if the matrix is sparse, you don't linearly index it nor do you call numel or who knows what other functions that I haven't thought about yet. - Jordi G. H. -- To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: Bug#706376: octave: sparse matrix n*2^16
On Tue, Apr 30, 2013 at 6:37 AM, Jordi Gutiérrez Hermoso jord...@octave.org wrote: On 29 April 2013 17:51, Ed Meyer eem2...@gmail.com wrote: I'm not saying numel() is to blame or should be changed, only that I see no reason to ever use numel when handling sparse matrices unless you are converting it to full in which case the current behavior is ok. Well, look at the current implementation of trace (), in which numel is a perfectly reasonable function to call. If you don't want to ever call numel () for sparse matrices, then you would need to rewrite this function to check for sparsity as well as any other function that calls numel (). I don't see numel() used in trace(), nor in diag() which it calls - what am I missing? Here again, why would you ever want A(idx) for a sparse matrix? Because this is the only way to do arbitrary indexing, say, indexing with a logical index: n = 2^16; A = sprandsym (n, 1e-7); idx = A 0.5; A(idx) *= -1; The alternative to arbitrary indexing is looping. or if octave_idx_type were a class hierarchy with a (row,col) class in addition to linear indexing I agree that a special index type is the wrong approach; what I'm saying is that with sparse matrices you should never run into this problem in the first place if you don't try to treat them the same as full. Otherwise why have sparse matrices at all? It is desirable to have sparse matrices behave like ordinary matrices under most circumstances, such as indexing and when writing the trace function. Otherwise, you would have to write special code all over the place to make sure that if the matrix is sparse, you don't linearly index it nor do you call numel or who knows what other functions that I haven't thought about yet. - Jordi G. H. Not only is it desirable to have sparse and full matrices behave similarly, I believe the user should not need to be aware of which storage format is used so functions like eig() would work for either. The key is to use the C++ class system to have different implementations for each storage format. -- Ed Meyer
Bug#706376: [Pkg-octave-devel] Bug#706376: Bug#706376: Bug#706376: octave: sparse matrix n*2^16
On 30 April 2013 12:56, Ed Meyer eem2...@gmail.com wrote: On Tue, Apr 30, 2013 at 6:37 AM, Jordi Gutiérrez Hermoso jord...@octave.org wrote: On 29 April 2013 17:51, Ed Meyer eem2...@gmail.com wrote: I'm not saying numel() is to blame or should be changed, only that I see no reason to ever use numel when handling sparse matrices unless you are converting it to full in which case the current behavior is ok. Well, look at the current implementation of trace (), in which numel is a perfectly reasonable function to call. If you don't want to ever call numel () for sparse matrices, then you would need to rewrite this function to check for sparsity as well as any other function that calls numel (). I don't see numel() used in trace(), nor in diag() which it calls - what am I missing? Sorry, I was going from memory of a gdb session yesterday. The C++ isempty function ends up calling numel. I really don't see another way to reliably implement this function. Here again, why would you ever want A(idx) for a sparse matrix? Because this is the only way to do arbitrary indexing, say, indexing with a logical index: n = 2^16; A = sprandsym (n, 1e-7); idx = A 0.5; A(idx) *= -1; The alternative to arbitrary indexing is looping. or if octave_idx_type were a class hierarchy with a (row,col) class in addition to linear indexing This seems again like that special indexing type that seems to me like a lot of boring work. I agree that a special index type is the wrong approach; what I'm saying is that with sparse matrices you should never run into this problem in the first place if you don't try to treat them the same as full. Otherwise why have sparse matrices at all? It is desirable to have sparse matrices behave like ordinary matrices under most circumstances, such as indexing and when writing the trace function. Otherwise, you would have to write special code all over the place to make sure that if the matrix is sparse, you don't linearly index it nor do you call numel or who knows what other functions that I haven't thought about yet. The key is to use the C++ class system to have different implementations for each storage format. This again sounds like support for a special indexing type just for sparse matrices. - Jordi G. H. -- To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: Bug#706376: octave: sparse matrix n*2^16
On 04/30/2013 12:56 PM, Ed Meyer wrote: Not only is it desirable to have sparse and full matrices behave similarly, I believe the user should not need to be aware of which storage format is used so functions like eig() would work for either. The key is to use the C++ class system to have different implementations for each storage format. I haven't been following this thread closely and I haven't thought much about the details but I have no objection to trying to do a better job with handling numel and dimensions/indices generally. Is there some way we can get the better behavior in a minimally invasive way? Even if it requires significant changes, maybe we should consider what the options are anyway. What changes are needed to make octave_idx_type behave the way you way you want? jwe -- To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: Bug#706376: octave: sparse matrix n*2^16
On 04/30/2013 01:19 PM, John W. Eaton wrote: On 04/30/2013 12:56 PM, Ed Meyer wrote: Not only is it desirable to have sparse and full matrices behave similarly, I believe the user should not need to be aware of which storage format is used so functions like eig() would work for either. The key is to use the C++ class system to have different implementations for each storage format. I haven't been following this thread closely and I haven't thought much about the details but I have no objection to trying to do a better job with handling numel and dimensions/indices generally. Is there some way we can get the better behavior in a minimally invasive way? Even if it requires significant changes, maybe we should consider what the options are anyway. What changes are needed to make octave_idx_type behave the way you way you want? jwe It would be good to give some thought to the trade-off of fixing the current system against completing the 64bit compiling system. The 64bit system will be useful in other ways and should be done in any case. Of course, any improvement in the current system is good, but the full implementation of large matrices could be a 64bit system feature. Jordi seemed to think that 2^63 indicies would work with no trouble in the 64bit system... Michael -- To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: octave: sparse matrix n*2^16
Package: octave Version: 3.6.2-5 Severity: normal Dear Maintainer, 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 == It isn't a bug in spdiags, but in the sparse matrices handling, the same results are given when A=spdiags(ones(n,1),0,n,n); in above code is replaced with (slower): B=sparse(n,n);for j=1:n;B(j,j)=1;endfor Miroslaw Kwasniak -- System Information: Debian Release: 7.0 APT prefers testing APT policy: (990, 'testing'), (500, 'testing-updates'), (500, 'testing-proposed-updates'), (500, 'unstable'), (1, 'experimental') Architecture: amd64 (x86_64) Foreign Architectures: i386 Kernel: Linux 3.2.0-4-amd64 (SMP w/2 CPU cores) Locale: LANG=pl_PL.UTF-8, LC_CTYPE=pl_PL.UTF-8 (charmap=UTF-8) Shell: /bin/sh linked to /bin/dash Versions of packages octave depends on: ii libamd2.2.0 1:3.4.0-3 ii libarpack2 3.1.1-2.1 ii libatlas3-base [liblapack.so.3] 3.8.4-9 ii libblas3 [libblas.so.3] 1.2.20110419-5 ii libc62.13-38 ii libcamd2.2.0 1:3.4.0-3 ii libccolamd2.7.1 1:3.4.0-3 ii libcholmod1.7.1 1:3.4.0-3 ii libcolamd2.7.1 1:3.4.0-3 ii libcurl3-gnutls 7.26.0-1+wheezy2 ii libcxsparse2.2.3 1:3.4.0-3 ii libfftw3-3 3.3.2-3.1 ii libfltk1.1 1.1.10-14 ii libfreetype6 2.4.9-1.1 ii libgcc1 1:4.7.2-5 ii libgl1-mesa-glx [libgl1] 8.0.5-4 ii libglpk0 4.45-1 ii libgomp1 4.7.2-5 ii libgraphicsmagick++3 1.3.16-1.1 ii libgraphicsmagick3 1.3.16-1.1 ii liblapack3 [liblapack.so.3] 3.4.1+dfsg-1 ii liboctave1 3.6.2-5 ii libpcre3 1:8.30-5 ii libqhull52009.1-3 ii libqrupdate1 1.1.1-1 ii libstdc++6 4.7.2-5 ii libumfpack5.4.0 1:3.4.0-3 ii octave-common3.6.2-5 ii texinfo 4.13a.dfsg.1-10 ii zlib1g 1:1.2.7.dfsg-13 Versions of packages octave recommends: ii gnuplot-x11 4.6.0-8 ii libatlas3-base 3.8.4-9 Versions of packages octave suggests: pn octave-doc none pn octave-htmldoc none ii octave-info 3.6.2-5 -- no debconf information -- To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: octave: sparse matrix n*2^16
On 29 April 2013 06:25, Miroslaw Kwasniak miroslaw.kwasn...@pwr.wroc.pl 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_limitsint::max(). Your matrix is simply too large to be indexed, and this breaks assumptions elsewhere in our code. - Jordi G. H. -- To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: octave: sparse matrix n*2^16
On Mon, Apr 29, 2013 at 6:10 AM, Jordi Gutiérrez Hermoso jord...@octave.org wrote: On 29 April 2013 06:25, Miroslaw Kwasniak miroslaw.kwasn...@pwr.wroc.pl 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_limitsint::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 -- Ed Meyer
Bug#706376: [Pkg-octave-devel] Bug#706376: octave: sparse matrix n*2^16
On 29 April 2013 12:40, Ed Meyer eem2...@gmail.com wrote: On Mon, Apr 29, 2013 at 6:10 AM, Jordi Gutiérrez Hermoso jord...@octave.org wrote: On 29 April 2013 06:25, Miroslaw Kwasniak miroslaw.kwasn...@pwr.wroc.pl 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_limitsint::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_limitsoctave_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. -- To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: octave: sparse matrix n*2^16
On Mon, Apr 29, 2013 at 9:47 AM, Jordi Gutiérrez Hermoso jord...@octave.org wrote: 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_limitsoctave_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. true for full matrices but sparse matrices are stored as three arrays and the nonzero and row index arrays are the only ones that need be limited. So you are saying that sparse matrices are treated as full in some places? -- Ed Meyer
Bug#706376: [Pkg-octave-devel] Bug#706376: octave: sparse matrix n*2^16
On 29 April 2013 13:21, Ed Meyer eem2...@gmail.com wrote: On Mon, Apr 29, 2013 at 9:47 AM, Jordi Gutiérrez Hermoso jord...@octave.org wrote: 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_limitsoctave_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. true for full matrices but sparse matrices are stored as three arrays and the nonzero and row index arrays are the only ones that need be limited. So you are saying that sparse matrices are treated as full in some places? No, the issue is that *all* indices are limited to octave_idx_type. We would have to do some sort of deep surgery to introduce a special indexing type for sparse matrices only. It seems like a lot of work for potentially little benefit. And yes, sparse matrices can be indexed by a single index instead of two, like any other matrix. Internally in Octave's source, the assumption that a single index of octave_idx_type is available is used throughout. - Jordi G. H. -- To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: octave: sparse matrix n*2^16
On 29 April 2013 14:00, Jordi Gutiérrez Hermoso jord...@octave.org wrote: And yes, sparse matrices can be indexed by a single index instead of two, like any other matrix. Internally in Octave's source, the assumption that a single index of octave_idx_type is available is used throughout. The specific case where this fails in this instance is octave_idx_type dim_vector::numel (), which is obtained simply by multiplying each of the dimensions of the matrix, even for sparse matrices (this is unlike nnz, so a workaround just for trace.m would be to use nnz instead of numel, but I think this would still leave some pretty broken sparse matrices lying around with other problems). We would have to change numel () to use some other type that can hold the result of a larger size, but this is a pretty fundamental function in Octave. The overall assumption is that you can linearly index up to numel (). Every place that calls this function would need to be checked to see what happens if we change its return type to be some special bigint. - Jordi G. H. -- To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: octave: sparse matrix n*2^16
On Mon, Apr 29, 2013 at 11:07 AM, Jordi Gutiérrez Hermoso jord...@octave.org wrote: On 29 April 2013 14:00, Jordi Gutiérrez Hermoso jord...@octave.org wrote: And yes, sparse matrices can be indexed by a single index instead of two, like any other matrix. Internally in Octave's source, the assumption that a single index of octave_idx_type is available is used throughout. The specific case where this fails in this instance is octave_idx_type dim_vector::numel (), which is obtained simply by multiplying each of the dimensions of the matrix, even for sparse matrices (this is unlike nnz, so a workaround just for trace.m would be to use nnz instead of numel, but I think this would still leave some pretty broken sparse matrices lying around with other problems). We would have to change numel () to use some other type that can hold the result of a larger size, but this is a pretty fundamental function in Octave. The overall assumption is that you can linearly index up to numel (). Every place that calls this function would need to be checked to see what happens if we change its return type to be some special bigint. - Jordi G. H. I'm not proposing using anything but octave_idx_type for indexing or changing the return type of numel() - I just question why numel() is used for sparse matrices. It should be irrelevant for anything but ccs2full(). Rather than restrict the size of sparse matrices I think it would make more sense to fix problems like this as they come up so that sparse storage is used as it was intended - to reduce storage op count. -- Ed Meyer
Bug#706376: [Pkg-octave-devel] Bug#706376: octave: sparse matrix n*2^16
On 29 April 2013 14:50, Ed Meyer eem2...@gmail.com wrote: I'm not proposing using anything but octave_idx_type for indexing or changing the return type of numel() - I just question why numel() is used for sparse matrices. It should be irrelevant for anything but ccs2full(). The numel function is just one example where this sparse matrix with big dimensions broke. Sparse matrices this large can still break in other ways. Furthermore, what do you propose to do if numel is called for sparse matrices, despite your suggestion to not call it for sparse matrices? A lot of code out there already assumes that you can treat a sparse matrix like any other matrix. I don't think numel is to blame. Another way I can think of would be that A(idx) would also break if idx is greater than the maximum index size. We would have to introduce special rules to handle that for sparse matrices of large dimensions. The problem isn't the storage size, it's the index size, which is why for other matrices Octave's error message says out of memory or index too big. I really don't see a way around this other than introducing a special index type for sparse matrices, and I don't see this as hugely useful. It also looks like a lot of boring work. But if you insist on doing this, don't let me discourage you. If you can figure out a consistent behaviour that doesn't break current code, you write the patch, and all tests pass, I'll happily apply it. - Jordi G. H. -- To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: octave: sparse matrix n*2^16
On 04/29/2013 03:12 PM, Jordi Gutiérrez Hermoso wrote: On 29 April 2013 14:50, Ed Meyereem2...@gmail.com wrote: I'm not proposing using anything but octave_idx_type for indexing or changing the return type of numel() - I just question why numel() is used for sparse matrices. It should be irrelevant for anything but ccs2full(). The numel function is just one example where this sparse matrix with big dimensions broke. Sparse matrices this large can still break in other ways. Furthermore, what do you propose to do if numel is called for sparse matrices, despite your suggestion to not call it for sparse matrices? A lot of code out there already assumes that you can treat a sparse matrix like any other matrix. I don't think numel is to blame. Another way I can think of would be that A(idx) would also break if idx is greater than the maximum index size. We would have to introduce special rules to handle that for sparse matrices of large dimensions. The problem isn't the storage size, it's the index size, which is why for other matrices Octave's error message says out of memory or index too big. I really don't see a way around this other than introducing a special index type for sparse matrices, and I don't see this as hugely useful. It also looks like a lot of boring work. But if you insist on doing this, don't let me discourage you. If you can figure out a consistent behaviour that doesn't break current code, you write the patch, and all tests pass, I'll happily apply it. - Jordi G. H. Will 64bit Octave automatically fix this? Michael
Bug#706376: [Pkg-octave-devel] Bug#706376: Bug#706376: octave: sparse matrix n*2^16
On 29 April 2013 16:46, Michael D. Godfrey michaeldgodf...@gmail.com wrote: On 04/29/2013 03:12 PM, Jordi Gutiérrez Hermoso wrote: On 29 April 2013 14:50, Ed Meyer eem2...@gmail.com wrote: I'm not proposing using anything but octave_idx_type for indexing or changing the return type of numel() - I just question why numel() is used for sparse matrices. It should be irrelevant for anything but ccs2full(). The numel function is just one example where this sparse matrix with big dimensions broke. Sparse matrices this large can still break in other ways. Furthermore, what do you propose to do if numel is called for sparse matrices, despite your suggestion to not call it for sparse matrices? A lot of code out there already assumes that you can treat a sparse matrix like any other matrix. I don't think numel is to blame. Will 64bit Octave automatically fix this? Yes, if you mean compiling Octave with 64-bit indexing. Or at least it will delay the problem until someone tries to create a matrix of size 2^64 by 2^64. - Jordi G. H. -- To UNSUBSCRIBE, email to debian-bugs-dist-requ...@lists.debian.org with a subject of unsubscribe. Trouble? Contact listmas...@lists.debian.org
Bug#706376: [Pkg-octave-devel] Bug#706376: octave: sparse matrix n*2^16
On Mon, Apr 29, 2013 at 12:12 PM, Jordi Gutiérrez Hermoso jord...@octave.org wrote: On 29 April 2013 14:50, Ed Meyer eem2...@gmail.com wrote: I'm not proposing using anything but octave_idx_type for indexing or changing the return type of numel() - I just question why numel() is used for sparse matrices. It should be irrelevant for anything but ccs2full(). The numel function is just one example where this sparse matrix with big dimensions broke. Sparse matrices this large can still break in other ways. Furthermore, what do you propose to do if numel is called for sparse matrices, despite your suggestion to not call it for sparse matrices? A lot of code out there already assumes that you can treat a sparse matrix like any other matrix. I don't think numel is to blame. I'm not saying numel() is to blame or should be changed, only that I see no reason to ever use numel when handling sparse matrices unless you are converting it to full in which case the current behavior is ok. Another way I can think of would be that A(idx) would also break if idx is greater than the maximum index size. We would have to introduce special rules to handle that for sparse matrices of large dimensions. Here again, why would you ever want A(idx) for a sparse matrix? The problem isn't the storage size, it's the index size, which is why for other matrices Octave's error message says out of memory or index too big. I really don't see a way around this other than introducing a special index type for sparse matrices, and I don't see this as hugely useful. It also looks like a lot of boring work. I agree that a special index type is the wrong approach; what I'm saying is that with sparse matrices you should never run into this problem in the first place if you don't try to treat them the same as full. Otherwise why have sparse matrices at all? But if you insist on doing this, don't let me discourage you. If you can figure out a consistent behaviour that doesn't break current code, you write the patch, and all tests pass, I'll happily apply it. this doesn't seem to be a burning issue so I'll shut up -- Ed Meyer
