Eigen sparse solver błędne wyniki
Próbuję rozwiązać rzadki system liniowy Ax = B za pomocą biblioteki Eigen w C ++, jednak następujący trywialny przykład wydaje się dawać nieprawidłowe rozwiązanie:
#include <Eigen/SparseCholesky>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <iostream>
#include <vector>
using namespace std;
using namespace Eigen;
int main(){
SimplicialLDLT<SparseMatrix<double>> solver;
SparseMatrix<double> A(9,9);
typedef Triplet<double> T;
vector<T> triplet;
VectorXd B(9);
for(int i=0; i<4; i++){
triplet.push_back(T(i,i,1));
triplet.push_back(T(i+5,i+5,1));
}
triplet.push_back(T(4,1,-1));
triplet.push_back(T(4,3,-1));
triplet.push_back(T(4,5,-1));
triplet.push_back(T(4,7,-1));
triplet.push_back(T(4,4,4));
A.setFromTriplets(triplet.begin(),triplet.end());
B << 0,0,0,0,0.387049,0,0,0,0;
solver.compute(A);
VectorXd x = solver.solve(B);
cout << "A\n" << A << "\n";
cout << "B\n" << B << "\n";
cout << "x\n" << x << "\n";
return 0;
}
Nie widzę żadnych błędów, algorytm zwraca „0” co oznacza „Sukces”, jednak rozwiązanie jakie otrzymuję to
x = 0 0.193524 0 0.193524 0.193524 0 0 0 0
co oczywiście nie jest rozwiązaniem tego systemu, jest właściwe
x = 0 0 0 0 0.0967621 0 0 0 0
Odpowiedzi
Oto dokumentacja dla SimplicialLDLTsolvera:
Ta klasa zapewnia czynniki LDL ^ T Cholesky'ego bez pierwiastka kwadratowego z rzadkich macierzy, które są samosprzężone i dodatnio określone .
Gdy macierz przechowuje liczby rzeczywiste w elementach, samosprzężenie == symetryczne. Twoja matryca najwyraźniej nie jest symetryczna. Ponadto nie każda macierz symetryczna ma charakter dodatni, patrz przykłady .
Krótko mówiąc, wybrany solver ma zastosowanie tylko do bardzo wąskiej klasy macierzy. Jak już zauważyłeś, SparseLUsolver działa dla twoich danych wejściowych.
ConjugateGradientSolver nie będzie działać albo, że nie wymaga matrycy do dodatnio określony, ale to nie wymaga to być samosprzężony.