[Numpy-discussion] Scipy dot

[Sebastian Berg]

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;
}

Regards,

Sebastian

> 
> 
> On Thu, Nov 8, 2012 at 12:06 PM, Nicolas SCHEFFER
> <scheffer.nicolas at 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 at 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 at 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 at gmail.com> wrote:
> >>>>
> >>>> On Thu, Nov 8, 2012 at 7:06 AM, David Cournapeau <cournape at gmail.com>
> >>>> wrote:
> >>>>>
> >>>>> On Thu, Nov 8, 2012 at 12:12 PM, Dag Sverre Seljebotn
> >>>>> <d.s.seljebotn at 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 at scipy.org
> >>>>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
> >>>>
> >>>>
> >>>>
> >>>> _______________________________________________
> >>>> NumPy-Discussion mailing list
> >>>> NumPy-Discussion at scipy.org
> >>>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
> >>>>
> >>>
> >>>
> >>> _______________________________________________
> >>> NumPy-Discussion mailing list
> >>> NumPy-Discussion at scipy.org
> >>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
> >>>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>