[Numpy-discussion] Scipy dot

Sebastian Berg sebastian@sipsolutions....
Thu Nov 8 17:58:29 CST 2012


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




More information about the NumPy-Discussion mailing list