3

My problem is very similar to that posed here, and so I am trying to implement what was proposed in the answer there, based on this example of calling blas bundled with python. My code is only a slightly modified version of that example.pyx, here's my relevant code (with file f2pyptr.h unchanged) :

#example.pyx (single precision linear system solver via lapack)
import numpy as np
import scipy.linalg.lapack
import cython
cimport numpy as np

cdef extern from "/(pwd)/f2pyptr.h":
void *f2py_pointer(object) except NULL

ctypedef int sgesv_t(
int *N, int *NRHS,
float *A, int* LDA,
int* IPIV,
float *B,
int *LDB,
int *INFO)

cdef sgesv_t *sgesv = <sgesv_t*>f2py_pointer(scipy.linalg.lapack.sgesv._cpointer)

def myfunc_sgesv():
  cdef int N, NRHS, LDA, LDB, INFO
  N = 10
  NRHS = 1
  LDA = 10
  LDB = 10
  cdef np.ndarray[np.float32_t, ndim=2] A = np.ones((N,N), np.float32, order = "F")
  cdef np.ndarray[np.float32_t, ndim=2] B = np.ones((N, NRHS),np.float32, order = "F")
  cdef np.ndarray[np.int_t, ndim=1] IPIV = np.empty((N,), dtype = np.int)
  sgesv(&N, &NRHS, &A[0,0], &LDA, &IPIV, &B[0,0], &LDB, &INFO)*

def myfunc_sgesv():
  cdef int N, NRHS, LDA, LDB, INFO
  N = 10
  NRHS = 1
  LDA = 10
  LDB = 10

  cdef np.ndarray[np.float32_t, ndim=2] A = np.ones((N,N), np.float32, order = "F")
  cdef np.ndarray[np.float32_t, ndim=2] B = np.ones((N, NRHS),np.float32, order = "F")
  cdef np.ndarray[np.int_t, ndim=1] IPIV = np.empty((N,), dtype = np.int)

  sgesv(&N, &NRHS, &A[0,0], &LDA, &IPIV, &B[0,0], &LDB, &INFO)

With the standard setup.py file, I call python setup.py build_ext --inplace

which results in

Compiling example.pyx because it changed.
Cythonizing example.pyx

Error compiling Cython file:
------------------------------------------------------------
...

    cdef np.ndarray[np.float32_t, ndim=2] A = np.ones((N,N), np.float32, order = "F")
    cdef np.ndarray[np.float32_t, ndim=2] B = np.ones((N, NRHS),np.float32, order = "F")
    cdef np.ndarray[np.int_t, ndim=1] IPIV = np.empty((N,), dtype = np.int)

    sgesv(&N, &NRHS, &A[0,0], &LDA, &IPIV, &B[0,0], &LDB, &INFO)
                            ^
------------------------------------------------------------

example.pyx:103:33: Cannot take address of Python variable
Traceback (most recent call last):
File "setup.py", line 5, in <module>
ext_modules = cythonize("example.pyx")
File "/idiap/home/jnewling/.local/lib/python2.7/site-packages/Cython-0.20.1/Cython/Buil/Dependencies.py", line 785, in cythonize
cythonize_one(*args[1:]) File "/idiap/home/jnewling/.local/lib/python2.7/site-packages/Cython-0.20.1/Cython/Buil/Dependencies.py", line 902, in cythonize_one raise CompileError(None, pyx_file) Cython.Compiler.Errors.CompileError: example.pyx

I have the original example.pyx running fine, as well as a modified single precision version (i.e. sgemm), and so cannot see why my sgesv implementation fails.

Any specific diagnosis or suggestions of alternative approaches would be most welcome.

Community
  • 1
  • 1
newling
  • 624
  • 6
  • 10
  • You may want to take a look at this [repository](https://github.com/insertinterestingnamehere/cylinalg). You should be able to install it, then cimport from it. The only catch is that, in your setup files, you'll have to add the include directory where the header and cython `pxd` files are installed. – IanH Sep 23 '14 at 19:57

1 Answers1

3

Not sure if this answer is correct, but I will tell you what I tried as it may help you:

I got the error Cannot take address of Python variable at a different place when I compiled your code. I got the error at & IPIV. I then noticed it is an array so I changed it to:

sgesv(&N, &NRHS, &A[0,0], &LDA, &IPIV[0], &B[0,0], &LDB, &INFO). I then got the error Cannot assign type 'int_t *' to 'int *' . Then found the following answer: Passing numpy integer array to c code and changed the code to:

cdef np.ndarray[int, ndim=1, mode='c'] IPIV = np.empty((N,), dtype = ctypes.c_int) sgesv(&N, &NRHS, &A[0,0], &LDA, &IPIV[0], &B[0,0], &LDB, &INFO)

and now it compiles without any errors.

Yuvika
  • 5,624
  • 2
  • 16
  • 21
  • Perfect answer! Strange that the ^ was not under & IPIV, but if I had counted 33 characters it would have been under & IPIV and I would have (I like to believe) found the careless error. The second compile error looks nastier, thank you for saving me the time of diagnosing that too! – newling Sep 18 '14 at 18:03