- 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/
如果我有一个 1,000 x 1,000 的方阵,Lapack 可以计算这个矩阵的特征向量和特征值吗?如果可以,需要多长时间?对于 10,000 x 10,000 矩阵甚至 1,000,000 x 1
我正在使用 LAPACK 的 zheev(在英特尔 MKL 中)。我得到了 int INFO=99。我一直在互联网上搜索这对应的内容,但找不到包含所有整数错误代码及其含义列表的文档。 有没有人有指向
通过以下方式安装 lapack 后: yum install lapack lapack-devel 并重新启动 httpd 服务我仍然得到 Fatal error: Class 'Lapack' n
Lapack 很可能没有任何计算行列式的例程。但是,我们可以使用 LU、QR 或 SVD 分解来计算它。我更喜欢使用 LU 分解。现在 lapack 使用一些 dgetrf 子程序将矩阵 A 分解为带
根据官方用户指南,sgelsd用于解决最小二乘问题 min_x || b - Ax ||_2 并允许矩阵 A 为矩形且秩亏。并且根据sgelsd源码中的接口(interface)描述,b作为输入输出参
我用 Eigen 和 实现了一段代码。我希望 Eigen 使用 BLAS 和 LAPACK 。 我看过here ,这是可能的,但我不知道如何或将这些值/指令放在代码中的位置。 我必须在某处指定值 EI
简介:我用 C++ 开发了一个应用程序,它在 Windows 上使用 LAPACK(LAPACKE) 和 MPI。在 Windows 中工作正常(编译和链接通过 Code::Blocks IDE 处理
我最近从 Linux 切换到 Mac OS。我需要 BLAS 和 LAPACK 来做一些计算。通过查看 BLAS 的维基百科,我了解到这两个库已经在 Mac OS 中实现。不过,据说 Apple's
我基于以下链接为我的 Visual Studio 2008 构建了 LAPACKE 的 DLL 和库: http://icl.cs.utk.edu/lapack-for-windows/lapack/
我很好奇用于在 MATLAB 中计算 SVD 的 DGESVD 函数。据我从 Gene H. Golub 和 Charles F. Van Loan 的“矩阵计算”中可以看出,使用了两种可能的双对角化
我是数值线性代数的新手,我刚刚开始使用 LAPACK 和 BLAS。 是否有可以在打包存储和完整存储之间复制/转换对称矩阵的例程? 我找到了 dtrttp ,我可以用它来将 double 全对称矩阵转
根据我的理解,需要进行分解/因式分解(LU、QR、Cholesky 等),然后基于因式分解进行矩阵逆计算。还有其他方法可以解决这个问题吗(我试图弄清楚我是否可以坚持使用 CULAtools 试用版中免
我正在尝试编写一个函数,该函数可以为代表性不足的方程组生成单一解(例如,描述该系统的矩阵宽大于高)。为了做到这一点,我一直在 LAPACK 文档中寻找一种将矩阵行归约到它的归约梯队形式的方法,类似于
我是使用 LAPACK 例程的新手,所以我对它们不是很了解,我想在并行化循环(openmp)中使用它们。 我使用 Ubuntu 14.04LTS 并使用我的包管理器安装了 LAPACK。安装的版本是:
我在允许我求逆矩阵的 c 代码中使用 LAPACK。更准确地说,我使用 dgetrf_ 然后使用 dgetri_ 进行反演。 但是当我处理大矩阵并且我不知道矩阵是否可逆时,我浪费了很多时间来反转不可逆
我想使用Fortran和LAPACK对角化一个实对称矩阵。 LAPACK基本上提供了两个例程,一个例程在完整矩阵上运行,另一个例程在打包存储中的矩阵上运行。虽然后者肯定会使用较少的内存,但我想知道关于
我正在使用 LAPACK 库中的 DSYEV 和 DSYEVD 来查找特征值和特征向量(编译语法: gfortran -llapack )。但是,我发现特定矩阵的错误特征值(-0.44,0.35,0.
我想计算对称矩阵的特征值,并希望使用 C++ 中英特尔 MKL 库中的 LAPACKE_dsyev 函数来计算。但我对矩阵需要传递的形式有点困惑。 来自文档https://software.intel
我使用链接:./configure --with-blas-incdir="-L/home/moritz/build/CoinIpopt_test/ThirdParty/openblas/includ
我想安装最新版本的numpy(Python的数值库),版本(v1.6.1)还没有在Ubuntu Oneiric repositories .当我继续手动安装它时,我阅读了 INSTALL numpy
我是一名优秀的程序员,十分优秀!