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

Я пытаюсь решить линейную систему уравнений с 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]);

Как задано в этом вопросе, необходимо сначала факторизовать матрицу. Предполагается, что dpotrf факторизует, dpotrs решает систему, используя факторизованную матрицу.

Однако мой результат

2.5   4.0   3.5

явно неверен, я проверил это здесь с WolframAlpha

Где моя ошибка?


c linear-algebra lapack
person bogus    schedule 06.02.2014    source источник

Ответы (1)


arrow_upward
2
arrow_downward

Любопытно, что 2.5 4.0 3.5 является решением правой стороны 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

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

До свидания,

Фрэнсис

person francis    schedule 06.02.2014