![]() |
Eigen-unsupported
3.3.3
|
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, Index& M, Index& N, Index& i, Index& 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, Index& M, Index& N, Index& i, Index& 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 typedef typename SparseMatrixType::Index Index; 00137 std::ifstream input(filename.c_str(),std::ios::in); 00138 if(!input) 00139 return false; 00140 00141 const int maxBuffersize = 2048; 00142 char buffer[maxBuffersize]; 00143 00144 bool readsizes = false; 00145 00146 typedef Triplet<Scalar,Index> T; 00147 std::vector<T> elements; 00148 00149 Index M(-1), N(-1), NNZ(-1); 00150 Index count = 0; 00151 while(input.getline(buffer, maxBuffersize)) 00152 { 00153 // skip comments 00154 //NOTE An appropriate test should be done on the header to get the symmetry 00155 if(buffer[0]=='%') 00156 continue; 00157 00158 std::stringstream line(buffer); 00159 00160 if(!readsizes) 00161 { 00162 line >> M >> N >> NNZ; 00163 if(M > 0 && N > 0 && NNZ > 0) 00164 { 00165 readsizes = true; 00166 //std::cout << "sizes: " << M << "," << N << "," << NNZ << "\n"; 00167 mat.resize(M,N); 00168 mat.reserve(NNZ); 00169 } 00170 } 00171 else 00172 { 00173 Index i(-1), j(-1); 00174 Scalar value; 00175 if( internal::GetMarketLine(line, M, N, i, j, value) ) 00176 { 00177 ++ count; 00178 elements.push_back(T(i,j,value)); 00179 } 00180 else 00181 std::cerr << "Invalid read: " << i << "," << j << "\n"; 00182 } 00183 } 00184 mat.setFromTriplets(elements.begin(), elements.end()); 00185 if(count!=NNZ) 00186 std::cerr << count << "!=" << NNZ << "\n"; 00187 00188 input.close(); 00189 return true; 00190 } 00191 00192 template<typename VectorType> 00193 bool loadMarketVector(VectorType& vec, const std::string& filename) 00194 { 00195 typedef typename VectorType::Scalar Scalar; 00196 std::ifstream in(filename.c_str(), std::ios::in); 00197 if(!in) 00198 return false; 00199 00200 std::string line; 00201 int n(0), col(0); 00202 do 00203 { // Skip comments 00204 std::getline(in, line); eigen_assert(in.good()); 00205 } while (line[0] == '%'); 00206 std::istringstream newline(line); 00207 newline >> n >> col; 00208 eigen_assert(n>0 && col>0); 00209 vec.resize(n); 00210 int i = 0; 00211 Scalar value; 00212 while ( std::getline(in, line) && (i < n) ){ 00213 internal::GetVectorElt(line, value); 00214 vec(i++) = value; 00215 } 00216 in.close(); 00217 if (i!=n){ 00218 std::cerr<< "Unable to read all elements from file " << filename << "\n"; 00219 return false; 00220 } 00221 return true; 00222 } 00223 00224 template<typename SparseMatrixType> 00225 bool saveMarket(const SparseMatrixType& mat, const std::string& filename, int sym = 0) 00226 { 00227 typedef typename SparseMatrixType::Scalar Scalar; 00228 std::ofstream out(filename.c_str(),std::ios::out); 00229 if(!out) 00230 return false; 00231 00232 out.flags(std::ios_base::scientific); 00233 out.precision(64); 00234 std::string header; 00235 internal::putMarketHeader<Scalar>(header, sym); 00236 out << header << std::endl; 00237 out << mat.rows() << " " << mat.cols() << " " << mat.nonZeros() << "\n"; 00238 int count = 0; 00239 for(int j=0; j<mat.outerSize(); ++j) 00240 for(typename SparseMatrixType::InnerIterator it(mat,j); it; ++it) 00241 { 00242 ++ count; 00243 internal::PutMatrixElt(it.value(), it.row()+1, it.col()+1, out); 00244 // out << it.row()+1 << " " << it.col()+1 << " " << it.value() << "\n"; 00245 } 00246 out.close(); 00247 return true; 00248 } 00249 00250 template<typename VectorType> 00251 bool saveMarketVector (const VectorType& vec, const std::string& filename) 00252 { 00253 typedef typename VectorType::Scalar Scalar; 00254 std::ofstream out(filename.c_str(),std::ios::out); 00255 if(!out) 00256 return false; 00257 00258 out.flags(std::ios_base::scientific); 00259 out.precision(64); 00260 if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value) 00261 out << "%%MatrixMarket matrix array complex general\n"; 00262 else 00263 out << "%%MatrixMarket matrix array real general\n"; 00264 out << vec.size() << " "<< 1 << "\n"; 00265 for (int i=0; i < vec.size(); i++){ 00266 internal::putVectorElt(vec(i), out); 00267 } 00268 out.close(); 00269 return true; 00270 } 00271 00272 } // end namespace Eigen 00273 00274 #endif // EIGEN_SPARSE_MARKET_IO_H