Решение линейной системы с dpotrs (факторизация Cholesky)

Question

Я пытаюсь решить линейную систему уравнений с клапкой.

Мой код выглядит следующим образом:

//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 является решением РИТ 1 2 3 ... если матрица

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

dpotrf и dpotrs сделаны для symetric положительно определенных матриц ...

Я хотел бы предложить использовать 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

Выбери свой любимый!

Bye,

Фрэнсис