![]() |
Eigen-unsupported
3.3.3
|
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