1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud (at) inria.fr> 5 // Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam (at) inria.fr> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #ifndef EIGEN_SPARSE_MARKET_IO_H 12 #define EIGEN_SPARSE_MARKET_IO_H 13 14 #include <iostream> 15 16 namespace Eigen { 17 18 namespace internal 19 { 20 template <typename Scalar> 21 inline bool GetMarketLine (std::stringstream& line, int& M, int& N, int& i, int& j, Scalar& value) 22 { 23 line >> i >> j >> value; 24 i--; 25 j--; 26 if(i>=0 && j>=0 && i<M && j<N) 27 { 28 return true; 29 } 30 else 31 return false; 32 } 33 template <typename Scalar> 34 inline bool GetMarketLine (std::stringstream& line, int& M, int& N, int& i, int& j, std::complex<Scalar>& value) 35 { 36 Scalar valR, valI; 37 line >> i >> j >> valR >> valI; 38 i--; 39 j--; 40 if(i>=0 && j>=0 && i<M && j<N) 41 { 42 value = std::complex<Scalar>(valR, valI); 43 return true; 44 } 45 else 46 return false; 47 } 48 49 template <typename RealScalar> 50 inline void GetVectorElt (const std::string& line, RealScalar& val) 51 { 52 std::istringstream newline(line); 53 newline >> val; 54 } 55 56 template <typename RealScalar> 57 inline void GetVectorElt (const std::string& line, std::complex<RealScalar>& val) 58 { 59 RealScalar valR, valI; 60 std::istringstream newline(line); 61 newline >> valR >> valI; 62 val = std::complex<RealScalar>(valR, valI); 63 } 64 65 template<typename Scalar> 66 inline void putMarketHeader(std::string& header,int sym) 67 { 68 header= "%%MatrixMarket matrix coordinate "; 69 if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value) 70 { 71 header += " complex"; 72 if(sym == Symmetric) header += " symmetric"; 73 else if (sym == SelfAdjoint) header += " Hermitian"; 74 else header += " general"; 75 } 76 else 77 { 78 header += " real"; 79 if(sym == Symmetric) header += " symmetric"; 80 else header += " general"; 81 } 82 } 83 84 template<typename Scalar> 85 inline void PutMatrixElt(Scalar value, int row, int col, std::ofstream& out) 86 { 87 out << row << " "<< col << " " << value << "\n"; 88 } 89 template<typename Scalar> 90 inline void PutMatrixElt(std::complex<Scalar> value, int row, int col, std::ofstream& out) 91 { 92 out << row << " " << col << " " << value.real() << " " << value.imag() << "\n"; 93 } 94 95 96 template<typename Scalar> 97 inline void putVectorElt(Scalar value, std::ofstream& out) 98 { 99 out << value << "\n"; 100 } 101 template<typename Scalar> 102 inline void putVectorElt(std::complex<Scalar> value, std::ofstream& out) 103 { 104 out << value.real << " " << value.imag()<< "\n"; 105 } 106 107 } // end namepsace internal 108 109 inline bool getMarketHeader(const std::string& filename, int& sym, bool& iscomplex, bool& isvector) 110 { 111 sym = 0; 112 isvector = false; 113 std::ifstream in(filename.c_str(),std::ios::in); 114 if(!in) 115 return false; 116 117 std::string line; 118 // The matrix header is always the first line in the file 119 std::getline(in, line); assert(in.good()); 120 121 std::stringstream fmtline(line); 122 std::string substr[5]; 123 fmtline>> substr[0] >> substr[1] >> substr[2] >> substr[3] >> substr[4]; 124 if(substr[2].compare("array") == 0) isvector = true; 125 if(substr[3].compare("complex") == 0) iscomplex = true; 126 if(substr[4].compare("symmetric") == 0) sym = Symmetric; 127 else if (substr[4].compare("Hermitian") == 0) sym = SelfAdjoint; 128 129 return true; 130 } 131 132 template<typename SparseMatrixType> 133 bool loadMarket(SparseMatrixType& mat, const std::string& filename) 134 { 135 typedef typename SparseMatrixType::Scalar Scalar; 136 std::ifstream input(filename.c_str(),std::ios::in); 137 if(!input) 138 return false; 139 140 const int maxBuffersize = 2048; 141 char buffer[maxBuffersize]; 142 143 bool readsizes = false; 144 145 typedef Triplet<Scalar,int> T; 146 std::vector<T> elements; 147 148 int M(-1), N(-1), NNZ(-1); 149 int count = 0; 150 while(input.getline(buffer, maxBuffersize)) 151 { 152 // skip comments 153 //NOTE An appropriate test should be done on the header to get the symmetry 154 if(buffer[0]=='%') 155 continue; 156 157 std::stringstream line(buffer); 158 159 if(!readsizes) 160 { 161 line >> M >> N >> NNZ; 162 if(M > 0 && N > 0 && NNZ > 0) 163 { 164 readsizes = true; 165 std::cout << "sizes: " << M << "," << N << "," << NNZ << "\n"; 166 mat.resize(M,N); 167 mat.reserve(NNZ); 168 } 169 } 170 else 171 { 172 int i(-1), j(-1); 173 Scalar value; 174 if( internal::GetMarketLine(line, M, N, i, j, value) ) 175 { 176 ++ count; 177 elements.push_back(T(i,j,value)); 178 } 179 else 180 std::cerr << "Invalid read: " << i << "," << j << "\n"; 181 } 182 } 183 mat.setFromTriplets(elements.begin(), elements.end()); 184 if(count!=NNZ) 185 std::cerr << count << "!=" << NNZ << "\n"; 186 187 input.close(); 188 return true; 189 } 190 191 template<typename VectorType> 192 bool loadMarketVector(VectorType& vec, const std::string& filename) 193 { 194 typedef typename VectorType::Scalar Scalar; 195 std::ifstream in(filename.c_str(), std::ios::in); 196 if(!in) 197 return false; 198 199 std::string line; 200 int n(0), col(0); 201 do 202 { // Skip comments 203 std::getline(in, line); assert(in.good()); 204 } while (line[0] == '%'); 205 std::istringstream newline(line); 206 newline >> n >> col; 207 assert(n>0 && col>0); 208 vec.resize(n); 209 int i = 0; 210 Scalar value; 211 while ( std::getline(in, line) && (i < n) ){ 212 internal::GetVectorElt(line, value); 213 vec(i++) = value; 214 } 215 in.close(); 216 if (i!=n){ 217 std::cerr<< "Unable to read all elements from file " << filename << "\n"; 218 return false; 219 } 220 return true; 221 } 222 223 template<typename SparseMatrixType> 224 bool saveMarket(const SparseMatrixType& mat, const std::string& filename, int sym = 0) 225 { 226 typedef typename SparseMatrixType::Scalar Scalar; 227 std::ofstream out(filename.c_str(),std::ios::out); 228 if(!out) 229 return false; 230 231 out.flags(std::ios_base::scientific); 232 out.precision(64); 233 std::string header; 234 internal::putMarketHeader<Scalar>(header, sym); 235 out << header << std::endl; 236 out << mat.rows() << " " << mat.cols() << " " << mat.nonZeros() << "\n"; 237 int count = 0; 238 for(int j=0; j<mat.outerSize(); ++j) 239 for(typename SparseMatrixType::InnerIterator it(mat,j); it; ++it) 240 { 241 ++ count; 242 internal::PutMatrixElt(it.value(), it.row()+1, it.col()+1, out); 243 // out << it.row()+1 << " " << it.col()+1 << " " << it.value() << "\n"; 244 } 245 out.close(); 246 return true; 247 } 248 249 template<typename VectorType> 250 bool saveMarketVector (const VectorType& vec, const std::string& filename) 251 { 252 typedef typename VectorType::Scalar Scalar; 253 std::ofstream out(filename.c_str(),std::ios::out); 254 if(!out) 255 return false; 256 257 out.flags(std::ios_base::scientific); 258 out.precision(64); 259 if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value) 260 out << "%%MatrixMarket matrix array complex general\n"; 261 else 262 out << "%%MatrixMarket matrix array real general\n"; 263 out << vec.size() << " "<< 1 << "\n"; 264 for (int i=0; i < vec.size(); i++){ 265 internal::putVectorElt(vec(i), out); 266 } 267 out.close(); 268 return true; 269 } 270 271 } // end namespace Eigen 272 273 #endif // EIGEN_SPARSE_MARKET_IO_H 274
