On Thu, Nov 8, 2012 at 8:03 PM, Nicolas SCHEFFER <scheffer.nicolas@gmail.com> 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 Thu, Nov 8, 2012 at 3:58 PM, Sebastian Berg
<sebastian@sipsolutions.net> wrote:
> 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#L757
>> > 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
>> > <scheffer.nicolas@gmail.com> wrote:
>> > > 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
>> > > <scheffer.nicolas@gmail.com> wrote:
>> > >> 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 <nouiz@nouiz.org> wrote:
>> > >>> Hi,
>> > >>>
>> > >>> I also think it should go into numpy.dot and that the output order should
>> > >>> not be changed.
>> > >>>
>> > >>> A new point, what about the additional overhead for small ndarray? To remove
>> > >>> this, I would suggest to put this code into the C function that do the
>> > >>> actual work (at least, from memory it is a c function, not a python one).
>> > >>>
>> > >>> 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
>> > >>>>> <d.s.seljebotn@astro.uio.no> wrote:
>> > >>>>> > On 11/08/2012 01:07 PM, Gael Varoquaux wrote:
>> > >>>>> >> On Thu, Nov 08, 2012 at 11:28:21AM +0000, Nathaniel Smith wrote:
>> > >>>>> >>> I think everyone would be very happy to see numpy.dot modified to do
>> > >>>>> >>> this automatically. But adding a scipy.dot IMHO would be fixing
>> > >>>>> >>> things
>> > >>>>> >>> in the wrong place and just create extra confusion.
>> > >>>>> >>
>> > >>>>> >> I am not sure I agree: numpy is often compiled without lapack support,
>> > >>>>> >> as
>> > >>>>> >> it is not necessary. On the other hand scipy is always compiled with
>> > >>>>> >> lapack. Thus this makes more sens in scipy.
>> > >>>>> >
>> > >>>>> > Well, numpy.dot already contains multiple fallback cases for when it is
>> > >>>>> > compiled with BLAS and not. So I'm +1 on just making this an
>> > >>>>> > improvement
>> > >>>>> > on numpy.dot. I don't think there's a time when you would not want to
>> > >>>>> > use this (provided the output order issue is fixed), and it doesn't
>> > >>>>> > make
>> > >>>>> > sense to not have old codes take advantage of the speed improvement.
>> > >>>>>
>> > >>>>> Indeed, there is no reason not to make this available in NumPy.
>> > >>>>>
>> > >>>>> Nicolas, can you prepare a patch for numpy ?
>> > >>>>
>> > >>>>
>> > >>>> +1, I agree, this should be a fix in numpy, not scipy.
>> > >>>>
>> > >>>> Be Well
>> > >>>> Anthony
>> > >>>>
>> > >>>>>
>> > >>>>>
>> > >>>>> David
>> > >>>>> _______________________________________________
>> > >>>>> 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
>> >
>>
>>
>> _______________________________________________
>> 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