MatrixMarketIterator.h
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 // Eigen is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU Lesser General Public
9 // License as published by the Free Software Foundation; either
10 // version 3 of the License, or (at your option) any later version.
11 //
12 // Alternatively, you can redistribute it and/or
13 // modify it under the terms of the GNU General Public License as
14 // published by the Free Software Foundation; either version 2 of
15 // the License, or (at your option) any later version.
16 //
17 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
18 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
20 // GNU General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License and a copy of the GNU General Public License along with
24 // Eigen. If not, see <http://www.gnu.org/licenses/>.
25 
26 #ifndef EIGEN_BROWSE_MATRICES_H
27 #define EIGEN_BROWSE_MATRICES_H
28 
29 namespace Eigen {
30 
31 enum {
32  SPD = 0x100,
33  NonSymmetric = 0x0
34 };
35 
56 template <typename Scalar>
58 {
59  public:
60  typedef Matrix<Scalar,Dynamic,1> VectorType;
61  typedef SparseMatrix<Scalar,ColMajor> MatrixType;
62 
63  public:
64  MatrixMarketIterator(const std::string folder):m_sym(0),m_isvalid(false),m_matIsLoaded(false),m_hasRhs(false),m_hasrefX(false),m_folder(folder)
65  {
66  m_folder_id = opendir(folder.c_str());
67  if (!m_folder_id){
68  m_isvalid = false;
69  std::cerr << "The provided Matrix folder could not be opened \n\n";
70  abort();
71  }
72  Getnextvalidmatrix();
73  }
74 
76  {
77  if (m_folder_id) closedir(m_folder_id);
78  }
79 
80  inline MatrixMarketIterator& operator++()
81  {
82  m_matIsLoaded = false;
83  m_hasrefX = false;
84  m_hasRhs = false;
85  Getnextvalidmatrix();
86  return *this;
87  }
88  inline operator bool() const { return m_isvalid;}
89 
91  inline MatrixType& matrix()
92  {
93  // Read the matrix
94  if (m_matIsLoaded) return m_mat;
95 
96  std::string matrix_file = m_folder + "/" + m_matname + ".mtx";
97  if ( !loadMarket(m_mat, matrix_file))
98  {
99  m_matIsLoaded = false;
100  return m_mat;
101  }
102  m_matIsLoaded = true;
103 
104  if (m_sym != NonSymmetric)
105  { // Store the upper part of the matrix. It is needed by the solvers dealing with nonsymmetric matrices ??
106  MatrixType B;
107  B = m_mat;
108  m_mat = B.template selfadjointView<Lower>();
109  }
110  return m_mat;
111  }
112 
116  inline VectorType& rhs()
117  {
118  // Get the right hand side
119  if (m_hasRhs) return m_rhs;
120 
121  std::string rhs_file;
122  rhs_file = m_folder + "/" + m_matname + "_b.mtx"; // The pattern is matname_b.mtx
123  m_hasRhs = Fileexists(rhs_file);
124  if (m_hasRhs)
125  {
126  m_rhs.resize(m_mat.cols());
127  m_hasRhs = loadMarketVector(m_rhs, rhs_file);
128  }
129  if (!m_hasRhs)
130  {
131  // Generate a random right hand side
132  if (!m_matIsLoaded) this->matrix();
133  m_refX.resize(m_mat.cols());
134  m_refX.setRandom();
135  m_rhs = m_mat * m_refX;
136  m_hasrefX = true;
137  m_hasRhs = true;
138  }
139  return m_rhs;
140  }
141 
148  inline VectorType& refX()
149  {
150  // Check if a reference solution is provided
151  if (m_hasrefX) return m_refX;
152 
153  std::string lhs_file;
154  lhs_file = m_folder + "/" + m_matname + "_x.mtx";
155  m_hasrefX = Fileexists(lhs_file);
156  if (m_hasrefX)
157  {
158  m_refX.resize(m_mat.cols());
159  m_hasrefX = loadMarketVector(m_refX, lhs_file);
160  }
161  return m_refX;
162  }
163 
164  inline std::string& matname() { return m_matname; }
165 
166  inline int sym() { return m_sym; }
167 
168  inline bool hasRhs() {return m_hasRhs; }
169  inline bool hasrefX() {return m_hasrefX; }
170 
171  protected:
172 
173  inline bool Fileexists(std::string file)
174  {
175  std::ifstream file_id(file.c_str());
176  if (!file_id.good() )
177  {
178  return false;
179  }
180  else
181  {
182  file_id.close();
183  return true;
184  }
185  }
186 
187  void Getnextvalidmatrix( )
188  {
189  m_isvalid = false;
190  // Here, we return with the next valid matrix in the folder
191  while ( (m_curs_id = readdir(m_folder_id)) != NULL) {
192  m_isvalid = false;
193  std::string curfile;
194  curfile = m_folder + "/" + m_curs_id->d_name;
195  // Discard if it is a folder
196  if (m_curs_id->d_type == DT_DIR) continue; //FIXME This may not be available on non BSD systems
197 // struct stat st_buf;
198 // stat (curfile.c_str(), &st_buf);
199 // if (S_ISDIR(st_buf.st_mode)) continue;
200 
201  // Determine from the header if it is a matrix or a right hand side
202  bool isvector,iscomplex;
203  if(!getMarketHeader(curfile,m_sym,iscomplex,isvector)) continue;
204  if(isvector) continue;
205 
206  // Get the matrix name
207  std::string filename = m_curs_id->d_name;
208  m_matname = filename.substr(0, filename.length()-4);
209 
210  // Find if the matrix is SPD
211  size_t found = m_matname.find("SPD");
212  if( (found!=std::string::npos) && (m_sym != NonSymmetric) )
213  m_sym = SPD;
214 
215  m_isvalid = true;
216  break;
217  }
218  }
219  int m_sym; // Symmetry of the matrix
220  MatrixType m_mat; // Current matrix
221  VectorType m_rhs; // Current vector
222  VectorType m_refX; // The reference solution, if exists
223  std::string m_matname; // Matrix Name
224  bool m_isvalid;
225  bool m_matIsLoaded; // Determine if the matrix has already been loaded from the file
226  bool m_hasRhs; // The right hand side exists
227  bool m_hasrefX; // A reference solution is provided
228  std::string m_folder;
229  DIR * m_folder_id;
230  struct dirent *m_curs_id;
231 
232 };
233 
234 } // end namespace Eigen
235 
236 #endif