使用LAPACKE在C中通过Cholesky分解求解Ax=b遇问题

核心问题：ldb参数设置错误 #

你的代码中调用LAPACKE_dpotrs时，ldb参数设置为m=3，这是错误的。在**行主序（ROW_MAJOR）**模式下：

• ldb表示右侧矩阵B的每行元素个数（即矩阵的列数）
• 你这里nrhs=1（仅一个右侧向量），因此B是一个3行1列的矩阵，每行只有1个元素，ldb应设为1，而非3。

错误的ldb会导致LAPACKE访问错误的内存位置，从而计算出错误的解。

修正后的代码 #

修改ldb的赋值即可：

// Solve Ax = b using Cholesky decomposed A from above
// Note: As per documentation the solution is written to vector b
int nhrs = 1;
int ldb = 1; // 这里改为1，而非m=3
info = LAPACKE_dpotrs(LAPACK_ROW_MAJOR, uplo, n, nhrs, A, lda, b, ldb);
if (info != 0) {
  printf("Failed to solve Ax = b!\n");
}

额外验证点 #

1. 矩阵正定检查：你的矩阵A是正定的（顺序主子式分别为2、3、1，均大于0），因此LAPACKE_dpotrf的info返回值应为0，分解过程无问题。
2. 参数一致性：确保uplo参数在dpotrf和dpotrs中保持一致（均为'U'），你当前的代码已经满足这一点。

编译运行验证 #

修改后重新编译运行：

gcc test_chol.c -lblas -llapack -llapacke -o test_chol
./test_chol

此时输出应为：

x: [1.000000, 1.000000, 1.000000]

内容的提问来源于stack exchange，提问作者chutsu