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, Index& M, Index& N, Index& i, Index& 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, Index& M, Index& N, Index& i, Index& 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); eigen_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 typedef typename SparseMatrixType::Index Index; 137 std::ifstream input(filename.c_str(),std::ios::in); 138 if(!input) 139 return false; 140 141 const int maxBuffersize = 2048; 142 char buffer[maxBuffersize]; 143 144 bool readsizes = false; 145 146 typedef Triplet<Scalar,Index> T; 147 std::vector<T> elements; 148 149 Index M(-1), N(-1), NNZ(-1); 150 Index count = 0; 151 while(input.getline(buffer, maxBuffersize)) 152 { 153 // skip comments 154 //NOTE An appropriate test should be done on the header to get the symmetry 155 if(buffer[0]=='%') 156 continue; 157 158 std::stringstream line(buffer); 159 160 if(!readsizes) 161 { 162 line >> M >> N >> NNZ; 163 if(M > 0 && N > 0 && NNZ > 0) 164 { 165 readsizes = true; 166 //std::cout << "sizes: " << M << "," << N << "," << NNZ << "\n"; 167 mat.resize(M,N); 168 mat.reserve(NNZ); 169 } 170 } 171 else 172 { 173 Index i(-1), j(-1); 174 Scalar value; 175 if( internal::GetMarketLine(line, M, N, i, j, value) ) 176 { 177 ++ count; 178 elements.push_back(T(i,j,value)); 179 } 180 else 181 std::cerr << "Invalid read: " << i << "," << j << "\n"; 182 } 183 } 184 mat.setFromTriplets(elements.begin(), elements.end()); 185 if(count!=NNZ) 186 std::cerr << count << "!=" << NNZ << "\n"; 187 188 input.close(); 189 return true; 190 } 191 192 template<typename VectorType> 193 bool loadMarketVector(VectorType& vec, const std::string& filename) 194 { 195 typedef typename VectorType::Scalar Scalar; 196 std::ifstream in(filename.c_str(), std::ios::in); 197 if(!in) 198 return false; 199 200 std::string line; 201 int n(0), col(0); 202 do 203 { // Skip comments 204 std::getline(in, line); eigen_assert(in.good()); 205 } while (line[0] == '%'); 206 std::istringstream newline(line); 207 newline >> n >> col; 208 eigen_assert(n>0 && col>0); 209 vec.resize(n); 210 int i = 0; 211 Scalar value; 212 while ( std::getline(in, line) && (i < n) ){ 213 internal::GetVectorElt(line, value); 214 vec(i++) = value; 215 } 216 in.close(); 217 if (i!=n){ 218 std::cerr<< "Unable to read all elements from file " << filename << "\n"; 219 return false; 220 } 221 return true; 222 } 223 224 template<typename SparseMatrixType> 225 bool saveMarket(const SparseMatrixType& mat, const std::string& filename, int sym = 0) 226 { 227 typedef typename SparseMatrixType::Scalar Scalar; 228 std::ofstream out(filename.c_str(),std::ios::out); 229 if(!out) 230 return false; 231 232 out.flags(std::ios_base::scientific); 233 out.precision(64); 234 std::string header; 235 internal::putMarketHeader<Scalar>(header, sym); 236 out << header << std::endl; 237 out << mat.rows() << " " << mat.cols() << " " << mat.nonZeros() << "\n"; 238 int count = 0; 239 for(int j=0; j<mat.outerSize(); ++j) 240 for(typename SparseMatrixType::InnerIterator it(mat,j); it; ++it) 241 { 242 ++ count; 243 internal::PutMatrixElt(it.value(), it.row()+1, it.col()+1, out); 244 // out << it.row()+1 << " " << it.col()+1 << " " << it.value() << "\n"; 245 } 246 out.close(); 247 return true; 248 } 249 250 template<typename VectorType> 251 bool saveMarketVector (const VectorType& vec, const std::string& filename) 252 { 253 typedef typename VectorType::Scalar Scalar; 254 std::ofstream out(filename.c_str(),std::ios::out); 255 if(!out) 256 return false; 257 258 out.flags(std::ios_base::scientific); 259 out.precision(64); 260 if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value) 261 out << "%%MatrixMarket matrix array complex general\n"; 262 else 263 out << "%%MatrixMarket matrix array real general\n"; 264 out << vec.size() << " "<< 1 << "\n"; 265 for (int i=0; i < vec.size(); i++){ 266 internal::putVectorElt(vec(i), out); 267 } 268 out.close(); 269 return true; 270 } 271 272 } // end namespace Eigen 273 274 #endif // EIGEN_SPARSE_MARKET_IO_H 275
