Improve loading of symmetric sparse matrices in MatrixMarketIterator

2 files changed, 23 insertions, 6 deletions

diff --git a/unsupported/Eigen/src/SparseExtra/MarketIO.h b/unsupported/Eigen/src/SparseExtra/MarketIO.h
index 100e617b2..cdc14f86e 100644
--- a/unsupported/Eigen/src/SparseExtra/MarketIO.h
+++ b/unsupported/Eigen/src/SparseExtra/MarketIO.h
@@ -163,7 +163,7 @@ bool loadMarket(SparseMatrixType& mat, const std::string& filename)
       if(M > 0 && N > 0 && NNZ > 0) 
       {
         readsizes = true;
-        std::cout << "sizes: " << M << "," << N << "," << NNZ << "\n";
+        //std::cout << "sizes: " << M << "," << N << "," << NNZ << "\n";
         mat.resize(M,N);
         mat.reserve(NNZ);
       }
diff --git a/unsupported/Eigen/src/SparseExtra/MatrixMarketIterator.h b/unsupported/Eigen/src/SparseExtra/MatrixMarketIterator.h
index bf13cf21f..4a40b26be 100644
--- a/unsupported/Eigen/src/SparseExtra/MatrixMarketIterator.h
+++ b/unsupported/Eigen/src/SparseExtra/MatrixMarketIterator.h
@@ -41,6 +41,7 @@ enum {
 template <typename Scalar>
 class MatrixMarketIterator 
 {
+    typedef typename NumTraits<Scalar>::Real RealScalar;
   public:
     typedef Matrix<Scalar,Dynamic,1> VectorType; 
     typedef SparseMatrix<Scalar,ColMajor> MatrixType; 
@@ -81,16 +82,30 @@ class MatrixMarketIterator
       std::string matrix_file = m_folder + "/" + m_matname + ".mtx";
       if ( !loadMarket(m_mat, matrix_file)) 
       {
+        std::cerr << "Warning loadMarket failed when loading \"" << matrix_file << "\"\n";
         m_matIsLoaded = false;
         return m_mat;
       }
       m_matIsLoaded = true; 
-      
+
       if (m_sym != NonSymmetric) 
-      { // Store the upper part of the matrix. It is needed by the solvers dealing with nonsymmetric matrices ??
-        MatrixType B; 
-        B = m_mat;
-        m_mat = B.template selfadjointView<Lower>();
+      {
+        // Check whether we need to restore a full matrix:
+        RealScalar diag_norm  = m_mat.diagonal().norm();
+        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; 
     }
@@ -143,6 +158,8 @@ class MatrixMarketIterator
         m_refX.resize(m_mat.cols());
         m_hasrefX = loadMarketVector(m_refX, lhs_file);
       }
+      else
+        m_refX.resize(0);
       return m_refX; 
     }