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 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00025 
00026 #ifndef EIGEN_BROWSE_MATRICES_H
00027 #define EIGEN_BROWSE_MATRICES_H
00028 
00029 namespace Eigen {
00030 
00031 enum {
00032   SPD = 0x100,
00033   NonSymmetric = 0x0
00034 }; 
00035 
00056 template <typename Scalar>
00057 class MatrixMarketIterator 
00058 {
00059   public:
00060     typedef Matrix<Scalar,Dynamic,1> VectorType; 
00061     typedef SparseMatrix<Scalar,ColMajor> MatrixType; 
00062   
00063   public:
00064     MatrixMarketIterator(const std::string folder):m_sym(0),m_isvalid(false),m_matIsLoaded(false),m_hasRhs(false),m_hasrefX(false),m_folder(folder)
00065     {
00066       m_folder_id = opendir(folder.c_str());
00067       if (!m_folder_id){
00068         m_isvalid = false;
00069         std::cerr << "The provided Matrix folder could not be opened \n\n";
00070         abort();
00071       }
00072       Getnextvalidmatrix();
00073     }
00074     
00075     ~MatrixMarketIterator()
00076     {
00077       if (m_folder_id) closedir(m_folder_id); 
00078     }
00079     
00080     inline MatrixMarketIterator& operator++()
00081     {
00082       m_matIsLoaded = false;
00083       m_hasrefX = false;
00084       m_hasRhs = false;
00085       Getnextvalidmatrix();
00086       return *this;
00087     }
00088     inline operator bool() { return m_isvalid;}
00089     
00091     inline MatrixType& matrix() 
00092     { 
00093       // Read the matrix
00094       if (m_matIsLoaded) return m_mat;
00095       
00096       std::string matrix_file = m_folder + "/" + m_matname + ".mtx";
00097       if ( !loadMarket(m_mat, matrix_file)) 
00098       {
00099         m_matIsLoaded = false;
00100         return m_mat;
00101       }
00102       m_matIsLoaded = true; 
00103       
00104       if (m_sym != NonSymmetric) 
00105       { // Store the upper part of the matrix. It is needed by the solvers dealing with nonsymmetric matrices ??
00106         MatrixType B; 
00107         B = m_mat;
00108         m_mat = B.template selfadjointView<Lower>();
00109       }
00110       return m_mat; 
00111     }
00112     
00116     inline VectorType& rhs() 
00117     { 
00118        // Get the right hand side
00119       if (m_hasRhs) return m_rhs;
00120       
00121       std::string rhs_file;
00122       rhs_file = m_folder + "/" + m_matname + "_b.mtx"; // The pattern is matname_b.mtx
00123       m_hasRhs = Fileexists(rhs_file);
00124       if (m_hasRhs)
00125       {
00126         m_rhs.resize(m_mat.cols());
00127         m_hasRhs = loadMarketVector(m_rhs, rhs_file);
00128       }
00129       if (!m_hasRhs)
00130       {
00131         // Generate a random right hand side
00132         if (!m_matIsLoaded) this->matrix(); 
00133         m_refX.resize(m_mat.cols());
00134         m_refX.setRandom();
00135         m_rhs = m_mat * m_refX;
00136         m_hasrefX = true;
00137         m_hasRhs = true;
00138       }
00139       return m_rhs; 
00140     }
00141     
00148     inline VectorType& refX() 
00149     { 
00150       // Check if a reference solution is provided
00151       if (m_hasrefX) return m_refX;
00152       
00153       std::string lhs_file;
00154       lhs_file = m_folder + "/" + m_matname + "_x.mtx"; 
00155       m_hasrefX = Fileexists(lhs_file);
00156       if (m_hasrefX)
00157       {
00158         m_refX.resize(m_mat.cols());
00159         m_hasrefX = loadMarketVector(m_refX, lhs_file);
00160       }
00161       return m_refX; 
00162     }
00163     
00164     inline std::string& matname() { return m_matname; }
00165     
00166     inline int sym() { return m_sym; }
00167     
00168     inline bool hasRhs() {return m_hasRhs; }
00169     inline bool hasrefX() {return m_hasrefX; }
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       // Here, we return with the next valid matrix in the folder
00190       while ( (m_curs_id = readdir(m_folder_id)) != NULL) {
00191         m_isvalid = false;
00192         std::string curfile;
00193         curfile = m_folder + "/" + m_curs_id->d_name;
00194         // Discard if it is a folder
00195         if (m_curs_id->d_type == DT_DIR) continue; //FIXME This may not be available on non BSD systems
00196 //         struct stat st_buf; 
00197 //         stat (curfile.c_str(), &st_buf);
00198 //         if (S_ISDIR(st_buf.st_mode)) continue;
00199         
00200         // Determine from the header if it is a matrix or a right hand side 
00201         bool isvector,iscomplex;
00202         if(!getMarketHeader(curfile,m_sym,iscomplex,isvector)) continue;
00203         if(isvector) continue;
00204         
00205         // Get the matrix name
00206         std::string filename = m_curs_id->d_name;
00207         m_matname = filename.substr(0, filename.length()-4); 
00208         
00209         // Find if the matrix is SPD 
00210         size_t found = m_matname.find("SPD");
00211         if( (found!=std::string::npos) && (m_sym != NonSymmetric) )
00212           m_sym = SPD;
00213        
00214         m_isvalid = true;
00215         break; 
00216       }
00217     }
00218     int m_sym; // Symmetry of the matrix
00219     MatrixType m_mat; // Current matrix  
00220     VectorType m_rhs;  // Current vector
00221     VectorType m_refX; // The reference solution, if exists
00222     std::string m_matname; // Matrix Name
00223     bool m_isvalid; 
00224     bool m_matIsLoaded; // Determine if the matrix has already been loaded from the file
00225     bool m_hasRhs; // The right hand side exists
00226     bool m_hasrefX; // A reference solution is provided
00227     std::string m_folder;
00228     DIR * m_folder_id;
00229     struct dirent *m_curs_id; 
00230     
00231 };
00232 
00233 } // end namespace Eigen
00234 
00235 #endif