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
