lapack_lite should use 64-bit integer indices

[njsmith]

Currently lapack-lite uses int for indices, so large arrays cause overflow and things crash (e.g. #5898). There's a separate need to handle this issue when using real BLAS/LAPACK (which may or may not handle arrays with >2**31 elements, depending on the vendor), but when falling back on our built-in code there's no reason we shouldn't get this right.

[argriffing]

I remember a mailing list thread suggesting the possibility of blis, a BSD-licensed BLAS-like framework. I just checked its config notes, and it explicitly allows separate internal vs. interface integer sizes as documented in the 'BLAS compatibility layer' section in
https://code.google.com/p/blis/wiki/ConfigurationHowTo
In particular, you could have

#define BLIS_BLAS2BLIS_INT_TYPE_SIZE 32
#define BLIS_INT_TYPE_SIZE 64

If I understand correctly, the current numpy code doesn't have this flexibility. On the other hand I'm not sure if 'blis' and its related projects include enough of the BLAS/LAPACK functionality yet, or if its Windows support is mature enough, or if enough of its plumbing is otherwise compatible with numpy.

[njsmith]

Blis is only a replacement for blas, not for lapack. But fortunately, the
same people also have a pure-C lapack replacement called flame, which is
feature complete at least. So yeah, this is one possible direction
direction could go. They're at UT Austin so we should talk to them at
SciPy. (And I'm in contact contact them in the next email thread over...)
That's a somewhat more ambitious and complicated strategy than just doing a
search-replace on our existing fallback code though :-).
On May 22, 2015 8:45 AM, "argriffing" notifications@github.com wrote:

I remember a mailing list thread
http://mail.scipy.org/pipermail/numpy-discussion/2014-April/069868.html
suggesting the possibility of blis https://code.google.com/p/blis/, a
BSD-licensed BLAS-like framework. I just checked its config notes, and it
explicitly allows separate internal vs. interface integer sizes as
documented in the 'BLAS compatibility layer' section in
https://code.google.com/p/blis/wiki/ConfigurationHowTo
In particular, you could have

#define BLIS_BLAS2BLIS_INT_TYPE_SIZE 32
#define BLIS_INT_TYPE_SIZE 64

If I understand correctly, the current numpy code doesn't have this
flexibility. On the other hand I'm not sure if 'blis' and its related
projects include enough of the BLAS/LAPACK functionality yet, or if its
Windows support is mature enough, or if enough of its plumbing is otherwise
compatible with numpy.

—
Reply to this email directly or view it on GitHub
#5906 (comment).

[argriffing]

just doing a search-replace on our existing fallback code though :-)

Which file should be search and replaced? I think the C file is somehow autogenerated. Also, both the interface and internal variables use type integer, so it would not be possible to just re-typedef integer or to search-replace integer in the C file.

[njsmith]

I think (say, 70% confidence, so don't trust me blindly :-)), that: (a)
those files were were autogenerated once long ago via means that may or may
not be reproducible now, and which may or may not have been hand-tweaked
since, so hand tweaking them now is not necessarily an issue. (b) The files
are used entirely internally inside numpy, so we can change the interface
however we like.
On May 22, 2015 10:52 AM, "argriffing" notifications@github.com wrote:

just doing a search-replace on our existing fallback code though :-)

Which file should be search and replaced? I think the C file is somehow
autogenerated. Also, both the interface and internal variables use type
integer, so it would not be possible to just re-typedef integer or to
search-replace integer in the C file.

—
Reply to this email directly or view it on GitHub
#5906 (comment).

[argriffing]

we can change the interface however we like.

I think it would make sense for the interface to remain 32 bit, for compatibility with the real BLASes.

[njsmith]

What compatibility issues do you foresee? It's purely an internal
interface, and having a 64 bit version gives us something to fall back on
when real blas/lapack don't work. Bug-for-bug compatibility isn't really
a goal here.
On May 22, 2015 11:06 AM, "argriffing" notifications@github.com wrote:

we can change the interface however we like.

I think it would make sense for the interface to remain 32 bit, for
compatibility with the real BLASes.

—
Reply to this email directly or view it on GitHub
#5906 (comment).

[argriffing]

@njsmith The 32-bitness of the interface is not the main cause of the bugs. For example, https://github.com/numpy/numpy/blob/master/numpy/linalg/umath_linalg.c.src uses int in a bunch of places, including in the interface. I'm suggesting to not change these types. On the other hand, https://raw.githubusercontent.com/numpy/numpy/master/numpy/linalg/lapack_lite/dlapack_lite.c has a bunch of integer declarations, some of which are for the interface and some of which are for internal indexing. I think the interface types should remain 32-bit and the integer types used for indexing in that file should be bumped to size_t.

[njsmith]

What happens when someone has an array with more than 2**31 elements along
one axis?
On May 22, 2015 11:20 AM, "argriffing" notifications@github.com wrote:

@njsmith https://github.com/njsmith The 32-bitness of the interface is
not the main cause of the bugs. For example,
https://github.com/numpy/numpy/blob/master/numpy/linalg/umath_linalg.c.src
uses int in a bunch of places, including in the interface. I'm suggesting
to not change these types. On the other hand,
https://raw.githubusercontent.com/numpy/numpy/master/numpy/linalg/lapack_lite/dlapack_lite.c
has a bunch of integer declarations, some of which are for the interface
and some of which are for internal indexing. I think the interface types
should remain 32-bit and the integer types used for indexing in that file
should be bumped to size_t.

—
Reply to this email directly or view it on GitHub
#5906 (comment).

[argriffing]

What happens when someone has an array with more than 2**31 elements along one axis?

Then everything is ruined. I think this is the case even for the default settings of the 'good' real BLASes but I'm not completely sure.

[njsmith]

It is the case for real blas (modulo non-portable extensions), but that
doesn't mean it has to be the case for numpy.
On May 22, 2015 12:03 PM, "argriffing" notifications@github.com wrote:

What happens when someone has an array with more than 2**31 elements along
one axis?

Then everything is ruined. I think this is the case even for the default
settings of the 'good' real BLASes but I'm not completely sure.

—
Reply to this email directly or view it on GitHub
#5906 (comment).

[argriffing]

I think it would make sense to fix these problems one at a time. For example, the first improvement would get the numpy blas to the point where it does not segfault when real blas does not segfault, then later maybe add the second improvement that lets the numpy blas work with larger arrays than the real blas can handle.

[pv]

The fix is anyway doable in a mechanical way; just replace "integer n -> interace_int interface_n" in the argument lists, add "integer n = interface_n" declarations to the body, and finally change the default integer type to ssize_t.

[pv]

... except it isn't because it's full of integer pointers :( Lots of work if done manually, it seems.

[argriffing]

Lots of work if done manually, it seems.

I'm vaguely aware of research projects that try to automate this kind of thing, but I don't know any specifics. Also, some time around 2009-2010 the netlib lapack switched from the f2c CLAPACK (which I think we are using in numpy) to including C implementations in their main lapack distribution, if I understand their webpage correctly.

[njsmith]

According to the README in lapack_lite/, we are using an independent
translation sorta similar to but different from CLAPACK. There's also code
in there to do the translation and automatically postprocess the code in
various ways.
On May 22, 2015 12:53 PM, "argriffing" notifications@github.com wrote:

Lots of work if done manually, it seems.

I'm vaguely aware of research projects that try to automate this kind of
thing, but I don't know any specifics. Also, some time around 2009-2010 the
netlib lapack switched from the f2c CLAPACK (which I think we are using in
numpy) to including C implementations in their main lapack distribution, if
I understand their webpage correctly.

—
Reply to this email directly or view it on GitHub
#5906 (comment).