共享库中使用的稀疏 BLAS 求解器不起作用 (returns '-1')

Question

我曾尝试用 C 编写一个求解器函数（在 Code::Blocks 中使用 Ubuntu）。我想从 Python 脚本调用该求解器，因此我将其构建为共享库并在 Python 中使用 'ctypes' 调用它。首先尝试了那里的示例 http://www.netlib.org/blas/blast-forum/chapter3.pdf（第 12 页）。这是一个简单的矩阵（稀疏）-向量（密集）乘法。一切正常。这里的乘法代码：

const int n =4;
const int nz=6;
double val[] = {1.1,2.2,3.3,4.1,2.4,4.4};
int indx[]={0,1,2,3,3,3};
int jndx[]={0,1,2,0,1,3};

double x[]={1.0,1.0,1.0,1.0};
double y[]={0.0,0.0,0.0,0.0};

blas_sparse_matrix A;
double alph=1.0;


A=BLAS_duscr_begin(n,n);
for(i=0;i<nz;i++){
    BLAS_duscr_insert_entry(A,val[i],indx[i],jndx[i]);
}

BLAS_duscr_end(A);

int err=BLAS_ussp(A,blas_lower_triangular);
printf("%i\n",err);

if(BLAS_usgp(A, blas_lower_triangular) || BLAS_usgp(A, blas_upper_triangular)){
    BLAS_dussv(blas_no_trans,alph,A,x,n);
    printf("SOLVE... ");
}
else{
    printf("FAIL...");
}

BLAS_usds(A);

for(i=0;i<n;i++){
    printf("%f\n",x[i]);
}

当一切都成功时，它会执行 y<-A*x 和 returns 0。

现在我刚刚更改了这一行：

BLAS_dusmv(blas_no_trans, alph,A,x,1,y,1);

收件人：

BLAS_dussv(blas_no_trans,alph,A,x,1);

这是一个执行 x <- A^(-1)*x 的三角求解器。但是这个函数只是returns -1（出了点问题）并且没有计算任何东西。

有人知道如何正确使用稀疏BLAS的求解器吗？

求解器是只求解三角矩阵还是将其解释为方阵？（取矩阵的下半部分，假设为对称矩阵）

Answer 1

您可以在第 121 页的 link 中阅读 dussv() 的文档。

http://www.netlib.org/blas/blast-forum/chapter3.pdf

dussv()要求你的A是三角矩阵，也就是说

either USGP(A, blas_lower_triangular) or USGP(A, blas_upper_triangular) must be true.

根据你的代码，你指定了两次A(3,3)分别是4.4和20.0，这可能是个问题。

另一方面，我建议您使用更多 powerful/mature 个稀疏求解器库，例如 MKL sparse solvers, or eaiser to use C++ Eigen。他们都可以解决任意形状的A。

最后，python 在 Scipy 中已经有了稀疏求解器，你为什么要实现自己的？

http://docs.scipy.org/doc/scipy/reference/sparse.linalg.html#solving-linear-problems

编辑

修改后，你的问题变为：

And now I think the problem lies in the function BLAS_ussp(A,blas_lower_triangular); which isn't able to set the 'blas_lower_triangular' flag (returns -1).

请阅读第 130 页 ussp() 的文档。您需要在任何 INSERT 例程之前调用 ussp()。

共享库中使用的稀疏 BLAS 求解器不起作用 (returns '-1')

sparse BLAS solver used in a shared library doesn't work (returns '-1')

c

python

solver

blas

sparse-matrix