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 
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:
MatrixMarketIterator(const std::string folder)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 
~MatrixMarketIterator()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 */
matrix()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      */
rhs()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      */
refX()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 
matname()149     inline std::string& matname() { return m_matname; }
150 
sym()151     inline int sym() { return m_sym; }
152 
hasRhs()153     inline bool hasRhs() {return m_hasRhs; }
hasrefX()154     inline bool hasrefX() {return m_hasrefX; }
155 
156   protected:
157 
Fileexists(std::string file)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 
Getnextvalidmatrix()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=false;
188         if(!getMarketHeader(curfile,m_sym,iscomplex,isvector)) continue;
189         if(isvector) continue;
190         if (!iscomplex)
191         {
192           if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value)
193             continue;
194         }
195         if (iscomplex)
196         {
197           if(internal::is_same<Scalar, float>::value || internal::is_same<Scalar, double>::value)
198             continue;
199         }
200 
201 
202         // Get the matrix name
203         std::string filename = m_curs_id->d_name;
204         m_matname = filename.substr(0, filename.length()-4);
205 
206         // Find if the matrix is SPD
207         size_t found = m_matname.find("SPD");
208         if( (found!=std::string::npos) && (m_sym != NonSymmetric) )
209           m_sym = SPD;
210 
211         m_isvalid = true;
212         break;
213       }
214     }
215     int m_sym; // Symmetry of the matrix
216     MatrixType m_mat; // Current matrix
217     VectorType m_rhs;  // Current vector
218     VectorType m_refX; // The reference solution, if exists
219     std::string m_matname; // Matrix Name
220     bool m_isvalid;
221     bool m_matIsLoaded; // Determine if the matrix has already been loaded from the file
222     bool m_hasRhs; // The right hand side exists
223     bool m_hasrefX; // A reference solution is provided
224     std::string m_folder;
225     DIR * m_folder_id;
226     struct dirent *m_curs_id;
227 
228 };
229 
230 } // end namespace Eigen
231 
232 #endif
233