Hi,
I suspect the current tests are not enought. You need to test all the
combination for the 3 inputs with thoses strides:
c-contiguous
f-contiguous
something else like strided.
Also, try with matrix with shape of 1 in each dimensions. Not all blas
libraries accept the strides that numpy use in that cases. Also, not all
blas version accept the same stuff, so if this isn't in the current
version, there will be probably some adjustment later on that side. What
blas do you use? I think ATLAS was one that was causing problem.
When we did this in Theano, it was more complicated then this diff... But
much of the code is boillerplate code.
Fred
On Thu, Nov 8, 2012 at 8:03 PM, Nicolas SCHEFFER wrote: Thanks Sebastien, didn't think of that. Well I went ahead and tried the change, and it's indeed straightforward. I've run some tests, among which:
nosetests numpy/numpy/core/tests/test_blasdot.py
and it looks ok. I'm assuming this is good news. I've copy-pasting the diff below, but I have that in my branch and can
create a PR if we agree on it.
I still cannot believe it's that easy (well this has been bugging me a
while... ;))
So I wouldn't mind waiting a day or two to see reactions on the list
before moving ahead. diff --git a/numpy/core/blasdot/_dotblas.c b/numpy/core/blasdot/_dotblas.c
index c73dd6a..2b4be7c 100644
--- a/numpy/core/blasdot/_dotblas.c
+++ b/numpy/core/blasdot/_dotblas.c
@@ -770,7 +770,7 @@ dotblas_matrixproduct(PyObject *NPY_UNUSED(dummy),
PyObject *args, PyObject* kwa
* using appropriate values of Order, Trans1, and Trans2.
*/ - if (!PyArray_ISCONTIGUOUS(ap2)) {
+ if (!PyArray_IS_C_CONTIGUOUS(ap2) &&
!PyArray_IS_F_CONTIGUOUS(ap2)) {
PyObject *new = PyArray_Copy(ap2); Py_DECREF(ap2);
@@ -779,7 +779,7 @@ dotblas_matrixproduct(PyObject *NPY_UNUSED(dummy),
PyObject *args, PyObject* kwa
goto fail;
}
}
- if (!PyArray_ISCONTIGUOUS(ap1)) {
+ if (!PyArray_IS_C_CONTIGUOUS(ap1) &&
!PyArray_IS_F_CONTIGUOUS(ap1)) {
PyObject *new = PyArray_Copy(ap1); Py_DECREF(ap1);
@@ -800,6 +800,19 @@ dotblas_matrixproduct(PyObject
*NPY_UNUSED(dummy), PyObject *args, PyObject* kwa
lda = (PyArray_DIM(ap1, 1) > 1 ? PyArray_DIM(ap1, 1) : 1);
ldb = (PyArray_DIM(ap2, 1) > 1 ? PyArray_DIM(ap2, 1) : 1);
ldc = (PyArray_DIM(ret, 1) > 1 ? PyArray_DIM(ret, 1) : 1);
+
+ /*
+ * Avoid temporary copies for arrays in Fortran order
+ */
+ if (PyArray_IS_F_CONTIGUOUS(ap1)) {
+ Trans1 = CblasTrans;
+ lda = (PyArray_DIM(ap1, 0) > 1 ? PyArray_DIM(ap1, 0) : 1);
+ }
+ if (PyArray_IS_F_CONTIGUOUS(ap2)) {
+ Trans2 = CblasTrans;
+ ldb = (PyArray_DIM(ap2, 0) > 1 ? PyArray_DIM(ap2, 0) : 1);
+ }
+
if (typenum == NPY_DOUBLE) {
cblas_dgemm(Order, Trans1, Trans2,
L, N, M, On Fri, 2012-11-09 at 00:24 +0100, Sebastian Berg wrote: Hey, On Thu, 2012-11-08 at 14:44 -0800, Nicolas SCHEFFER wrote: Well, hinted by what Fabien said, I looked at the C level dot
function.
Quite verbose! But starting line 757, we can see that it shouldn't be too much work
to fix that bug (well there is even a comment there that states just
that) https://github.com/numpy/numpy/blob/master/numpy/core/blasdot/_dotblas.c#L75... I now think that should be the cleanest. This would only work for gemm though. I don't know what the benefit is
for gemv for instance, but we should make that kind of changes
everywhere we can.
The evil PyArray_Copy is there twice and that's what we want to get
rid of. I'm not sure, but it looks to me that removing the copy and doing the
following would do the work:
Order = CblasRowMajor;
Trans1 = CblasNoTrans;
Trans2 = CblasNoTrans;
if (!PyArray_ISCONTIGUOUS(ap1)) {
Trans1 = CblasTrans;
}
if (!PyArray_ISCONTIGUOUS(ap2)) {
Trans2 = CblasTrans;
}
might be too easy to be true. Sounds nice, though don't forget that the array may also be neither C-
or F-Contiguous, in which case you need a copy in any case. So it would
probably be more like: if (PyArray_IS_C_CONTIGUOUS(ap1)) {
Trans1 = CblasNoTrans;
}
else if (PyArray_IS_F_CONTIGUOUS(ap1)) {
Trans1 = CblasTrans;
}
else {
Trans1 = CblasNoTrans;
PyObject *new = PyArray_Copy(ap1);
Py_DECREF(ap1);
ap1 = (PyArrayObject *)new;
} Well, of course I forgot error checking there, and maybe you need to set
some of the other parameters differently, but it looks like its probably
that easy, and I am sure everyone will welcome a PR with such changes. Regards, Sebastian On Thu, Nov 8, 2012 at 12:06 PM, Nicolas SCHEFFER
I've made the necessary changes to get the proper order for the output array. Also, a pass of pep8 and some tests (fixmes are in failing tests)
http://pastebin.com/M8TfbURi -n On Thu, Nov 8, 2012 at 11:38 AM, Nicolas SCHEFFER
Thanks for all the responses folks. This is indeed a nice problem
to solve. Few points:
I. Change the order from 'F' to 'C': I'll look into it.
II. Integration with scipy / numpy: opinions are diverging here.
Let's wait a bit to get more responses on what people think.
One thing though: I'd need the same functionality as
get_blas_funcs in numpy.
Since numpy does not require lapack, what functions can I get?
III. Complex arrays
I unfortunately don't have enough knowledge here. If someone could
propose a fix, that'd be great.
IV. C
Writing this in C sounds like a good idea. I'm not sure I'd be the
right person to this though.
V. Patch in numpy
I'd love to do that and learn to do it as a byproduct.
Let's make sure we agree this can go in numpy first and that all
FIXME
can be fixed.
Although I guess we can resolve fixmes using git. Let me know how you'd like to proceed, Thanks! FIXMEs:
- Fix for ndim != 2
- Fix for dtype == np.complex*
- Fix order of output array On Thu, Nov 8, 2012 at 9:42 AM, Frédéric Bastien On Thu, Nov 8, 2012 at 3:58 PM, Sebastian Berg
>
> HTH
>
> Fred
>
>
>
> On Thu, Nov 8, 2012 at 12:29 PM, Anthony Scopatz <
scopatz@gmail.com> wrote:
>>
>> On Thu, Nov 8, 2012 at 7:06 AM, David Cournapeau <
cournape@gmail.com>
>> wrote:
>>>
>>> On Thu, Nov 8, 2012 at 12:12 PM, Dag Sverre Seljebotn
>>> NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion _______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion _______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion