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

Question

我尝试用 C 语言编写求解器函数(在 Code::Blocks 中使用 Ubuntu)。我想从 Python 脚本调用该求解器，因此我将其构建为共享库，并在 Python 中使用“ctypes”进行调用。 首先尝试了那里的例子 http://www.netlib.org/blas/blast-forum/chapter3.pdf (第 12 页)。这是一个简单的矩阵(稀疏)- vector (密集)乘法。一切都运转良好。 这是乘法代码:

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并返回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 。但这个函数只是返回 -1 (出了问题)并且不计算任何东西。

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

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

Answer 1

您可以在第 121 页的链接中阅读 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，这可能是一个问题。

<小时/>

另一方面，我建议您使用更强大/成熟的稀疏求解器库，例如 MKL sparse solvers ，或者更容易使用 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()。