- c - 在位数组中找到第一个零
- linux - Unix 显示有关匹配两种模式之一的文件的信息
- 正则表达式替换多个文件
- linux - 隐藏来自 xtrace 的命令
我正在尝试包装 LAPACK 函数 dgtsv
(三对角方程组的求解器)使用 Cython。
我遇到了 this previous answer ,但由于 dgtsv
不是包装在 scipy.linalg
中的 LAPACK 函数之一,我认为我无法使用这种特殊方法。相反,我一直在尝试关注 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 包装器:
#!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 包装器和测试脚本:
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
被链接到正确的共享库在运行时。
我不太确定如何从这里开始 - 有什么想法吗?
简要更新:
对于较小的 n
值,我往往不会立即得到段错误,但我确实得到了无意义的结果(||x - x_hat|| 应该非常接近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
很小,我最终也会遇到段错误。
最佳答案
好吧,我最终弄明白了 - 看来我误解了在这种情况下行和列主要指的是什么。
由于 C 连续数组遵循行优先顺序,我假设我应该将 LAPACK_ROW_MAJOR
指定为 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
这对我来说似乎很违反直觉 - 谁能解释为什么会这样?
关于python - 使用 Cython 包装 LAPACKE 函数,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/23200085/
c 不做边界检查。那么cython是如何检查是否编译成c的呢? %%cython --annotate cimport cython @cython.boundscheck(True) cpdef m
可以直接声明用于 Cython 构造函数? 据我了解,这是可能的: # Cython cdef int[3] li = [1, 2, 3] # C++ int[3] li = {1, 2, 3} 但
所以,如果你有一个头文件。 %%file test.h struct mystruct{ int i; int j; }; 然后你将它包装在 Cython 中: cdef extern fr
我正在构建一个独立于平台的 cython 项目,我想根据正在使用的编译器传递编译器参数。我可以猜测基于平台的编译器,或者假设它与用于 Python 的编译器相同,但不能保证匹配。通常我注入(injec
我使用诗歌构建我的 cython 包。我在所有函数和类中都有 NumPy 风格的文档字符串。我现在要做的是添加 Sphinx 自动文档并发布在 Read the Docs。 我已阅读此主题 How d
赛通 libcpp模块包含 priority_queue 的模板,这很好,除了一件事:我不能通过自定义比较器(或者,至少,我不知道如何)。 我需要这个,因为我需要 priority_queue做一个a
以下代码定义了一个简单的 Cython 函数(为方便起见,使用 Ipython 魔法)。 %load_ext cython %%cython def f(float x, float y=2):
我正在尝试使用 cython 进行复数计算。在示例代码中,我想计算复数的复指数函数。问题是我不知道如何将我的整数乘以虚数单位。python的虚数单位1.0j乘以cython执行时报错。 这是我的代码:
在这里停留在一些基本的 Cython 上 - 在 Cython 中定义字符串数组的规范且有效的方法是什么? 具体来说,我想定义一个定长常量数组char . (请注意,此时我不想引入 NumPy。) 在
是否有可能,如果是,如何确定 Cython 中整数数据类型的大小(以位为单位)? 我正在尝试做这样的事情,以获得整数大小: cdef WORD_BITS = 0 IF sizeof(unsigned
我只是想打印 cython 变量的地址,但我无法绕过错误消息: cdef int myvar print &myvar 抛出 Cannot convert 'int *' to Python obje
我有一个 C 头文件,它在宏中定义了一个函数。我需要从 Cython 调用它。有没有办法在 Cython 中使用宏并使其完全扩展?我已经有了 C 类型的参数。 我尝试像使用函数一样使用 cdef,我认
令人惊讶的是,我似乎找不到通过名称获取结构体元素的单个示例(无论是在网络上还是在 cython 示例中)。 所以我收到了一个指向 C 函数结构体的指针,并且想要一一访问这些元素并将它们重新打包到 py
我尝试围绕 C++ 库编写一个 Cython 包装器 http://primesieve.org/ 它包装了一个函数count。到目前为止,它可以正确安装 python setup.py instal
我正在尝试将 cython 模块 data.pyx 导入另一个 cython 模块 user.pyx。一切都编译得很好,但是当我尝试在 python 模块中调用 user.pyx 时,我收到错误“Im
更新:内存 View 获胜。Cython 使用类型化内存 View :0.0253449 特别感谢 lothario,他指出了几个关键的变化。 荒谬。当然现在的问题是,似乎不能对它们做太多算术(加法和
我有一个使用 memoryview 数组的 cython 模块,即... double[:,:] foo 我想使用多处理并行运行这个模块。但是我得到了错误: PicklingError: Can't
我正在尝试使用 Cython 加速 PEP 484 类型的 python 脚本。我想保持一些语义和可读性。 之前,我有一个 Flags = int def difference(f1: Flags,
这个问题已经有答案了: Collapse multiple submodules to one Cython extension (5 个回答) 已关闭 3 年前。 我在一个包中有多个 .py 文件
我已经能够在我的 .pyx 脚本上使用 cython 在 linux 上创建一个 .so 文件。我也可以成功地在我的 python 解释器上进行导入。 我的问题是如何在不使用 cython 的情况下将
我是一名优秀的程序员,十分优秀!