使用 dpotrs 求解线性系统（Cholesky 分解）

【问题标题】：Solving a linear system with dpotrs (Cholesky factorization)使用 dpotrs 求解线性系统（Cholesky 分解）
【发布时间】：2014-02-06 20:06:06
【问题描述】：

我正在尝试使用 clapack 求解线性方程组。

我的代码如下：

//ATTENTION: matrix in column-major
double A[3*3]={  2.0, -1.0,  0.0,
                0.0,  2.0, -1.0,
                0.0,  0.0,  2.0},
b[3]={1.0,2.0,3.0};

integer n=3,info,nrhs=1; char uplo='L';

dpotrf_("L", &n, A, &n, &info);
dpotrs_("L", &n, &nrhs, A, &n, b, &n, &info);

printf("Solution: %10.4f %10.4f %10.4f",b[0], b[1], b[2]);

正如this question 中所问的，有必要首先对矩阵进行因式分解。 dpotrf 应该是因式分解，dpotrs 使用分解后的矩阵求解系统。

但是，我的结果

2.5   4.0   3.5

明显错了，我这里查了with WolframAlpha

我的错误在哪里？

【问题讨论】：

标签： c linear-algebra lapack

【解决方案1】：

奇怪的是，2.5 4.0 3.5 是 rhs 的解 1 2 3...如果矩阵是

 2   -1  0
 -1  2   -1
 0   -1  2

dpotrf 和 dpotrs 用于对称正定矩阵...

我建议改用dgetrf 和dgetrs。在你的情况下：

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <lapacke.h>

int main()
{
    //ATTENTION: matrix in column-major
    double A[3*3]={  2.0, -1.0,  0.0,
            0,  2.0, -1.0,
            0.0,  0,  2.0},
            b[3]={1.0,2,3.0};

    int n=3,info,nrhs=1; char uplo='L';
    int ipvs[3];
    dgetrf_(&n, &n, A, &n,ipvs, &info);
    printf("info %d\n",info);
    dgetrs_("T", &n, &nrhs, A, &n,ipvs, b, &n, &info);
    printf("info %d\n",info);
    printf("Solution: %10.4f %10.4f %10.4f\n",b[0], b[1], b[2]);

    return 0;
}

我得到的解决方案是1.3750 1.75 1.5。如果不是列主序，切换到"N"而不是"T"。那么，解决方案就变成了0.5 1.25 2.125

选择你最喜欢的！

再见，

弗朗西斯

【讨论】：