用dpotrs解決線性系統（Cholesky分解）

Question

我想用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

哪裏是我的錯？

Answer 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

選擇您最喜歡的！

再見，

弗朗西斯