MatrixMarketIterator.h
00001 
00002 // This file is part of Eigen, a lightweight C++ template library
00003 // for linear algebra.
00004 //
00005 // Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam@inria.fr>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_BROWSE_MATRICES_H
00012 #define EIGEN_BROWSE_MATRICES_H
00013 
00014 namespace Eigen {
00015 
00016 enum {
00017   SPD = 0x100,
00018   NonSymmetric = 0x0
00019 }; 
00020 
00041 template <typename Scalar>
00042 class MatrixMarketIterator 
00043 {
00044     typedef typename NumTraits<Scalar>::Real RealScalar;
00045   public:
00046     typedef Matrix<Scalar,Dynamic,1> VectorType; 
00047     typedef SparseMatrix<Scalar,ColMajor> MatrixType; 
00048   
00049   public:
00050     MatrixMarketIterator(const std::string &folder)
00051       : m_sym(0), m_isvalid(false), m_matIsLoaded(false), m_hasRhs(false), m_hasrefX(false), m_folder(folder)
00052     {
00053       m_folder_id = opendir(folder.c_str());
00054       if(m_folder_id)
00055         Getnextvalidmatrix();
00056     }
00057     
00058     ~MatrixMarketIterator()
00059     {
00060       if (m_folder_id) closedir(m_folder_id); 
00061     }
00062     
00063     inline MatrixMarketIterator& operator++()
00064     {
00065       m_matIsLoaded = false;
00066       m_hasrefX = false;
00067       m_hasRhs = false;
00068       Getnextvalidmatrix();
00069       return *this;
00070     }
00071     inline operator bool() const { return m_isvalid;}
00072     
00074     inline MatrixType& matrix() 
00075     { 
00076       // Read the matrix
00077       if (m_matIsLoaded) return m_mat;
00078       
00079       std::string matrix_file = m_folder + "/" + m_matname + ".mtx";
00080       if ( !loadMarket(m_mat, matrix_file)) 
00081       {
00082         std::cerr << "Warning loadMarket failed when loading \"" << matrix_file << "\"" << std::endl;
00083         m_matIsLoaded = false;
00084         return m_mat;
00085       }
00086       m_matIsLoaded = true; 
00087 
00088       if (m_sym != NonSymmetric) 
00089       {
00090         // Check whether we need to restore a full matrix:
00091         RealScalar diag_norm  = m_mat.diagonal().norm();
00092         RealScalar lower_norm = m_mat.template triangularView<Lower>().norm();
00093         RealScalar upper_norm = m_mat.template triangularView<Upper>().norm();
00094         if(lower_norm>diag_norm && upper_norm==diag_norm)
00095         {
00096           // only the lower part is stored
00097           MatrixType tmp(m_mat);
00098           m_mat = tmp.template selfadjointView<Lower>();
00099         }
00100         else if(upper_norm>diag_norm && lower_norm==diag_norm)
00101         {
00102           // only the upper part is stored
00103           MatrixType tmp(m_mat);
00104           m_mat = tmp.template selfadjointView<Upper>();
00105         }
00106       }
00107       return m_mat; 
00108     }
00109     
00113     inline VectorType& rhs() 
00114     { 
00115        // Get the right hand side
00116       if (m_hasRhs) return m_rhs;
00117       
00118       std::string rhs_file;
00119       rhs_file = m_folder + "/" + m_matname + "_b.mtx"; // The pattern is matname_b.mtx
00120       m_hasRhs = Fileexists(rhs_file);
00121       if (m_hasRhs)
00122       {
00123         m_rhs.resize(m_mat.cols());
00124         m_hasRhs = loadMarketVector(m_rhs, rhs_file);
00125       }
00126       if (!m_hasRhs)
00127       {
00128         // Generate a random right hand side
00129         if (!m_matIsLoaded) this->matrix(); 
00130         m_refX.resize(m_mat.cols());
00131         m_refX.setRandom();
00132         m_rhs = m_mat * m_refX;
00133         m_hasrefX = true;
00134         m_hasRhs = true;
00135       }
00136       return m_rhs; 
00137     }
00138     
00145     inline VectorType& refX() 
00146     { 
00147       // Check if a reference solution is provided
00148       if (m_hasrefX) return m_refX;
00149       
00150       std::string lhs_file;
00151       lhs_file = m_folder + "/" + m_matname + "_x.mtx"; 
00152       m_hasrefX = Fileexists(lhs_file);
00153       if (m_hasrefX)
00154       {
00155         m_refX.resize(m_mat.cols());
00156         m_hasrefX = loadMarketVector(m_refX, lhs_file);
00157       }
00158       else
00159         m_refX.resize(0);
00160       return m_refX; 
00161     }
00162     
00163     inline std::string& matname() { return m_matname; }
00164     
00165     inline int sym() { return m_sym; }
00166     
00167     bool hasRhs() {return m_hasRhs; }
00168     bool hasrefX() {return m_hasrefX; }
00169     bool isFolderValid() { return bool(m_folder_id); }
00170     
00171   protected:
00172     
00173     inline bool Fileexists(std::string file)
00174     {
00175       std::ifstream file_id(file.c_str());
00176       if (!file_id.good() ) 
00177       {
00178         return false;
00179       }
00180       else 
00181       {
00182         file_id.close();
00183         return true;
00184       }
00185     }
00186     
00187     void Getnextvalidmatrix( )
00188     {
00189       m_isvalid = false;
00190       // Here, we return with the next valid matrix in the folder
00191       while ( (m_curs_id = readdir(m_folder_id)) != NULL) {
00192         m_isvalid = false;
00193         std::string curfile;
00194         curfile = m_folder + "/" + m_curs_id->d_name;
00195         // Discard if it is a folder
00196         if (m_curs_id->d_type == DT_DIR) continue; //FIXME This may not be available on non BSD systems
00197 //         struct stat st_buf; 
00198 //         stat (curfile.c_str(), &st_buf);
00199 //         if (S_ISDIR(st_buf.st_mode)) continue;
00200         
00201         // Determine from the header if it is a matrix or a right hand side 
00202         bool isvector,iscomplex=false;
00203         if(!getMarketHeader(curfile,m_sym,iscomplex,isvector)) continue;
00204         if(isvector) continue;
00205         if (!iscomplex)
00206         {
00207           if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value)
00208             continue; 
00209         }
00210         if (iscomplex)
00211         {
00212           if(internal::is_same<Scalar, float>::value || internal::is_same<Scalar, double>::value)
00213             continue; 
00214         }
00215         
00216         
00217         // Get the matrix name
00218         std::string filename = m_curs_id->d_name;
00219         m_matname = filename.substr(0, filename.length()-4); 
00220         
00221         // Find if the matrix is SPD 
00222         size_t found = m_matname.find("SPD");
00223         if( (found!=std::string::npos) && (m_sym != NonSymmetric) )
00224           m_sym = SPD;
00225        
00226         m_isvalid = true;
00227         break; 
00228       }
00229     }
00230     int m_sym; // Symmetry of the matrix
00231     MatrixType m_mat; // Current matrix  
00232     VectorType m_rhs;  // Current vector
00233     VectorType m_refX; // The reference solution, if exists
00234     std::string m_matname; // Matrix Name
00235     bool m_isvalid; 
00236     bool m_matIsLoaded; // Determine if the matrix has already been loaded from the file
00237     bool m_hasRhs; // The right hand side exists
00238     bool m_hasrefX; // A reference solution is provided
00239     std::string m_folder;
00240     DIR * m_folder_id;
00241     struct dirent *m_curs_id; 
00242     
00243 };
00244 
00245 } // end namespace Eigen
00246 
00247 #endif
 All Classes Functions Variables Typedefs Enumerator