CODEWITHSUNDEEP
pythonparallel-processingcython
ABZANMASTER
ABZANMASTER

Reputation: 105

Parallel programming in Cython

from numpy.linalg import inv
from numpy import array
def test (A,U,m):
    Lambda=array(list(map(lambda x:inv(U)@A[x]@U,(range(m*m)))))

This is a simple code that calculates an array product. How can I right this in cython with parallel programming in the loops with prange? I tried many times but in order to do so, I need to use nogil for prange. But inv() needs the gil. How can this be done efficiently in order to be faster than my original code?

Upvotes: 0

Views: 164

Answers (1)

DavidW
DavidW

Reputation: 30899

  1. You rarely actually want to compute the inverse. inv(U)@A[x] is better computed as np.linalg.solve(U, A[x]). That is almost certainly quicker and more numerically accurate. From a quick test on my PC you might expect a 40% speedup from that alone.

  2. You calculate the inverse inside a loop, but it does not depend on the looped variable. Consider taking the calculation of the inverse outside the loop and calculating it once! (or using solve...). You should fix this type of basic mistake before worrying about parallelizing it.

  3. The last time you asked this question you were pointed towards a different question. Essentially Scipy (not Numpy) comes with wrappers for BLAS and LAPACK. You can cimport these to get direct, GIL-free access to these low-level matrix algebra functions.

    I don't plan on giving a full worked example here (laziness...) but you can use dgesv for solve (assuming that both U and A[x] are double precision floating point matrices... you don't say) and dgemm for the matrix multiplication (again, assuming double-precision floating point matrices). If you insist on calculating the inverse then I think you need to go through its LU factorization.

    These functions are obviously considerable more complex than just using the Numpy functions. To get the pointers to pass to them you'll probably want to take a contiguous memoryview of your array cdef double[:,::1] Umview = U and then get the pointer of the 0,0 element (&U[0,0]).

Upvotes: 1

Related Questions