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
16namespace Eigen {
17
18namespace 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
109inline 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
132template<typename SparseMatrixType>
133bool 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
191template<typename VectorType>
192bool 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); eigen_assert(in.good());
204  } while (line[0] == '%');
205  std::istringstream newline(line);
206  newline  >> n >> col;
207  eigen_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
223template<typename SparseMatrixType>
224bool 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
249template<typename VectorType>
250bool 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