Adding support to read and write complex matrices in Matrix Market format

unsupported/Eigen/src/SparseExtra/MarketIO.h

@@ -2,6 +2,7 @@
 // for linear algebra.
 //
 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam@inria.fr>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -24,7 +25,119 @@
 
 #ifndef EIGEN_SPARSE_MARKET_IO_H
 #define EIGEN_SPARSE_MARKET_IO_H
+namespace internal 
+{
+  template <typename Scalar>
+  inline bool GetMarketLine (std::stringstream& line, int& M, int& N, int& i, int& j, Scalar& value)
+  {
+    line >> i >> j >> value;
+    i--;
+    j--;
+    if(i>=0 && j>=0 && i<M && j<N)
+    {
+      return true; 
+    }
+    else
+      return false;
+  }
+  template <typename Scalar>
+  inline bool GetMarketLine (std::stringstream& line, int& M, int& N, int& i, int& j, std::complex<Scalar>& value)
+  {
+    Scalar valR, valI;
+    line >> i >> j >> valR >> valI;
+    i--;
+    j--;
+    if(i>=0 && j>=0 && i<M && j<N)
+    {
+      value = std::complex<Scalar>(valR, valI);
+      return true; 
+    }
+    else
+      return false;
+  }
 
+  template <typename RealScalar>
+  inline void  GetVectorElt (const std::string& line, RealScalar& val)
+  {
+    std::istringstream newline(line);
+    newline >> val;  
+  }
+
+  template <typename RealScalar>
+  inline void GetVectorElt (const std::string& line, std::complex<RealScalar>& val)
+  {
+    RealScalar valR, valI; 
+    std::istringstream newline(line);
+    newline >> valR >> valI; 
+    val = std::complex<RealScalar>(valR, valI);
+  }
+  
+  template<typename Scalar>
+  inline void putMarketHeader(std::string& header,int sym)
+  {
+    header= "%%MatrixMarket matrix coordinate ";
+    if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value)
+    {
+      header += " complex"; 
+      if(sym == Symmetric) header += " symmetric";
+      else if (sym == SelfAdjoint) header += " hermitian";
+      else header += " general";
+    }
+    else
+    {
+      header += " real"; 
+      if(sym == Symmetric) header += " symmetric";
+      else header += " general";
+    }
+  }
+
+  template<typename Scalar>
+  inline void PutMatrixElt(Scalar value, int row, int col, std::ofstream& out)
+  {
+    out << row << " "<< col << " " << value << "\n";
+  }
+  template<typename Scalar>
+  inline void PutMatrixElt(std::complex<Scalar> value, int row, int col, std::ofstream& out)
+  {
+    out << row << " " << col << " " << value.real() << " " << value.imag() << "\n";
+  }
+
+
+  template<typename Scalar>
+  inline void putVectorElt(Scalar value, std::ofstream& out)
+  {
+    out << value << "\n"; 
+  }
+  template<typename Scalar>
+  inline void putVectorElt(std::complex<Scalar> value, std::ofstream& out)
+  {
+    out << value.real << " " << value.imag()<< "\n"; 
+  }
+
+}
+
+inline bool getMarketHeader(const std::string& filename, int& sym, bool& iscomplex, bool& isvector)
+{
+  sym = 0; 
+  isvector = false;
+  std::ifstream in(filename.c_str(),std::ios::in);
+  if(!in)
+    return false;
+  
+  std::string line; 
+  // The matrix header is always the first line in the file 
+  std::getline(in, line); assert(in.good());
+  
+  std::stringstream fmtline(line); 
+  std::string substr[5];
+  fmtline>> substr[0] >> substr[1] >> substr[2] >> substr[3] >> substr[4];
+  if(substr[2].compare("array") == 0) isvector = true;
+  if(substr[3].compare("complex") == 0) iscomplex = true;
+  if(substr[4].compare("symmetric") == 0) sym = Symmetric;
+  else if (substr[4].compare("hermitian") == 0) sym = SelfAdjoint;
+  
+  return true;
+}
   
 template<typename SparseMatrixType>
 bool loadMarket(SparseMatrixType& mat, const std::string& filename)
@@ -41,10 +154,10 @@ bool loadMarket(SparseMatrixType& mat, const std::string& filename)
   
   int M(-1), N(-1), NNZ(-1);
   int count = 0;
-  
   while(input.getline(buffer, maxBuffersize))
   {
     // skip comments   
+    //NOTE An appropriate test should be done on the header to get the  symmetry
     if(buffer[0]=='%')
       continue;
     
@@ -61,15 +174,11 @@ bool loadMarket(SparseMatrixType& mat, const std::string& filename)
     else
     { 
       int i(-1), j(-1);
-      Scalar v;
+      Scalar value; 
-      line >> i >> j >> v;
+      if( internal::GetMarketLine(line, M, N, i, j, value) ) 
-      i--;
-      j--;
-      if(i>=0 && j>=0 && i<M && j<N)
       {
         ++ count;
-        //std::cout << "M[" << i << "," << j << "] = " << v << "\n";
+        mat.insert(i,j) = value;
-        mat.insert(i,j) = v;
       }
       else 
         std::cerr << "Invalid read: " << i << "," << j << "\n";        
@@ -83,25 +192,83 @@ bool loadMarket(SparseMatrixType& mat, const std::string& filename)
   return true;
 }
 
-template<typename SparseMatrixType>
+template<typename VectorType>
-bool saveMarket(const SparseMatrixType& mat, const std::string& filename)
+bool loadMarketVector(VectorType& vec, const std::string& filename)
 {
+   typedef typename VectorType::Scalar Scalar;
+  std::ifstream in(filename.c_str(), std::ios::in);
+  if(!in)
+    return false;
+  
+  std::string line; 
+  int n(0), col(0); 
+  do 
+  { // Skip comments
+    std::getline(in, line); assert(in.good());
+  } while (line[0] == '%');
+  std::istringstream newline(line);
+  newline  >> n >> col; 
+  assert(n>0 && col>0);
+  vec.resize(n);
+  int i = 0; 
+  Scalar value; 
+  while ( std::getline(in, line) && (i < n) ){
+    internal::GetVectorElt(line, value); 
+    vec(i++) = value; 
+  }
+  in.close();
+  if (i!=n){
+    std::cerr<< "Unable to read all elements from file " << filename << "\n";
+    return false;
+  }
+  return true;
+}
+
+template<typename SparseMatrixType>
+bool saveMarket(const SparseMatrixType& mat, const std::string& filename, int sym = 0)
+{
+  typedef typename SparseMatrixType::Scalar Scalar;
   std::ofstream out(filename.c_str(),std::ios::out);
   if(!out)
     return false;
   
   out.flags(std::ios_base::scientific);
   out.precision(64);
+  std::string header; 
+  internal::putMarketHeader<Scalar>(header, sym); 
+  out << header << std::endl; 
   out << mat.rows() << " " << mat.cols() << " " << mat.nonZeros() << "\n";
   int count = 0;
   for(int j=0; j<mat.outerSize(); ++j)
     for(typename SparseMatrixType::InnerIterator it(mat,j); it; ++it)
     {
 	++ count;
-      out << it.row()+1 << " " << it.col()+1 << " " << it.value() << "\n";
+	internal::PutMatrixElt(it.value(), it.row()+1, it.col()+1, out);
+	// out << it.row()+1 << " " << it.col()+1 << " " << it.value() << "\n";
     }
   out.close();
   return true;
 }
 
+template<typename VectorType>
+bool saveMarketVector (const VectorType& vec, const std::string& filename)
+{
+ typedef typename VectorType::Scalar Scalar; 
+ std::ofstream out(filename.c_str(),std::ios::out);
+  if(!out)
+    return false;
+  
+  out.flags(std::ios_base::scientific);
+  out.precision(64);
+  if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value)
+      out << "%%MatrixMarket matrix array complex general\n"; 
+  else
+    out << "%%MatrixMarket matrix array real general\n"; 
+  out << vec.size() << " "<< 1 << "\n";
+  for (int i=0; i < vec.size(); i++){
+    internal::putVectorElt(vec(i), out); 
+  }
+  out.close();
+  return true; 
+}
 #endif // EIGEN_SPARSE_MARKET_IO_H