使用Cython包装LAPACKE函数 [英] Wrapping a LAPACKE function using Cython
问题描述
我正在尝试包装LAPACK函数
I'm trying to wrap the LAPACK function dgtsv (a solver for tridiagonal systems of equations) using Cython.
我遇到了此先前的答案，但是由于dgtsv不是包装在其中的LAPACK函数之一scipy.linalg我认为我不能使用这种特殊方法.相反，我一直在尝试遵循此示例.
I came across this previous answer, but since dgtsv is not one of the LAPACK functions that are wrapped in scipy.linalg I don't think I can use this particular approach. Instead I've been trying to follow this example.
这是我的lapacke.pxd文件的内容:
ctypedef int lapack_int
cdef extern from "lapacke.h" nogil:
int LAPACK_ROW_MAJOR
int LAPACK_COL_MAJOR
lapack_int LAPACKE_dgtsv(int matrix_order,
lapack_int n,
lapack_int nrhs,
double * dl,
double * d,
double * du,
double * b,
lapack_int ldb)
...这是我在_solvers.pyx中的薄Cython包装纸:
...here's my thin Cython wrapper in _solvers.pyx:
#!python
cimport cython
from lapacke cimport *
cpdef TDMA_lapacke(double[::1] DL, double[::1] D, double[::1] DU,
double[:, ::1] B):
cdef:
lapack_int n = D.shape[0]
lapack_int nrhs = B.shape[1]
lapack_int ldb = B.shape[0]
double * dl = &DL[0]
double * d = &D[0]
double * du = &DU[0]
double * b = &B[0, 0]
lapack_int info
info = LAPACKE_dgtsv(LAPACK_ROW_MAJOR, n, nrhs, dl, d, du, b, ldb)
return info
...这是一个Python包装器和测试脚本:
...and here's a Python wrapper and test script:
import numpy as np
from scipy import sparse
from cymodules import _solvers
def trisolve_lapacke(dl, d, du, b, inplace=False):
if (dl.shape[0] != du.shape[0] or dl.shape[0] != d.shape[0] - 1
or b.shape != d.shape):
raise ValueError('Invalid diagonal shapes')
if b.ndim == 1:
# b is (LDB, NRHS)
b = b[:, None]
# be sure to force a copy of d and b if we're not solving in place
if not inplace:
d = d.copy()
b = b.copy()
# this may also force copies if arrays are improperly typed/noncontiguous
dl, d, du, b = (np.ascontiguousarray(v, dtype=np.float64)
for v in (dl, d, du, b))
# b will now be modified in place to contain the solution
info = _solvers.TDMA_lapacke(dl, d, du, b)
print info
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_lapacke(dl, d, du, b)
print "||x - x_hat|| = ", np.linalg.norm(x - x_hat)
不幸的是，test_trisolve只是在调用_solvers.TDMA_lapacke时出现段错误.
我很确定我的setup.py是正确的-ldd _solvers.so显示_solvers.so在运行时已链接到正确的共享库.
Unfortunately, test_trisolve just segfaults on the call to _solvers.TDMA_lapacke.
I'm pretty sure my setup.py is correct - ldd _solvers.so shows that _solvers.so is being linked to the correct shared libraries at runtime.
我不太确定如何从这里开始-有什么想法吗?
I'm not really sure how to proceed from here - any ideas?
简要更新:
对于较小的n值，我往往不会立即得到段错误，但我确实会得到废话的结果( || x-x_hat || 应该非常接近0):
for smaller values of n I tend not to get segfaults immediately, but I do get nonsense results (||x - x_hat|| ought to be very close to 0):
In [28]: test_trisolve2.test_trisolve(10)
0
||x - x_hat|| = 6.23202576396
In [29]: test_trisolve2.test_trisolve(10)
-7
||x - x_hat|| = 3.88623414288
In [30]: test_trisolve2.test_trisolve(10)
0
||x - x_hat|| = 2.60190676562
In [31]: test_trisolve2.test_trisolve(10)
0
||x - x_hat|| = 3.86631743386
In [32]: test_trisolve2.test_trisolve(10)
Segmentation fault
通常，LAPACKE_dgtsv返回代码为0(应该表示成功)，但是偶尔我会收到-7，这意味着参数7(b)具有非法值.实际情况是，实际上仅修改了b的第一个值.如果我继续拨打test_trisolve，即使n很小，我最终也会遇到段错误.
Usually LAPACKE_dgtsv returns with code 0 (which should indicate success), but occasionally I get -7, which means that argument 7 (b) had an illegal value. What's happening is that only the first value of b is actually being modified in place. If I keep on calling test_trisolve I will eventually hit a segfault even when n is small.
推荐答案
好的，我最终弄清楚了-似乎我误解了在这种情况下行和列主要指的是什么.
OK, I figured it out eventually - it seems I've misunderstood what row- and column-major refer to in this case.
由于C连续数组遵循行优先顺序，因此我假设我应该将LAPACK_ROW_MAJOR指定为LAPACKE_dgtsv的第一个参数.
Since C-contiguous arrays follow row-major order, I assumed that I ought to specify LAPACK_ROW_MAJOR as the first argument to LAPACKE_dgtsv.
事实上，如果我改变了
info = LAPACKE_dgtsv(LAPACK_ROW_MAJOR, ...)
到
info = LAPACKE_dgtsv(LAPACK_COL_MAJOR, ...)
然后我的函数起作用:
test_trisolve2.test_trisolve()
0
||x - x_hat|| = 6.67064747632e-12
这对我来说似乎很违反直觉-有人可以解释为什么会这样吗?
This seems pretty counter-intuitive to me - can anyone explain why this is the case?
这篇关于使用Cython包装LAPACKE函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
