python - Wrapping a LAPACKE function using Cython -

I am trying to wrap the LAPACK function (a solver for the triaday system of equations) using the statement.

I came across, but since dgtsv is not one of the LAPACK functions wrapped in scipy.linalg , I do not think I would like this particular I can use the approach. Instead I am trying to follow

Here is my content of lapacke.pxd file:

  ctypedef integer lapack_int cdef "lapacke .h "Excluded from nogil: integer LAPACK_ROW_MAJOR integer LAPACK_COL_MAJOR lapack_int LAPACKE_dgtsv (integer matrix_order, lapack_int n, lapack_int nrhs, double * dl, double * d, double * du, double * b, lapack_int ldb)   
 ... here is my thin Cython wrapper in  _solvers.pyx :  
  #! python cimport cython lapacke cimport * cpdef TDMA_lapacke (double [:: 1] dl, double [:: 1] d, double [:: 1] du, double [:, :: 1] b): cd: lapak_int n = D. Espe [0] lepac_INT NRHS = biship [1] lepac_ind LDB = bishap [0] double * DL = & amp; DL [0] double * d = & amp; D [0] double * du = & amp; DU [0] double * b = and b [0, 0] lapack_int information info = LAPACKE_dgtsv (LAPACK_ROW_MAJOR, n, nrhs, dl, d, du, b, ldb) information back   
 ... and there is a python wrapper and test script:  
  imported from CP imports as import numpy imported from SPM simulations, DSDF DRF Trisolve_Lapac (DL, D, Do, B, Inplay = Fills): If (DL.spa [0]! = Du shape [0] or DL. Size [0]! = D. size [0] - 1 or b .shape! = D.shape): ValueError ( Increase 'invalid diagonal shape') if b. ndim == 1: # B (LDB, NRHS) B = B [, none] A copy of # D and B. If we are not sorting in place, then there is no place: d = d.copy ( ) B = b.copy () # This can also be copied if arrays are typed incorrectly / noncontiguous DL, D, du, B = (np.ascontiguousarray v (v, dtype = np.float64) v (Dl, d, du, b)) # b Now the solution information = _solvers.TDMA_lapacke (dl, d, du, b) print information return b.ravel () def test_trisolve (n = 20000): dl = np.random .randn (n - 1) d = np.random.randn (n) du = np.random.randn (n - 1) M = sparse.diags ((DL, D, DU), (-1, 0, 1 ), format = 'csc') x = Np.random.randn (n) b = M.dot (x) x_hat = Trisolve_lapake (DL, D, Do, B) print "|| x-x_hat || = ", Np.linalg.norm (x - x_hat)   
 Unfortunately, for calling  test_trisolve , just call  _solvers.TDMA_lapacke  Call segfaults I am absolutely sure that my  setup.py  is correct -  ldd _solvers.so  shows that the correct sharing of  _solvers.so   
   << >>  n  for small values ​​I do not immediately get segfaults, but the result of me nonsense ( = x - x_hat ||  should be very close to 0):  
  in [28]: test_trisolve2.test_trisolv e (10) 0 || X - x_hat || = 6.232025763 in 96 [29]: test_trisolve2.test_trisolve (10) -7 || x-x_hat || = 3.88623414288 in [30]: test_trisolve2.test_trisolve (10 ) 0} X = x_hat = 2.601 9 0676562 in [31]: test_trisolve2.test_trisolve (10) 0 || X-x_hat || = 3.86631743386 in [32]: test_trisolve2.test_trisolve (10) partition fault < / code>  
 usually returns with the  LAPACKE_dgtsv  code  0  (which shows success), but sometimes I'll see the  - 7, which means that there was an illegal value of Logic 7 ( b ). What's happening is that the first value of  b  is actually being modified in place. If I keep calling on  test_trisolve , I will eventually hit a segfault, whenever  n  is smaller.