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@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
14namespace Eigen {
15
16enum {
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 */
41template <typename Scalar>
42class 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