$treeview $search $mathjax
Eigen-unsupported
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> 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_SPARSE_MARKET_IO_H 00012 #define EIGEN_SPARSE_MARKET_IO_H 00013 00014 #include <iostream> 00015 00016 namespace Eigen { 00017 00018 namespace internal 00019 { 00020 template <typename Scalar> 00021 inline bool GetMarketLine (std::stringstream& line, int& M, int& N, int& i, int& j, Scalar& value) 00022 { 00023 line >> i >> j >> value; 00024 i--; 00025 j--; 00026 if(i>=0 && j>=0 && i<M && j<N) 00027 { 00028 return true; 00029 } 00030 else 00031 return false; 00032 } 00033 template <typename Scalar> 00034 inline bool GetMarketLine (std::stringstream& line, int& M, int& N, int& i, int& j, std::complex<Scalar>& value) 00035 { 00036 Scalar valR, valI; 00037 line >> i >> j >> valR >> valI; 00038 i--; 00039 j--; 00040 if(i>=0 && j>=0 && i<M && j<N) 00041 { 00042 value = std::complex<Scalar>(valR, valI); 00043 return true; 00044 } 00045 else 00046 return false; 00047 } 00048 00049 template <typename RealScalar> 00050 inline void GetVectorElt (const std::string& line, RealScalar& val) 00051 { 00052 std::istringstream newline(line); 00053 newline >> val; 00054 } 00055 00056 template <typename RealScalar> 00057 inline void GetVectorElt (const std::string& line, std::complex<RealScalar>& val) 00058 { 00059 RealScalar valR, valI; 00060 std::istringstream newline(line); 00061 newline >> valR >> valI; 00062 val = std::complex<RealScalar>(valR, valI); 00063 } 00064 00065 template<typename Scalar> 00066 inline void putMarketHeader(std::string& header,int sym) 00067 { 00068 header= "%%MatrixMarket matrix coordinate "; 00069 if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value) 00070 { 00071 header += " complex"; 00072 if(sym == Symmetric) header += " symmetric"; 00073 else if (sym == SelfAdjoint) header += " Hermitian"; 00074 else header += " general"; 00075 } 00076 else 00077 { 00078 header += " real"; 00079 if(sym == Symmetric) header += " symmetric"; 00080 else header += " general"; 00081 } 00082 } 00083 00084 template<typename Scalar> 00085 inline void PutMatrixElt(Scalar value, int row, int col, std::ofstream& out) 00086 { 00087 out << row << " "<< col << " " << value << "\n"; 00088 } 00089 template<typename Scalar> 00090 inline void PutMatrixElt(std::complex<Scalar> value, int row, int col, std::ofstream& out) 00091 { 00092 out << row << " " << col << " " << value.real() << " " << value.imag() << "\n"; 00093 } 00094 00095 00096 template<typename Scalar> 00097 inline void putVectorElt(Scalar value, std::ofstream& out) 00098 { 00099 out << value << "\n"; 00100 } 00101 template<typename Scalar> 00102 inline void putVectorElt(std::complex<Scalar> value, std::ofstream& out) 00103 { 00104 out << value.real << " " << value.imag()<< "\n"; 00105 } 00106 00107 } // end namepsace internal 00108 00109 inline bool getMarketHeader(const std::string& filename, int& sym, bool& iscomplex, bool& isvector) 00110 { 00111 sym = 0; 00112 isvector = false; 00113 std::ifstream in(filename.c_str(),std::ios::in); 00114 if(!in) 00115 return false; 00116 00117 std::string line; 00118 // The matrix header is always the first line in the file 00119 std::getline(in, line); eigen_assert(in.good()); 00120 00121 std::stringstream fmtline(line); 00122 std::string substr[5]; 00123 fmtline>> substr[0] >> substr[1] >> substr[2] >> substr[3] >> substr[4]; 00124 if(substr[2].compare("array") == 0) isvector = true; 00125 if(substr[3].compare("complex") == 0) iscomplex = true; 00126 if(substr[4].compare("symmetric") == 0) sym = Symmetric; 00127 else if (substr[4].compare("Hermitian") == 0) sym = SelfAdjoint; 00128 00129 return true; 00130 } 00131 00132 template<typename SparseMatrixType> 00133 bool loadMarket(SparseMatrixType& mat, const std::string& filename) 00134 { 00135 typedef typename SparseMatrixType::Scalar Scalar; 00136 std::ifstream input(filename.c_str(),std::ios::in); 00137 if(!input) 00138 return false; 00139 00140 const int maxBuffersize = 2048; 00141 char buffer[maxBuffersize]; 00142 00143 bool readsizes = false; 00144 00145 typedef Triplet<Scalar,int> T; 00146 std::vector<T> elements; 00147 00148 int M(-1), N(-1), NNZ(-1); 00149 int count = 0; 00150 while(input.getline(buffer, maxBuffersize)) 00151 { 00152 // skip comments 00153 //NOTE An appropriate test should be done on the header to get the symmetry 00154 if(buffer[0]=='%') 00155 continue; 00156 00157 std::stringstream line(buffer); 00158 00159 if(!readsizes) 00160 { 00161 line >> M >> N >> NNZ; 00162 if(M > 0 && N > 0 && NNZ > 0) 00163 { 00164 readsizes = true; 00165 std::cout << "sizes: " << M << "," << N << "," << NNZ << "\n"; 00166 mat.resize(M,N); 00167 mat.reserve(NNZ); 00168 } 00169 } 00170 else 00171 { 00172 int i(-1), j(-1); 00173 Scalar value; 00174 if( internal::GetMarketLine(line, M, N, i, j, value) ) 00175 { 00176 ++ count; 00177 elements.push_back(T(i,j,value)); 00178 } 00179 else 00180 std::cerr << "Invalid read: " << i << "," << j << "\n"; 00181 } 00182 } 00183 mat.setFromTriplets(elements.begin(), elements.end()); 00184 if(count!=NNZ) 00185 std::cerr << count << "!=" << NNZ << "\n"; 00186 00187 input.close(); 00188 return true; 00189 } 00190 00191 template<typename VectorType> 00192 bool loadMarketVector(VectorType& vec, const std::string& filename) 00193 { 00194 typedef typename VectorType::Scalar Scalar; 00195 std::ifstream in(filename.c_str(), std::ios::in); 00196 if(!in) 00197 return false; 00198 00199 std::string line; 00200 int n(0), col(0); 00201 do 00202 { // Skip comments 00203 std::getline(in, line); eigen_assert(in.good()); 00204 } while (line[0] == '%'); 00205 std::istringstream newline(line); 00206 newline >> n >> col; 00207 eigen_assert(n>0 && col>0); 00208 vec.resize(n); 00209 int i = 0; 00210 Scalar value; 00211 while ( std::getline(in, line) && (i < n) ){ 00212 internal::GetVectorElt(line, value); 00213 vec(i++) = value; 00214 } 00215 in.close(); 00216 if (i!=n){ 00217 std::cerr<< "Unable to read all elements from file " << filename << "\n"; 00218 return false; 00219 } 00220 return true; 00221 } 00222 00223 template<typename SparseMatrixType> 00224 bool saveMarket(const SparseMatrixType& mat, const std::string& filename, int sym = 0) 00225 { 00226 typedef typename SparseMatrixType::Scalar Scalar; 00227 std::ofstream out(filename.c_str(),std::ios::out); 00228 if(!out) 00229 return false; 00230 00231 out.flags(std::ios_base::scientific); 00232 out.precision(64); 00233 std::string header; 00234 internal::putMarketHeader<Scalar>(header, sym); 00235 out << header << std::endl; 00236 out << mat.rows() << " " << mat.cols() << " " << mat.nonZeros() << "\n"; 00237 int count = 0; 00238 for(int j=0; j<mat.outerSize(); ++j) 00239 for(typename SparseMatrixType::InnerIterator it(mat,j); it; ++it) 00240 { 00241 ++ count; 00242 internal::PutMatrixElt(it.value(), it.row()+1, it.col()+1, out); 00243 // out << it.row()+1 << " " << it.col()+1 << " " << it.value() << "\n"; 00244 } 00245 out.close(); 00246 return true; 00247 } 00248 00249 template<typename VectorType> 00250 bool saveMarketVector (const VectorType& vec, const std::string& filename) 00251 { 00252 typedef typename VectorType::Scalar Scalar; 00253 std::ofstream out(filename.c_str(),std::ios::out); 00254 if(!out) 00255 return false; 00256 00257 out.flags(std::ios_base::scientific); 00258 out.precision(64); 00259 if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value) 00260 out << "%%MatrixMarket matrix array complex general\n"; 00261 else 00262 out << "%%MatrixMarket matrix array real general\n"; 00263 out << vec.size() << " "<< 1 << "\n"; 00264 for (int i=0; i < vec.size(); i++){ 00265 internal::putVectorElt(vec(i), out); 00266 } 00267 out.close(); 00268 return true; 00269 } 00270 00271 } // end namespace Eigen 00272 00273 #endif // EIGEN_SPARSE_MARKET_IO_H