add a benchmark routine for all sparse linear solvers in Eigen

2 files changed, 242 insertions, 4 deletions

diff --git a/unsupported/Eigen/src/SparseExtra/MarketIO.h b/unsupported/Eigen/src/SparseExtra/MarketIO.h
index 4e2d5dc4e..3b103914a 100644
--- a/unsupported/Eigen/src/SparseExtra/MarketIO.h
+++ b/unsupported/Eigen/src/SparseExtra/MarketIO.h
@@ -166,10 +166,13 @@ bool loadMarket(SparseMatrixType& mat, const std::string& filename)
     if(!readsizes)
     {
       line >> M >> N >> NNZ;
-      readsizes = true;
-      std::cout << "sizes: " << M << "," << N << "," << NNZ << "\n";
-      mat.resize(M,N);
-      mat.reserve(NNZ);
+      if(M > 0 && N > 0 && NNZ > 0) 
+      {
+        readsizes = true;
+        std::cout << "sizes: " << M << "," << N << "," << NNZ << "\n";
+        mat.resize(M,N);
+        mat.reserve(NNZ);
+      }
     }
     else
     { 
diff --git a/unsupported/Eigen/src/SparseExtra/MatrixMarketIterator.h b/unsupported/Eigen/src/SparseExtra/MatrixMarketIterator.h
new file mode 100644
index 000000000..fe8406ea4
--- /dev/null
+++ b/unsupported/Eigen/src/SparseExtra/MatrixMarketIterator.h
@@ -0,0 +1,235 @@
+ 
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2012
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+#ifndef EIGEN_BROWSE_MATRICES_H
+#define EIGEN_BROWSE_MATRICES_H
+
+#include <dirent.h>
+#include <unsupported/Eigen/SparseExtra>
+using namespace Eigen;
+using std::string; 
+
+enum {
+  SPD = 0x100,
+  NonSymmetric = 0x0
+}; 
+
+/** 
+ * @brief Iterator to browse matrices from a specified folder
+ * 
+ * This is used to load all the matrices from a folder. 
+ * The matrices should be in Matrix Market format
+ * It is assumed that the matrices are named as matname.mtx
+ * and matname_SPD.mtx if the matrix is Symmetric and positive definite (or Hermitian)
+ * The right hand side vectors are loaded as well, if they exist.
+ * They should be named as matname_b.mtx. 
+ * Note that the right hand side for a SPD matrix is named as matname_SPD_b.mtx
+ * 
+ * Sometimes a reference solution is available. In this case, it should be named as matname_x.mtx
+ * 
+ * Sample code
+ * \code
+ * 
+ * \endcode
+ * 
+ * \tparam Scalar The scalar type 
+ */
+template <typename Scalar>
+class MatrixMarketIterator 
+{
+  public:
+    typedef Matrix<Scalar,Dynamic,1> VectorType; 
+    typedef SparseMatrix<Scalar,ColMajor> MatrixType; 
+  
+  public:
+    MatrixMarketIterator(const string folder):m_sym(0),m_isvalid(false),m_matIsLoaded(false),m_hasRhs(false),m_hasrefX(false),m_folder(folder)
+    {
+      m_folder_id = opendir(folder.c_str());
+      if (!m_folder_id){
+        m_isvalid = false;
+        std::cerr << "The provided Matrix folder could not be opened \n\n";
+        abort();
+      }
+      Getnextvalidmatrix();
+    }
+    
+    ~MatrixMarketIterator()
+    {
+      if (m_folder_id) closedir(m_folder_id); 
+    }
+    
+    inline MatrixMarketIterator& operator++()
+    {
+      m_matIsLoaded = false;
+      m_hasrefX = false;
+      m_hasRhs = false;
+      Getnextvalidmatrix();
+      return *this;
+    }
+    inline operator bool() { return m_isvalid;}
+    
+    /** Return the sparse matrix corresponding to the current file */
+    inline MatrixType& matrix() 
+    { 
+      // Read the matrix
+      if (m_matIsLoaded) return m_mat;
+      
+      string matrix_file = m_folder + "/" + m_matname + ".mtx";
+      if ( !loadMarket(m_mat, matrix_file)) 
+      {
+        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>();
+      }
+      return m_mat; 
+    }
+    
+    /** Return the right hand side corresponding to the current matrix. 
+     * If the rhs file is not provided, a random rhs is generated
+     */
+    inline VectorType& rhs() 
+    { 
+       // Get the right hand side
+      if (m_hasRhs) return m_rhs;
+      
+      string rhs_file;
+      rhs_file = m_folder + "/" + m_matname + "_b.mtx"; // The pattern is matname_b.mtx
+      m_hasRhs = Fileexists(rhs_file);
+      if (m_hasRhs)
+      {
+        m_rhs.resize(m_mat.cols());
+        m_hasRhs = loadMarketVector(m_rhs, rhs_file);
+      }
+      if (!m_hasRhs)
+      {
+        // Generate a random right hand side
+        if (!m_matIsLoaded) this->matrix(); 
+        m_refX.resize(m_mat.cols());
+        m_refX.setRandom();
+        m_rhs = m_mat * m_refX;
+        m_hasrefX = true;
+        m_hasRhs = true;
+      }
+      return m_rhs; 
+    }
+    
+    /** Return a reference solution
+     * If it is not provided and if the right hand side is not available
+     * then refX is randomly generated such that A*refX = b 
+     * where A and b are the matrix and the rhs. 
+     * Note that when a rhs is provided, refX is not available 
+     */
+    inline VectorType& refX() 
+    { 
+      // Check if a reference solution is provided
+      if (m_hasrefX) return m_refX;
+      
+      string lhs_file;
+      lhs_file = m_folder + "/" + m_matname + "_x.mtx"; 
+      m_hasrefX = Fileexists(lhs_file);
+      if (m_hasrefX)
+      {
+        m_refX.resize(m_mat.cols());
+        m_hasrefX = loadMarketVector(m_refX, lhs_file);
+      }
+      return m_refX; 
+    }
+    
+    inline string& matname() { return m_matname; }
+    
+    inline int sym() { return m_sym; }
+    
+    inline bool hasRhs() {return m_hasRhs; }
+    inline bool hasrefX() {return m_hasrefX; }
+    
+  protected:
+    
+    inline bool Fileexists(string file)
+    {
+      std::ifstream file_id(file.c_str());
+      if (!file_id.good() ) 
+      {
+        return false;
+      }
+      else 
+      {
+        file_id.close();
+        return true;
+      }
+    }
+    
+    void Getnextvalidmatrix( )
+    {
+      // Here, we return with the next valid matrix in the folder
+      while ( (m_curs_id = readdir(m_folder_id)) != NULL) {
+        m_isvalid = false;
+        string curfile;
+        curfile = m_folder + "/" + m_curs_id->d_name;
+        // Discard if it is a folder
+        if (m_curs_id->d_type == DT_DIR) continue; //FIXME This may not be available on non BSD systems
+//         struct stat st_buf; 
+//         stat (curfile.c_str(), &st_buf);
+//         if (S_ISDIR(st_buf.st_mode)) continue;
+        
+        // Determine from the header if it is a matrix or a right hand side 
+        bool isvector,iscomplex;
+        if(!getMarketHeader(curfile,m_sym,iscomplex,isvector)) continue;
+        if(isvector) continue;
+        
+        // Get the matrix name
+        string filename = m_curs_id->d_name;
+        m_matname = filename.substr(0, filename.length()-4); 
+        
+        // Find if the matrix is SPD 
+        size_t found = m_matname.find("SPD");
+        if( (found!=string::npos) && (m_sym == Symmetric) )
+          m_sym = SPD;
+       
+        m_isvalid = true;
+        break; 
+      }
+    }
+    int m_sym; // Symmetry of the matrix
+    MatrixType m_mat; // Current matrix  
+    VectorType m_rhs;  // Current vector
+    VectorType m_refX; // The reference solution, if exists
+    string m_matname; // Matrix Name
+    bool m_isvalid; 
+    bool m_matIsLoaded; // Determine if the matrix has already been loaded from the file
+    bool m_hasRhs; // The right hand side exists
+    bool m_hasrefX; // A reference solution is provided
+    string m_folder;
+    DIR * m_folder_id;
+    struct dirent *m_curs_id; 
+    
+};
+
+#endif
\ No newline at end of file