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 typedef typename NumTraits<Scalar>::Real RealScalar; 45 public: 46 typedef Matrix<Scalar,Dynamic,1> VectorType; 47 typedef SparseMatrix<Scalar,ColMajor> MatrixType; 48 49 public: 50 MatrixMarketIterator(const std::string &folder) 51 : m_sym(0), m_isvalid(false), m_matIsLoaded(false), m_hasRhs(false), m_hasrefX(false), m_folder(folder) 52 { 53 m_folder_id = opendir(folder.c_str()); 54 if(m_folder_id) 55 Getnextvalidmatrix(); 56 } 57 58 ~MatrixMarketIterator() 59 { 60 if (m_folder_id) closedir(m_folder_id); 61 } 62 63 inline MatrixMarketIterator& operator++() 64 { 65 m_matIsLoaded = false; 66 m_hasrefX = false; 67 m_hasRhs = false; 68 Getnextvalidmatrix(); 69 return *this; 70 } 71 inline operator bool() const { return m_isvalid;} 72 73 /** Return the sparse matrix corresponding to the current file */ 74 inline MatrixType& matrix() 75 { 76 // Read the matrix 77 if (m_matIsLoaded) return m_mat; 78 79 std::string matrix_file = m_folder + "/" + m_matname + ".mtx"; 80 if ( !loadMarket(m_mat, matrix_file)) 81 { 82 std::cerr << "Warning loadMarket failed when loading \"" << matrix_file << "\"" << std::endl; 83 m_matIsLoaded = false; 84 return m_mat; 85 } 86 m_matIsLoaded = true; 87 88 if (m_sym != NonSymmetric) 89 { 90 // Check whether we need to restore a full matrix: 91 RealScalar diag_norm = m_mat.diagonal().norm(); 92 RealScalar lower_norm = m_mat.template triangularView<Lower>().norm(); 93 RealScalar upper_norm = m_mat.template triangularView<Upper>().norm(); 94 if(lower_norm>diag_norm && upper_norm==diag_norm) 95 { 96 // only the lower part is stored 97 MatrixType tmp(m_mat); 98 m_mat = tmp.template selfadjointView<Lower>(); 99 } 100 else if(upper_norm>diag_norm && lower_norm==diag_norm) 101 { 102 // only the upper part is stored 103 MatrixType tmp(m_mat); 104 m_mat = tmp.template selfadjointView<Upper>(); 105 } 106 } 107 return m_mat; 108 } 109 110 /** Return the right hand side corresponding to the current matrix. 111 * If the rhs file is not provided, a random rhs is generated 112 */ 113 inline VectorType& rhs() 114 { 115 // Get the right hand side 116 if (m_hasRhs) return m_rhs; 117 118 std::string rhs_file; 119 rhs_file = m_folder + "/" + m_matname + "_b.mtx"; // The pattern is matname_b.mtx 120 m_hasRhs = Fileexists(rhs_file); 121 if (m_hasRhs) 122 { 123 m_rhs.resize(m_mat.cols()); 124 m_hasRhs = loadMarketVector(m_rhs, rhs_file); 125 } 126 if (!m_hasRhs) 127 { 128 // Generate a random right hand side 129 if (!m_matIsLoaded) this->matrix(); 130 m_refX.resize(m_mat.cols()); 131 m_refX.setRandom(); 132 m_rhs = m_mat * m_refX; 133 m_hasrefX = true; 134 m_hasRhs = true; 135 } 136 return m_rhs; 137 } 138 139 /** Return a reference solution 140 * If it is not provided and if the right hand side is not available 141 * then refX is randomly generated such that A*refX = b 142 * where A and b are the matrix and the rhs. 143 * Note that when a rhs is provided, refX is not available 144 */ 145 inline VectorType& refX() 146 { 147 // Check if a reference solution is provided 148 if (m_hasrefX) return m_refX; 149 150 std::string lhs_file; 151 lhs_file = m_folder + "/" + m_matname + "_x.mtx"; 152 m_hasrefX = Fileexists(lhs_file); 153 if (m_hasrefX) 154 { 155 m_refX.resize(m_mat.cols()); 156 m_hasrefX = loadMarketVector(m_refX, lhs_file); 157 } 158 else 159 m_refX.resize(0); 160 return m_refX; 161 } 162 163 inline std::string& matname() { return m_matname; } 164 165 inline int sym() { return m_sym; } 166 167 bool hasRhs() {return m_hasRhs; } 168 bool hasrefX() {return m_hasrefX; } 169 bool isFolderValid() { return bool(m_folder_id); } 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=false; 203 if(!getMarketHeader(curfile,m_sym,iscomplex,isvector)) continue; 204 if(isvector) continue; 205 if (!iscomplex) 206 { 207 if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value) 208 continue; 209 } 210 if (iscomplex) 211 { 212 if(internal::is_same<Scalar, float>::value || internal::is_same<Scalar, double>::value) 213 continue; 214 } 215 216 217 // Get the matrix name 218 std::string filename = m_curs_id->d_name; 219 m_matname = filename.substr(0, filename.length()-4); 220 221 // Find if the matrix is SPD 222 size_t found = m_matname.find("SPD"); 223 if( (found!=std::string::npos) && (m_sym != NonSymmetric) ) 224 m_sym = SPD; 225 226 m_isvalid = true; 227 break; 228 } 229 } 230 int m_sym; // Symmetry of the matrix 231 MatrixType m_mat; // Current matrix 232 VectorType m_rhs; // Current vector 233 VectorType m_refX; // The reference solution, if exists 234 std::string m_matname; // Matrix Name 235 bool m_isvalid; 236 bool m_matIsLoaded; // Determine if the matrix has already been loaded from the file 237 bool m_hasRhs; // The right hand side exists 238 bool m_hasrefX; // A reference solution is provided 239 std::string m_folder; 240 DIR * m_folder_id; 241 struct dirent *m_curs_id; 242 243 }; 244 245 } // end namespace Eigen 246 247 #endif 248
