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@inria.fr>
5 // Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam@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,typename IndexType>
GetMarketLine(std::stringstream & line,IndexType & M,IndexType & N,IndexType & i,IndexType & j,Scalar & value)21 inline bool GetMarketLine (std::stringstream& line, IndexType& M, IndexType& N, IndexType& i, IndexType& 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,typename IndexType>
GetMarketLine(std::stringstream & line,IndexType & M,IndexType & N,IndexType & i,IndexType & j,std::complex<Scalar> & value)34 inline bool GetMarketLine (std::stringstream& line, IndexType& M, IndexType& N, IndexType& i, IndexType& 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>
GetVectorElt(const std::string & line,RealScalar & val)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>
GetVectorElt(const std::string & line,std::complex<RealScalar> & val)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>
putMarketHeader(std::string & header,int sym)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>
PutMatrixElt(Scalar value,int row,int col,std::ofstream & out)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>
PutMatrixElt(std::complex<Scalar> value,int row,int col,std::ofstream & out)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>
putVectorElt(Scalar value,std::ofstream & out)97 inline void putVectorElt(Scalar value, std::ofstream& out)
98 {
99 out << value << "\n";
100 }
101 template<typename Scalar>
putVectorElt(std::complex<Scalar> value,std::ofstream & out)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
getMarketHeader(const std::string & filename,int & sym,bool & iscomplex,bool & isvector)109 inline bool getMarketHeader(const std::string& filename, int& sym, bool& iscomplex, bool& isvector)
110 {
111 sym = 0;
112 iscomplex = false;
113 isvector = false;
114 std::ifstream in(filename.c_str(),std::ios::in);
115 if(!in)
116 return false;
117
118 std::string line;
119 // The matrix header is always the first line in the file
120 std::getline(in, line); eigen_assert(in.good());
121
122 std::stringstream fmtline(line);
123 std::string substr[5];
124 fmtline>> substr[0] >> substr[1] >> substr[2] >> substr[3] >> substr[4];
125 if(substr[2].compare("array") == 0) isvector = true;
126 if(substr[3].compare("complex") == 0) iscomplex = true;
127 if(substr[4].compare("symmetric") == 0) sym = Symmetric;
128 else if (substr[4].compare("Hermitian") == 0) sym = SelfAdjoint;
129
130 return true;
131 }
132
133 template<typename SparseMatrixType>
loadMarket(SparseMatrixType & mat,const std::string & filename)134 bool loadMarket(SparseMatrixType& mat, const std::string& filename)
135 {
136 typedef typename SparseMatrixType::Scalar Scalar;
137 typedef typename SparseMatrixType::StorageIndex StorageIndex;
138 std::ifstream input(filename.c_str(),std::ios::in);
139 if(!input)
140 return false;
141
142 const int maxBuffersize = 2048;
143 char buffer[maxBuffersize];
144
145 bool readsizes = false;
146
147 typedef Triplet<Scalar,StorageIndex> T;
148 std::vector<T> elements;
149
150 StorageIndex M(-1), N(-1), NNZ(-1);
151 StorageIndex count = 0;
152 while(input.getline(buffer, maxBuffersize))
153 {
154 // skip comments
155 //NOTE An appropriate test should be done on the header to get the symmetry
156 if(buffer[0]=='%')
157 continue;
158
159 std::stringstream line(buffer);
160
161 if(!readsizes)
162 {
163 line >> M >> N >> NNZ;
164 if(M > 0 && N > 0 && NNZ > 0)
165 {
166 readsizes = true;
167 //std::cout << "sizes: " << M << "," << N << "," << NNZ << "\n";
168 mat.resize(M,N);
169 mat.reserve(NNZ);
170 }
171 }
172 else
173 {
174 StorageIndex i(-1), j(-1);
175 Scalar value;
176 if( internal::GetMarketLine(line, M, N, i, j, value) )
177 {
178 ++ count;
179 elements.push_back(T(i,j,value));
180 }
181 else
182 std::cerr << "Invalid read: " << i << "," << j << "\n";
183 }
184 }
185 mat.setFromTriplets(elements.begin(), elements.end());
186 if(count!=NNZ)
187 std::cerr << count << "!=" << NNZ << "\n";
188
189 input.close();
190 return true;
191 }
192
193 template<typename VectorType>
loadMarketVector(VectorType & vec,const std::string & filename)194 bool loadMarketVector(VectorType& vec, const std::string& filename)
195 {
196 typedef typename VectorType::Scalar Scalar;
197 std::ifstream in(filename.c_str(), std::ios::in);
198 if(!in)
199 return false;
200
201 std::string line;
202 int n(0), col(0);
203 do
204 { // Skip comments
205 std::getline(in, line); eigen_assert(in.good());
206 } while (line[0] == '%');
207 std::istringstream newline(line);
208 newline >> n >> col;
209 eigen_assert(n>0 && col>0);
210 vec.resize(n);
211 int i = 0;
212 Scalar value;
213 while ( std::getline(in, line) && (i < n) ){
214 internal::GetVectorElt(line, value);
215 vec(i++) = value;
216 }
217 in.close();
218 if (i!=n){
219 std::cerr<< "Unable to read all elements from file " << filename << "\n";
220 return false;
221 }
222 return true;
223 }
224
225 template<typename SparseMatrixType>
226 bool saveMarket(const SparseMatrixType& mat, const std::string& filename, int sym = 0)
227 {
228 typedef typename SparseMatrixType::Scalar Scalar;
229 std::ofstream out(filename.c_str(),std::ios::out);
230 if(!out)
231 return false;
232
233 out.flags(std::ios_base::scientific);
234 out.precision(64);
235 std::string header;
236 internal::putMarketHeader<Scalar>(header, sym);
237 out << header << std::endl;
238 out << mat.rows() << " " << mat.cols() << " " << mat.nonZeros() << "\n";
239 int count = 0;
240 for(int j=0; j<mat.outerSize(); ++j)
241 for(typename SparseMatrixType::InnerIterator it(mat,j); it; ++it)
242 {
243 ++ count;
244 internal::PutMatrixElt(it.value(), it.row()+1, it.col()+1, out);
245 // out << it.row()+1 << " " << it.col()+1 << " " << it.value() << "\n";
246 }
247 out.close();
248 return true;
249 }
250
251 template<typename VectorType>
saveMarketVector(const VectorType & vec,const std::string & filename)252 bool saveMarketVector (const VectorType& vec, const std::string& filename)
253 {
254 typedef typename VectorType::Scalar Scalar;
255 std::ofstream out(filename.c_str(),std::ios::out);
256 if(!out)
257 return false;
258
259 out.flags(std::ios_base::scientific);
260 out.precision(64);
261 if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value)
262 out << "%%MatrixMarket matrix array complex general\n";
263 else
264 out << "%%MatrixMarket matrix array real general\n";
265 out << vec.size() << " "<< 1 << "\n";
266 for (int i=0; i < vec.size(); i++){
267 internal::putVectorElt(vec(i), out);
268 }
269 out.close();
270 return true;
271 }
272
273 } // end namespace Eigen
274
275 #endif // EIGEN_SPARSE_MARKET_IO_H
276