eigen/unsupported/Eigen/src/SparseExtra/MatrixMarketIterator.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_BROWSE_MATRICES_H
#define EIGEN_BROWSE_MATRICES_H

namespace Eigen {

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 std::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() const { return m_isvalid;}
    
    /** Return the sparse matrix corresponding to the current file */
    inline MatrixType& matrix() 
    { 
      // Read the matrix
      if (m_matIsLoaded) return m_mat;
      
      std::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;
      
      std::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;
      
      std::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 std::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(std::string file)
    {
      std::ifstream file_id(file.c_str());
      if (!file_id.good() ) 
      {
        return false;
      }
      else 
      {
        file_id.close();
        return true;
      }
    }
    
    void Getnextvalidmatrix( )
    {
      m_isvalid = false;
      // Here, we return with the next valid matrix in the folder
      while ( (m_curs_id = readdir(m_folder_id)) != NULL) {
        m_isvalid = false;
        std::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=false;
        if(!getMarketHeader(curfile,m_sym,iscomplex,isvector)) continue;
        if(isvector) continue;
        if (!iscomplex)
        {
          if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value)
            continue; 
        }
        if (iscomplex)
        {
          if(internal::is_same<Scalar, float>::value || internal::is_same<Scalar, double>::value)
            continue; 
        }
        
        
        // Get the matrix name
        std::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!=std::string::npos) && (m_sym != NonSymmetric) )
          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
    std::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
    std::string m_folder;
    DIR * m_folder_id;
    struct dirent *m_curs_id; 
    
};

} // end namespace Eigen

#endif