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