Indeed for matrix multiplication and many other L3 BLAS functions, we are
lucky however for linear solve function ?getrs unfortunately no avail.
On Thu, Nov 11, 2021 at 12:31 AM Benjamin Root <ben.v.r...@gmail.com> wrote:
> I have found that a bunch of lapack functions seem to have arguments for
> stating whether or not the given arrays are C or F ordered. Then you
> wouldn't need to worry about handling the layout yourself. For example, I
> have some C++ code like so:
>
> extern "C" {
>
> /**
> * Forward declaration for LAPACK's Fortran dgemm function to allow use in
> C/C++ code.
> *
> * This function is used for matrix multiplication between two arrays of
> doubles.
> *
> * For complete reference:
> http://www.netlib.org/lapack/explore-html/d1/d54/group__double__blas__level3_gaeda3cbd99c8fb834a60a6412878226e1.html
> */
> void dgemm_(const char* TRANSA, const char* TRANSB, const int* M, const
> int* N, const int* K,
> const double* ALPHA, const double* A, const int* LDA, const double* B,
> const int* LDB,
> const double* BETA, double* C, const int* LDC);
> }
>
> ...
>
> dgemm_("C", "C", &nLayers, &N, &nVariables, &alpha, matrices.IW->data(),
> &nVariables,
> inputs.data(), &N, &beta, intermediate.data(), &nLayers);
>
> (in this case, I was using boost multiarrays, but the basic idea is the
> same). IIRC, a bunch of other lapack functions had similar features.
>
> I hope this is helpful.
>
> Ben Root
>
>
>
> On Wed, Nov 10, 2021 at 6:02 PM Ilhan Polat <ilhanpo...@gmail.com> wrote:
>
>> I've asked this in Cython mailing list but probably I should also get
>> some feedback here too.
>>
>> I have the following function defined in Cython and using flat memory
>> pointers to hold n by n array data.
>>
>>
>> cdef some_C_layout_func(double[:, :, ::1] Am) nogil: # ... cdef double *
>> work1 = <double*>malloc(n*n*sizeof(double)) cdef double *work2 = <double
>> *>malloc(n*n*sizeof(double)) # ... # Lots of C-layout operations here #
>> ... dgetrs(<char*>'T', &n, &n, &work1[0], &n, &ipiv[0], &work2[0], &n, &
>> info ) dcopy(&n2, &work2[0], &int1, &Am[0, 0, 0], &int1) free(...)
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> Here, I have done everything in C layout with work1 and work2 but I have
>> to convert work2 into Fortran layout to be able to solve AX = B. A can be
>> transposed in Lapack internally via the flag 'T' so the only obstacle I
>> have now is to shuffle work2 which holds B transpose in the eyes of Fortran
>> since it is still in C layout.
>>
>> If I go naively and make loops to get one layout to the other that
>> actually spoils all the speed benefits from this Cythonization due to cache
>> misses. In fact 60% of the time is spent in that naive loop across the
>> whole function. Same goes for the copy_fortran() of memoryviews.
>>
>> I have measured the regular NumPy np.asfortranarray() and the
>> performance is quite good enough compared to the actual linear solve. Hence
>> whatever it is doing underneath I would like to reach out and do the same
>> possibly via the C-API. But my C knowledge basically failed me around this
>> line
>> https://github.com/numpy/numpy/blob/8dbd507fb6c854b362c26a0dd056cd04c9c10f25/numpy/core/src/multiarray/multiarraymodule.c#L1817
>>
>> I have found the SO post from
>> https://stackoverflow.com/questions/45143381/making-a-memoryview-c-contiguous-fortran-contiguous
>> but I am not sure if that is the canonical way to do it in newer Python
>> versions.
>>
>> Can anyone show me how to go about it without interacting with Python
>> objects?
>>
>> Best,
>> ilhan
>> _______________________________________________
>> NumPy-Discussion mailing list -- numpy-discussion@python.org
>> To unsubscribe send an email to numpy-discussion-le...@python.org
>> https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
>> Member address: ben.v.r...@gmail.com
>>
> _______________________________________________
> NumPy-Discussion mailing list -- numpy-discussion@python.org
> To unsubscribe send an email to numpy-discussion-le...@python.org
> https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
> Member address: ilhanpo...@gmail.com
>
_______________________________________________
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
Member address: arch...@mail-archive.com
