eigen3, решение разреженных СЛАУ

olega

Хочу решить с помощью библиотеки Eigen 3 следующую разреженную систему:

/1  0  0\   /x1\     /1\

|1 -2  1| * |x2|  =  |0|

\0  0  1/   \x3/     \1/

Написал вот такой код:

#include <Eigen/Sparse>

#include <Eigen/IterativeLinearSolvers>

#include <Eigen/SuperLUSupport>

#include <Eigen/SparseCholesky>

#include <vector>

#include <iostream>





namespace E=Eigen;



typedef E::SparseMatrix<double> SpMat;

typedef E::Triplet<double> Td;



int main {

    E::Vector rhs = E::Vector(3);

    rhs(0) = 4;

    rhs(1) = 0;

    rhs(2) = 4;



    std::vector<Td> T;

    T.push_back(Td(0,0,1.0;

    T.push_back(Td(1,0,1.0;

    T.push_back(Td(1,1,-2.0;

    T.push_back(Td(1,2,1.0;

    T.push_back(Td(2,2,1.0;



    SpMat M = SpMat(3,3);

    M.setFromTriplets(T.begin T.end;



    // эти три варианта дают странный ответ

    // E::SimplicialCholesky<SpMat> solver; //(1)

    // E::SimplicialLDLT<SpMat> solver;     //(2)

    // E::ConjugateGradient<SpMat> solver;  //(3)



    // а так получается правильно

    E::SuperLU<SpMat> solver;               //(4)



    solver.compute(M);



    E::Vector x;

    x = solver.solve(rhs);



    std::cout<<x<<std::endl;

}

Компилирую следующей командой:
g++ -I /usr/include/eigen3/ -I /usr/include/superlu/ -lsuperlu test3.cc -o test3
В варианте, когда раскомментировано (4 получаю правильный ответ (1,1,1). А во всех остальных случаях ответ равен (0.666667,0.333333,1).
Я явно что-то не понимаю, но вот что?..

Serab

А, для холецкого написано в доках:

The sparse matrix A must be selfadjoint and positive definite.

Для остальных может быть так же.
Это же известный косяк, они решают через мет.опты, там не любую матрицу можно решить так.

olega

О, похоже еще BiCGSTAB правильно работает. Наверное, надо лучше читать документацию...

Serab

хотя ладно, Холецкий не через мет.опты, но сопряженный градиент через мет.опты и там эта причина, а для Холецкого просто математически.
Может в твоем случае можно доказать, что будет хорошая матрица?

olega

Хмм, я подумаю. Пока что матрица, очевидно, плохая. Может удастся как-то переформулировать задачу. В любом случае, спасибо за толковый ответ!