src/SparseExtra/MatrixMarketIterator.h

0001
0002 // This file is part of Eigen, a lightweight C++ template library
0003 // for linear algebra.
0004 //
0005 // Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam@inria.fr>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010
0011 #ifndef EIGEN_BROWSE_MATRICES_H
0012 #define EIGEN_BROWSE_MATRICES_H
0013
0014 namespace Eigen {
0015
0016 enum {
0017   SPD = 0x100,
0018   NonSymmetric = 0x0
0019 };
0020
0021 /**
0022  * @brief Iterator to browse matrices from a specified folder
0023  *
0024  * This is used to load all the matrices from a folder.
0025  * The matrices should be in Matrix Market format
0026  * It is assumed that the matrices are named as matname.mtx
0027  * and matname_SPD.mtx if the matrix is Symmetric and positive definite (or Hermitian)
0028  * The right hand side vectors are loaded as well, if they exist.
0029  * They should be named as matname_b.mtx.
0030  * Note that the right hand side for a SPD matrix is named as matname_SPD_b.mtx
0031  *
0032  * Sometimes a reference solution is available. In this case, it should be named as matname_x.mtx
0033  *
0034  * Sample code
0035  * \code
0036  *
0037  * \endcode
0038  *
0039  * \tparam Scalar The scalar type
0040  */
0041 template <typename Scalar>
0042 class MatrixMarketIterator
0043 {
0044     typedef typename NumTraits<Scalar>::Real RealScalar;
0045   public:
0046     typedef Matrix<Scalar,Dynamic,1> VectorType;
0047     typedef SparseMatrix<Scalar,ColMajor> MatrixType;
0048
0049   public:
0050     MatrixMarketIterator(const std::string &folder)
0051       : m_sym(0), m_isvalid(false), m_matIsLoaded(false), m_hasRhs(false), m_hasrefX(false), m_folder(folder)
0052     {
0053       m_folder_id = opendir(folder.c_str());
0054       if(m_folder_id)
0055         Getnextvalidmatrix();
0056     }
0057
0058     ~MatrixMarketIterator()
0059     {
0060       if (m_folder_id) closedir(m_folder_id);
0061     }
0062
0063     inline MatrixMarketIterator& operator++()
0064     {
0065       m_matIsLoaded = false;
0066       m_hasrefX = false;
0067       m_hasRhs = false;
0068       Getnextvalidmatrix();
0069       return *this;
0070     }
0071     inline operator bool() const { return m_isvalid;}
0072
0073     /** Return the sparse matrix corresponding to the current file */
0074     inline MatrixType& matrix()
0075     {
0076       // Read the matrix
0077       if (m_matIsLoaded) return m_mat;
0078
0079       std::string matrix_file = m_folder + "/" + m_matname + ".mtx";
0080       if ( !loadMarket(m_mat, matrix_file))
0081       {
0082         std::cerr << "Warning loadMarket failed when loading \"" << matrix_file << "\"" << std::endl;
0083         m_matIsLoaded = false;
0084         return m_mat;
0085       }
0086       m_matIsLoaded = true;
0087
0088       if (m_sym != NonSymmetric)
0089       {
0090         // Check whether we need to restore a full matrix:
0091         RealScalar diag_norm  = m_mat.diagonal().norm();
0092         RealScalar lower_norm = m_mat.template triangularView<Lower>().norm();
0093         RealScalar upper_norm = m_mat.template triangularView<Upper>().norm();
0094         if(lower_norm>diag_norm && upper_norm==diag_norm)
0095         {
0096           // only the lower part is stored
0097           MatrixType tmp(m_mat);
0098           m_mat = tmp.template selfadjointView<Lower>();
0099         }
0100         else if(upper_norm>diag_norm && lower_norm==diag_norm)
0101         {
0102           // only the upper part is stored
0103           MatrixType tmp(m_mat);
0104           m_mat = tmp.template selfadjointView<Upper>();
0105         }
0106       }
0107       return m_mat;
0108     }
0109
0110     /** Return the right hand side corresponding to the current matrix.
0111      * If the rhs file is not provided, a random rhs is generated
0112      */
0113     inline VectorType& rhs()
0114     {
0115        // Get the right hand side
0116       if (m_hasRhs) return m_rhs;
0117
0118       std::string rhs_file;
0119       rhs_file = m_folder + "/" + m_matname + "_b.mtx"; // The pattern is matname_b.mtx
0120       m_hasRhs = Fileexists(rhs_file);
0121       if (m_hasRhs)
0122       {
0123         m_rhs.resize(m_mat.cols());
0124         m_hasRhs = loadMarketVector(m_rhs, rhs_file);
0125       }
0126       if (!m_hasRhs)
0127       {
0128         // Generate a random right hand side
0129         if (!m_matIsLoaded) this->matrix();
0130         m_refX.resize(m_mat.cols());
0131         m_refX.setRandom();
0132         m_rhs = m_mat * m_refX;
0133         m_hasrefX = true;
0134         m_hasRhs = true;
0135       }
0136       return m_rhs;
0137     }
0138
0139     /** Return a reference solution
0140      * If it is not provided and if the right hand side is not available
0141      * then refX is randomly generated such that A*refX = b
0142      * where A and b are the matrix and the rhs.
0143      * Note that when a rhs is provided, refX is not available
0144      */
0145     inline VectorType& refX()
0146     {
0147       // Check if a reference solution is provided
0148       if (m_hasrefX) return m_refX;
0149
0150       std::string lhs_file;
0151       lhs_file = m_folder + "/" + m_matname + "_x.mtx";
0152       m_hasrefX = Fileexists(lhs_file);
0153       if (m_hasrefX)
0154       {
0155         m_refX.resize(m_mat.cols());
0156         m_hasrefX = loadMarketVector(m_refX, lhs_file);
0157       }
0158       else
0159         m_refX.resize(0);
0160       return m_refX;
0161     }
0162
0163     inline std::string& matname() { return m_matname; }
0164
0165     inline int sym() { return m_sym; }
0166
0167     bool hasRhs() {return m_hasRhs; }
0168     bool hasrefX() {return m_hasrefX; }
0169     bool isFolderValid() { return bool(m_folder_id); }
0170
0171   protected:
0172
0173     inline bool Fileexists(std::string file)
0174     {
0175       std::ifstream file_id(file.c_str());
0176       if (!file_id.good() )
0177       {
0178         return false;
0179       }
0180       else
0181       {
0182         file_id.close();
0183         return true;
0184       }
0185     }
0186
0187     void Getnextvalidmatrix( )
0188     {
0189       m_isvalid = false;
0190       // Here, we return with the next valid matrix in the folder
0191       while ( (m_curs_id = readdir(m_folder_id)) != NULL) {
0192         m_isvalid = false;
0193         std::string curfile;
0194         curfile = m_folder + "/" + m_curs_id->d_name;
0195         // Discard if it is a folder
0196         if (m_curs_id->d_type == DT_DIR) continue; //FIXME This may not be available on non BSD systems
0197 //         struct stat st_buf;
0198 //         stat (curfile.c_str(), &st_buf);
0199 //         if (S_ISDIR(st_buf.st_mode)) continue;
0200
0201         // Determine from the header if it is a matrix or a right hand side
0202         bool isvector,iscomplex=false;
0203         if(!getMarketHeader(curfile,m_sym,iscomplex,isvector)) continue;
0204         if(isvector) continue;
0205         if (!iscomplex)
0206         {
0207           if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value)
0208             continue;
0209         }
0210         if (iscomplex)
0211         {
0212           if(internal::is_same<Scalar, float>::value || internal::is_same<Scalar, double>::value)
0213             continue;
0214         }
0215
0216
0217         // Get the matrix name
0218         std::string filename = m_curs_id->d_name;
0219         m_matname = filename.substr(0, filename.length()-4);
0220
0221         // Find if the matrix is SPD
0222         size_t found = m_matname.find("SPD");
0223         if( (found!=std::string::npos) && (m_sym != NonSymmetric) )
0224           m_sym = SPD;
0225
0226         m_isvalid = true;
0227         break;
0228       }
0229     }
0230     int m_sym; // Symmetry of the matrix
0231     MatrixType m_mat; // Current matrix
0232     VectorType m_rhs;  // Current vector
0233     VectorType m_refX; // The reference solution, if exists
0234     std::string m_matname; // Matrix Name
0235     bool m_isvalid;
0236     bool m_matIsLoaded; // Determine if the matrix has already been loaded from the file
0237     bool m_hasRhs; // The right hand side exists
0238     bool m_hasrefX; // A reference solution is provided
0239     std::string m_folder;
0240     DIR * m_folder_id;
0241     struct dirent *m_curs_id;
0242
0243 };
0244
0245 } // end namespace Eigen
0246
0247 #endif