// This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2012 Desire NUENTSA WAKAM // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #ifndef EIGEN_BROWSE_MATRICES_H #define EIGEN_BROWSE_MATRICES_H namespace Eigen { enum { SPD = 0x100, NonSymmetric = 0x0 }; /** * @brief Iterator to browse matrices from a specified folder * * This is used to load all the matrices from a folder. * The matrices should be in Matrix Market format * It is assumed that the matrices are named as matname.mtx * and matname_SPD.mtx if the matrix is Symmetric and positive definite (or Hermitian) * The right hand side vectors are loaded as well, if they exist. * They should be named as matname_b.mtx. * Note that the right hand side for a SPD matrix is named as matname_SPD_b.mtx * * Sometimes a reference solution is available. In this case, it should be named as matname_x.mtx * * Sample code * \code * * \endcode * * \tparam Scalar The scalar type */ template class MatrixMarketIterator { typedef typename NumTraits::Real RealScalar; public: typedef Matrix VectorType; typedef SparseMatrix MatrixType; public: MatrixMarketIterator(const std::string &folder) : m_sym(0), m_isvalid(false), m_matIsLoaded(false), m_hasRhs(false), m_hasrefX(false), m_folder(folder) { m_folder_id = opendir(folder.c_str()); if(m_folder_id) Getnextvalidmatrix(); } ~MatrixMarketIterator() { if (m_folder_id) closedir(m_folder_id); } inline MatrixMarketIterator& operator++() { m_matIsLoaded = false; m_hasrefX = false; m_hasRhs = false; Getnextvalidmatrix(); return *this; } inline operator bool() const { return m_isvalid;} /** Return the sparse matrix corresponding to the current file */ inline MatrixType& matrix() { // Read the matrix if (m_matIsLoaded) return m_mat; std::string matrix_file = m_folder + "/" + m_matname + ".mtx"; if ( !loadMarket(m_mat, matrix_file)) { std::cerr << "Warning loadMarket failed when loading \"" << matrix_file << "\"" << std::endl; m_matIsLoaded = false; return m_mat; } m_matIsLoaded = true; if (m_sym != NonSymmetric) { // Check whether we need to restore a full matrix: RealScalar diag_norm = m_mat.diagonal().norm(); RealScalar lower_norm = m_mat.template triangularView().norm(); RealScalar upper_norm = m_mat.template triangularView().norm(); if(lower_norm>diag_norm && upper_norm==diag_norm) { // only the lower part is stored MatrixType tmp(m_mat); m_mat = tmp.template selfadjointView(); } else if(upper_norm>diag_norm && lower_norm==diag_norm) { // only the upper part is stored MatrixType tmp(m_mat); m_mat = tmp.template selfadjointView(); } } return m_mat; } /** Return the right hand side corresponding to the current matrix. * If the rhs file is not provided, a random rhs is generated */ inline VectorType& rhs() { // Get the right hand side if (m_hasRhs) return m_rhs; std::string rhs_file; rhs_file = m_folder + "/" + m_matname + "_b.mtx"; // The pattern is matname_b.mtx m_hasRhs = Fileexists(rhs_file); if (m_hasRhs) { m_rhs.resize(m_mat.cols()); m_hasRhs = loadMarketVector(m_rhs, rhs_file); } if (!m_hasRhs) { // Generate a random right hand side if (!m_matIsLoaded) this->matrix(); m_refX.resize(m_mat.cols()); m_refX.setRandom(); m_rhs = m_mat * m_refX; m_hasrefX = true; m_hasRhs = true; } return m_rhs; } /** Return a reference solution * If it is not provided and if the right hand side is not available * then refX is randomly generated such that A*refX = b * where A and b are the matrix and the rhs. * Note that when a rhs is provided, refX is not available */ inline VectorType& refX() { // Check if a reference solution is provided if (m_hasrefX) return m_refX; std::string lhs_file; lhs_file = m_folder + "/" + m_matname + "_x.mtx"; m_hasrefX = Fileexists(lhs_file); if (m_hasrefX) { m_refX.resize(m_mat.cols()); m_hasrefX = loadMarketVector(m_refX, lhs_file); } else m_refX.resize(0); return m_refX; } inline std::string& matname() { return m_matname; } inline int sym() { return m_sym; } bool hasRhs() {return m_hasRhs; } bool hasrefX() {return m_hasrefX; } bool isFolderValid() { return bool(m_folder_id); } protected: inline bool Fileexists(std::string file) { std::ifstream file_id(file.c_str()); if (!file_id.good() ) { return false; } else { file_id.close(); return true; } } void Getnextvalidmatrix( ) { m_isvalid = false; // Here, we return with the next valid matrix in the folder while ( (m_curs_id = readdir(m_folder_id)) != NULL) { m_isvalid = false; std::string curfile; curfile = m_folder + "/" + m_curs_id->d_name; // Discard if it is a folder if (m_curs_id->d_type == DT_DIR) continue; //FIXME This may not be available on non BSD systems // struct stat st_buf; // stat (curfile.c_str(), &st_buf); // if (S_ISDIR(st_buf.st_mode)) continue; // Determine from the header if it is a matrix or a right hand side bool isvector,iscomplex=false; if(!getMarketHeader(curfile,m_sym,iscomplex,isvector)) continue; if(isvector) continue; if (!iscomplex) { if(internal::is_same >::value || internal::is_same >::value) continue; } if (iscomplex) { if(internal::is_same::value || internal::is_same::value) continue; } // Get the matrix name std::string filename = m_curs_id->d_name; m_matname = filename.substr(0, filename.length()-4); // Find if the matrix is SPD size_t found = m_matname.find("SPD"); if( (found!=std::string::npos) && (m_sym != NonSymmetric) ) m_sym = SPD; m_isvalid = true; break; } } int m_sym; // Symmetry of the matrix MatrixType m_mat; // Current matrix VectorType m_rhs; // Current vector VectorType m_refX; // The reference solution, if exists std::string m_matname; // Matrix Name bool m_isvalid; bool m_matIsLoaded; // Determine if the matrix has already been loaded from the file bool m_hasRhs; // The right hand side exists bool m_hasrefX; // A reference solution is provided std::string m_folder; DIR * m_folder_id; struct dirent *m_curs_id; }; } // end namespace Eigen #endif