GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-07-19 19:34:29 +08:00
Code Issues Packages Projects Releases Wiki Activity

Improve loading of symmetric sparse matrices in MatrixMarketIterator

Browse Source
This commit is contained in:
Gael Guennebaud 2015-06-05 10:16:29 +02:00
parent acc761cf0c
commit 3e7bc8d686
2 changed files with 23 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
unsupported/Eigen/src/SparseExtra/MarketIO.h
View File

@ -163,7 +163,7 @@ bool loadMarket(SparseMatrixType& mat, const std::string& filename)
if(M > 0 && N > 0 && NNZ > 0) if(M > 0 && N > 0 && NNZ > 0)
{ {
readsizes = true; readsizes = true;
std::cout << "sizes: " << M << "," << N << "," << NNZ << "\n"; //std::cout << "sizes: " << M << "," << N << "," << NNZ << "\n";
mat.resize(M,N); mat.resize(M,N);
mat.reserve(NNZ); mat.reserve(NNZ);
} }

25
unsupported/Eigen/src/SparseExtra/MatrixMarketIterator.h
View File

@ -41,6 +41,7 @@ enum {
template <typename Scalar> template <typename Scalar>
class MatrixMarketIterator class MatrixMarketIterator
{ {
typedef typename NumTraits<Scalar>::Real RealScalar;
public: public:
typedef Matrix<Scalar,Dynamic,1> VectorType; typedef Matrix<Scalar,Dynamic,1> VectorType;
typedef SparseMatrix<Scalar,ColMajor> MatrixType; typedef SparseMatrix<Scalar,ColMajor> MatrixType;
@ -81,16 +82,30 @@ class MatrixMarketIterator
std::string matrix_file = m_folder + "/" + m_matname + ".mtx"; std::string matrix_file = m_folder + "/" + m_matname + ".mtx";
if ( !loadMarket(m_mat, matrix_file)) if ( !loadMarket(m_mat, matrix_file))
{ {
std::cerr << "Warning loadMarket failed when loading \"" << matrix_file << "\"\n";
m_matIsLoaded = false; m_matIsLoaded = false;
return m_mat; return m_mat;
} }
m_matIsLoaded = true; m_matIsLoaded = true;
if (m_sym != NonSymmetric) if (m_sym != NonSymmetric)
{ // Store the upper part of the matrix. It is needed by the solvers dealing with nonsymmetric matrices ?? {
MatrixType B; // Check whether we need to restore a full matrix:
B = m_mat; RealScalar diag_norm = m_mat.diagonal().norm();
m_mat = B.template selfadjointView<Lower>(); RealScalar lower_norm = m_mat.template triangularView<Lower>().norm();
RealScalar upper_norm = m_mat.template triangularView<Upper>().norm();
if(lower_norm>diag_norm && upper_norm==diag_norm)
{
// only the lower part is stored
MatrixType tmp(m_mat);
m_mat = tmp.template selfadjointView<Lower>();
}
else if(upper_norm>diag_norm && lower_norm==diag_norm)
{
// only the upper part is stored
MatrixType tmp(m_mat);
m_mat = tmp.template selfadjointView<Upper>();
}
} }
return m_mat; return m_mat;
} }
@ -143,6 +158,8 @@ class MatrixMarketIterator
m_refX.resize(m_mat.cols()); m_refX.resize(m_mat.cols());
m_hasrefX = loadMarketVector(m_refX, lhs_file); m_hasrefX = loadMarketVector(m_refX, lhs_file);
} }
else
m_refX.resize(0);
return m_refX; return m_refX;
} }