Решение большой разреженной линейной системы с помощью SimplicialCholesky Eigen - PullRequest
Новая
       55

Решение большой разреженной линейной системы с помощью SimplicialCholesky Eigen

0 голосов
/ 09 апреля 2019

Я на самом деле пытаюсь решить большие разреженные линейные системы, используя C ++ lib Eigen. Разреженные матрицы взяты из этой страницы . Каждая система имеет такую ​​структуру: Ax = b, где A - это разреженная матрица (n x n), b вычисляется как A*xe с xe вектором измерения n, содержащим только нули. После вычисления x мне нужно вычислить относительную ошибку между xe и x. Я написал некоторый код, но не могу понять, почему относительная ошибка настолько высока (1.49853e + 08) в конце вычисления. 

#include <iostream>
#include <Eigen/Dense>
#include <unsupported/Eigen/SparseExtra>
#include<Eigen/SparseCholesky>
#include <sys/time.h>
#include <sys/resource.h>


using namespace std;
using namespace Eigen;


int main()
{

    SparseMatrix<double> mat;
    loadMarket(mat, "/Users/anto/Downloads/ex15/ex15.mtx");

	VectorXd xe = VectorXd::Constant(mat.rows(), 1);
	VectorXd b = mat*xe;

    
    SimplicialCholesky<Eigen::SparseMatrix<double> > chol(mat);
    VectorXd x = chol.solve(b); 

    double relative_error = (x-xe).norm()/(xe).norm(); 
    cout << relative_error << endl;
    
}

Матрица ex15 может быть загружена с этой страницы . Это симметричная положительно определенная матрица. Может ли кто-нибудь помочь мне решить проблему? Заранее благодарю за помощь.

1 Ответ

0 голосов
/ 09 апреля 2019

Согласно эта страница , ex15 не является полным рангом. Вы должны проверить, что каждый шаг прошел хорошо: 

SimplicialLDLT<Eigen::SparseMatrix<double> > chol(mat);
if(chol.info()!=Eigen::Success)
  return;
VectorXd x = chol.solve(b); 
if(chol.info()!=Eigen::Success)
  return;

и затем проверьте, что у вас есть одно решение (если оно не является полным рангом и существует хотя бы одно решение, то существует целое подпространство решений): 

cout << (mat*x-b).norm()/b.norm() << "\n";
...