Home | History | Annotate | Download | only in SparseExtra
      1 
      2 // This file is part of Eigen, a lightweight C++ template library
      3 // for linear algebra.
      4 //
      5 // Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam (at) inria.fr>
      6 //
      7 // This Source Code Form is subject to the terms of the Mozilla
      8 // Public License v. 2.0. If a copy of the MPL was not distributed
      9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
     10 
     11 #ifndef EIGEN_BROWSE_MATRICES_H
     12 #define EIGEN_BROWSE_MATRICES_H
     13 
     14 namespace Eigen {
     15 
     16 enum {
     17   SPD = 0x100,
     18   NonSymmetric = 0x0
     19 };
     20 
     21 /**
     22  * @brief Iterator to browse matrices from a specified folder
     23  *
     24  * This is used to load all the matrices from a folder.
     25  * The matrices should be in Matrix Market format
     26  * It is assumed that the matrices are named as matname.mtx
     27  * and matname_SPD.mtx if the matrix is Symmetric and positive definite (or Hermitian)
     28  * The right hand side vectors are loaded as well, if they exist.
     29  * They should be named as matname_b.mtx.
     30  * Note that the right hand side for a SPD matrix is named as matname_SPD_b.mtx
     31  *
     32  * Sometimes a reference solution is available. In this case, it should be named as matname_x.mtx
     33  *
     34  * Sample code
     35  * \code
     36  *
     37  * \endcode
     38  *
     39  * \tparam Scalar The scalar type
     40  */
     41 template <typename Scalar>
     42 class MatrixMarketIterator
     43 {
     44   public:
     45     typedef Matrix<Scalar,Dynamic,1> VectorType;
     46     typedef SparseMatrix<Scalar,ColMajor> MatrixType;
     47 
     48   public:
     49     MatrixMarketIterator(const std::string folder):m_sym(0),m_isvalid(false),m_matIsLoaded(false),m_hasRhs(false),m_hasrefX(false),m_folder(folder)
     50     {
     51       m_folder_id = opendir(folder.c_str());
     52       if (!m_folder_id){
     53         m_isvalid = false;
     54         std::cerr << "The provided Matrix folder could not be opened \n\n";
     55         abort();
     56       }
     57       Getnextvalidmatrix();
     58     }
     59 
     60     ~MatrixMarketIterator()
     61     {
     62       if (m_folder_id) closedir(m_folder_id);
     63     }
     64 
     65     inline MatrixMarketIterator& operator++()
     66     {
     67       m_matIsLoaded = false;
     68       m_hasrefX = false;
     69       m_hasRhs = false;
     70       Getnextvalidmatrix();
     71       return *this;
     72     }
     73     inline operator bool() const { return m_isvalid;}
     74 
     75     /** Return the sparse matrix corresponding to the current file */
     76     inline MatrixType& matrix()
     77     {
     78       // Read the matrix
     79       if (m_matIsLoaded) return m_mat;
     80 
     81       std::string matrix_file = m_folder + "/" + m_matname + ".mtx";
     82       if ( !loadMarket(m_mat, matrix_file))
     83       {
     84         m_matIsLoaded = false;
     85         return m_mat;
     86       }
     87       m_matIsLoaded = true;
     88 
     89       if (m_sym != NonSymmetric)
     90       { // Store the upper part of the matrix. It is needed by the solvers dealing with nonsymmetric matrices ??
     91         MatrixType B;
     92         B = m_mat;
     93         m_mat = B.template selfadjointView<Lower>();
     94       }
     95       return m_mat;
     96     }
     97 
     98     /** Return the right hand side corresponding to the current matrix.
     99      * If the rhs file is not provided, a random rhs is generated
    100      */
    101     inline VectorType& rhs()
    102     {
    103        // Get the right hand side
    104       if (m_hasRhs) return m_rhs;
    105 
    106       std::string rhs_file;
    107       rhs_file = m_folder + "/" + m_matname + "_b.mtx"; // The pattern is matname_b.mtx
    108       m_hasRhs = Fileexists(rhs_file);
    109       if (m_hasRhs)
    110       {
    111         m_rhs.resize(m_mat.cols());
    112         m_hasRhs = loadMarketVector(m_rhs, rhs_file);
    113       }
    114       if (!m_hasRhs)
    115       {
    116         // Generate a random right hand side
    117         if (!m_matIsLoaded) this->matrix();
    118         m_refX.resize(m_mat.cols());
    119         m_refX.setRandom();
    120         m_rhs = m_mat * m_refX;
    121         m_hasrefX = true;
    122         m_hasRhs = true;
    123       }
    124       return m_rhs;
    125     }
    126 
    127     /** Return a reference solution
    128      * If it is not provided and if the right hand side is not available
    129      * then refX is randomly generated such that A*refX = b
    130      * where A and b are the matrix and the rhs.
    131      * Note that when a rhs is provided, refX is not available
    132      */
    133     inline VectorType& refX()
    134     {
    135       // Check if a reference solution is provided
    136       if (m_hasrefX) return m_refX;
    137 
    138       std::string lhs_file;
    139       lhs_file = m_folder + "/" + m_matname + "_x.mtx";
    140       m_hasrefX = Fileexists(lhs_file);
    141       if (m_hasrefX)
    142       {
    143         m_refX.resize(m_mat.cols());
    144         m_hasrefX = loadMarketVector(m_refX, lhs_file);
    145       }
    146       return m_refX;
    147     }
    148 
    149     inline std::string& matname() { return m_matname; }
    150 
    151     inline int sym() { return m_sym; }
    152 
    153     inline bool hasRhs() {return m_hasRhs; }
    154     inline bool hasrefX() {return m_hasrefX; }
    155 
    156   protected:
    157 
    158     inline bool Fileexists(std::string file)
    159     {
    160       std::ifstream file_id(file.c_str());
    161       if (!file_id.good() )
    162       {
    163         return false;
    164       }
    165       else
    166       {
    167         file_id.close();
    168         return true;
    169       }
    170     }
    171 
    172     void Getnextvalidmatrix( )
    173     {
    174       m_isvalid = false;
    175       // Here, we return with the next valid matrix in the folder
    176       while ( (m_curs_id = readdir(m_folder_id)) != NULL) {
    177         m_isvalid = false;
    178         std::string curfile;
    179         curfile = m_folder + "/" + m_curs_id->d_name;
    180         // Discard if it is a folder
    181         if (m_curs_id->d_type == DT_DIR) continue; //FIXME This may not be available on non BSD systems
    182 //         struct stat st_buf;
    183 //         stat (curfile.c_str(), &st_buf);
    184 //         if (S_ISDIR(st_buf.st_mode)) continue;
    185 
    186         // Determine from the header if it is a matrix or a right hand side
    187         bool isvector,iscomplex;
    188         if(!getMarketHeader(curfile,m_sym,iscomplex,isvector)) continue;
    189         if(isvector) continue;
    190 
    191         // Get the matrix name
    192         std::string filename = m_curs_id->d_name;
    193         m_matname = filename.substr(0, filename.length()-4);
    194 
    195         // Find if the matrix is SPD
    196         size_t found = m_matname.find("SPD");
    197         if( (found!=std::string::npos) && (m_sym != NonSymmetric) )
    198           m_sym = SPD;
    199 
    200         m_isvalid = true;
    201         break;
    202       }
    203     }
    204     int m_sym; // Symmetry of the matrix
    205     MatrixType m_mat; // Current matrix
    206     VectorType m_rhs;  // Current vector
    207     VectorType m_refX; // The reference solution, if exists
    208     std::string m_matname; // Matrix Name
    209     bool m_isvalid;
    210     bool m_matIsLoaded; // Determine if the matrix has already been loaded from the file
    211     bool m_hasRhs; // The right hand side exists
    212     bool m_hasrefX; // A reference solution is provided
    213     std::string m_folder;
    214     DIR * m_folder_id;
    215     struct dirent *m_curs_id;
    216 
    217 };
    218 
    219 } // end namespace Eigen
    220 
    221 #endif
    222