Hasil salah pemecah eigen jarang
Saya mencoba memecahkan sistem linier jarang Ax = B dengan pustaka Eigen di C ++, namun contoh sepele berikut tampaknya memberikan solusi yang salah:
#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;
}
Saya tidak melihat kesalahan apa pun, algoritme mengembalikan "0" yang berarti "Sukses", namun solusi yang saya dapatkan adalah
x = 0 0.193524 0 0.193524 0.193524 0 0 0 0
yang jelas bukan solusi untuk sistem ini, yang benar adalah
x = 0 0 0 0 0.0967621 0 0 0 0
Jawaban
Berikut dokumentasi untuk SimplicialLDLTpemecah:
Kelas ini menyediakan faktorisasi LDL ^ T Cholesky tanpa akar kuadrat dari matriks renggang yang merupakan titik penjemputan sendiri dan pasti positif .
Ketika matriks menyimpan bilangan real dalam elemen, self-adjoint == simetris. Matriks Anda jelas tidak simetris. Juga, tidak setiap matriks simetris pasti positif-pasti, lihat contoh .
Singkatnya, pemecah yang Anda pilih hanya berlaku untuk kelas matriks yang sangat sempit. Seperti yang telah Anda temukan, SparseLUpemecah bekerja untuk data masukan Anda.
ConjugateGradientsolver tidak akan bekerja baik, tidak memerlukan matriks positif-yang pasti tapi itu tidak memerlukannya untuk menjadi diri-adjoint.