src/SparseExtra/MarketIO.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 // Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam@inria.fr>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010
0011 #ifndef EIGEN_SPARSE_MARKET_IO_H
0012 #define EIGEN_SPARSE_MARKET_IO_H
0013
0014 #include <iostream>
0015 #include <vector>
0016
0017 namespace Eigen {
0018
0019 namespace internal
0020 {
0021   template <typename Scalar, typename StorageIndex>
0022   inline void GetMarketLine (const char* line, StorageIndex& i, StorageIndex& j, Scalar& value)
0023   {
0024     std::stringstream sline(line);
0025     sline >> i >> j >> value;
0026   }
0027
0028   template<> inline void GetMarketLine (const char* line, int& i, int& j, float& value)
0029   { std::sscanf(line, "%d %d %g", &i, &j, &value); }
0030
0031   template<> inline void GetMarketLine (const char* line, int& i, int& j, double& value)
0032   { std::sscanf(line, "%d %d %lg", &i, &j, &value); }
0033
0034   template<> inline void GetMarketLine (const char* line, int& i, int& j, std::complex<float>& value)
0035   { std::sscanf(line, "%d %d %g %g", &i, &j, &numext::real_ref(value), &numext::imag_ref(value)); }
0036
0037   template<> inline void GetMarketLine (const char* line, int& i, int& j, std::complex<double>& value)
0038   { std::sscanf(line, "%d %d %lg %lg", &i, &j, &numext::real_ref(value), &numext::imag_ref(value)); }
0039
0040   template <typename Scalar, typename StorageIndex>
0041   inline void GetMarketLine (const char* line, StorageIndex& i, StorageIndex& j, std::complex<Scalar>& value)
0042   {
0043     std::stringstream sline(line);
0044     Scalar valR, valI;
0045     sline >> i >> j >> valR >> valI;
0046     value = std::complex<Scalar>(valR,valI);
0047   }
0048
0049   template <typename RealScalar>
0050   inline void  GetVectorElt (const std::string& line, RealScalar& val)
0051   {
0052     std::istringstream newline(line);
0053     newline >> val;
0054   }
0055
0056   template <typename RealScalar>
0057   inline void GetVectorElt (const std::string& line, std::complex<RealScalar>& val)
0058   {
0059     RealScalar valR, valI;
0060     std::istringstream newline(line);
0061     newline >> valR >> valI;
0062     val = std::complex<RealScalar>(valR, valI);
0063   }
0064
0065   template<typename Scalar>
0066   inline void putMarketHeader(std::string& header,int sym)
0067   {
0068     header= "%%MatrixMarket matrix coordinate ";
0069     if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value)
0070     {
0071       header += " complex";
0072       if(sym == Symmetric) header += " symmetric";
0073       else if (sym == SelfAdjoint) header += " Hermitian";
0074       else header += " general";
0075     }
0076     else
0077     {
0078       header += " real";
0079       if(sym == Symmetric) header += " symmetric";
0080       else header += " general";
0081     }
0082   }
0083
0084   template<typename Scalar, typename StorageIndex>
0085   inline void PutMatrixElt(Scalar value, StorageIndex row, StorageIndex col, std::ofstream& out)
0086   {
0087     out << row << " "<< col << " " << value << "\n";
0088   }
0089   template<typename Scalar, typename StorageIndex>
0090   inline void PutMatrixElt(std::complex<Scalar> value, StorageIndex row, StorageIndex col, std::ofstream& out)
0091   {
0092     out << row << " " << col << " " << value.real() << " " << value.imag() << "\n";
0093   }
0094
0095
0096   template<typename Scalar>
0097   inline void putVectorElt(Scalar value, std::ofstream& out)
0098   {
0099     out << value << "\n";
0100   }
0101   template<typename Scalar>
0102   inline void putVectorElt(std::complex<Scalar> value, std::ofstream& out)
0103   {
0104     out << value.real() << " " << value.imag()<< "\n";
0105   }
0106
0107 } // end namespace internal
0108
0109 inline bool getMarketHeader(const std::string& filename, int& sym, bool& iscomplex, bool& isvector)
0110 {
0111   sym = 0;
0112   iscomplex = false;
0113   isvector = false;
0114   std::ifstream in(filename.c_str(),std::ios::in);
0115   if(!in)
0116     return false;
0117
0118   std::string line;
0119   // The matrix header is always the first line in the file
0120   std::getline(in, line); eigen_assert(in.good());
0121
0122   std::stringstream fmtline(line);
0123   std::string substr[5];
0124   fmtline>> substr[0] >> substr[1] >> substr[2] >> substr[3] >> substr[4];
0125   if(substr[2].compare("array") == 0) isvector = true;
0126   if(substr[3].compare("complex") == 0) iscomplex = true;
0127   if(substr[4].compare("symmetric") == 0) sym = Symmetric;
0128   else if (substr[4].compare("Hermitian") == 0) sym = SelfAdjoint;
0129
0130   return true;
0131 }
0132
0133 template<typename SparseMatrixType>
0134 bool loadMarket(SparseMatrixType& mat, const std::string& filename)
0135 {
0136   typedef typename SparseMatrixType::Scalar Scalar;
0137   typedef typename SparseMatrixType::StorageIndex StorageIndex;
0138   std::ifstream input(filename.c_str(),std::ios::in);
0139   if(!input)
0140     return false;
0141
0142   char rdbuffer[4096];
0143   input.rdbuf()->pubsetbuf(rdbuffer, 4096);
0144
0145   const int maxBuffersize = 2048;
0146   char buffer[maxBuffersize];
0147
0148   bool readsizes = false;
0149
0150   typedef Triplet<Scalar,StorageIndex> T;
0151   std::vector<T> elements;
0152
0153   Index M(-1), N(-1), NNZ(-1);
0154   Index count = 0;
0155   while(input.getline(buffer, maxBuffersize))
0156   {
0157     // skip comments
0158     //NOTE An appropriate test should be done on the header to get the  symmetry
0159     if(buffer[0]=='%')
0160       continue;
0161
0162     if(!readsizes)
0163     {
0164       std::stringstream line(buffer);
0165       line >> M >> N >> NNZ;
0166       if(M > 0 && N > 0)
0167       {
0168         readsizes = true;
0169         mat.resize(M,N);
0170         mat.reserve(NNZ);
0171       }
0172     }
0173     else
0174     {
0175       StorageIndex i(-1), j(-1);
0176       Scalar value;
0177       internal::GetMarketLine(buffer, i, j, value);
0178
0179       i--;
0180       j--;
0181       if(i>=0 && j>=0 && i<M && j<N)
0182       {
0183         ++count;
0184         elements.push_back(T(i,j,value));
0185       }
0186       else
0187         std::cerr << "Invalid read: " << i << "," << j << "\n";
0188     }
0189   }
0190
0191   mat.setFromTriplets(elements.begin(), elements.end());
0192   if(count!=NNZ)
0193     std::cerr << count << "!=" << NNZ << "\n";
0194
0195   input.close();
0196   return true;
0197 }
0198
0199 template<typename VectorType>
0200 bool loadMarketVector(VectorType& vec, const std::string& filename)
0201 {
0202    typedef typename VectorType::Scalar Scalar;
0203   std::ifstream in(filename.c_str(), std::ios::in);
0204   if(!in)
0205     return false;
0206
0207   std::string line;
0208   int n(0), col(0);
0209   do
0210   { // Skip comments
0211     std::getline(in, line); eigen_assert(in.good());
0212   } while (line[0] == '%');
0213   std::istringstream newline(line);
0214   newline  >> n >> col;
0215   eigen_assert(n>0 && col>0);
0216   vec.resize(n);
0217   int i = 0;
0218   Scalar value;
0219   while ( std::getline(in, line) && (i < n) ){
0220     internal::GetVectorElt(line, value);
0221     vec(i++) = value;
0222   }
0223   in.close();
0224   if (i!=n){
0225     std::cerr<< "Unable to read all elements from file " << filename << "\n";
0226     return false;
0227   }
0228   return true;
0229 }
0230
0231 template<typename SparseMatrixType>
0232 bool saveMarket(const SparseMatrixType& mat, const std::string& filename, int sym = 0)
0233 {
0234   typedef typename SparseMatrixType::Scalar Scalar;
0235   typedef typename SparseMatrixType::RealScalar RealScalar;
0236   std::ofstream out(filename.c_str(),std::ios::out);
0237   if(!out)
0238     return false;
0239
0240   out.flags(std::ios_base::scientific);
0241   out.precision(std::numeric_limits<RealScalar>::digits10 + 2);
0242   std::string header;
0243   internal::putMarketHeader<Scalar>(header, sym);
0244   out << header << std::endl;
0245   out << mat.rows() << " " << mat.cols() << " " << mat.nonZeros() << "\n";
0246   int count = 0;
0247   for(int j=0; j<mat.outerSize(); ++j)
0248     for(typename SparseMatrixType::InnerIterator it(mat,j); it; ++it)
0249     {
0250       ++ count;
0251       internal::PutMatrixElt(it.value(), it.row()+1, it.col()+1, out);
0252     }
0253   out.close();
0254   return true;
0255 }
0256
0257 template<typename VectorType>
0258 bool saveMarketVector (const VectorType& vec, const std::string& filename)
0259 {
0260  typedef typename VectorType::Scalar Scalar;
0261  typedef typename VectorType::RealScalar RealScalar;
0262  std::ofstream out(filename.c_str(),std::ios::out);
0263   if(!out)
0264     return false;
0265
0266   out.flags(std::ios_base::scientific);
0267   out.precision(std::numeric_limits<RealScalar>::digits10 + 2);
0268   if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value)
0269       out << "%%MatrixMarket matrix array complex general\n";
0270   else
0271     out << "%%MatrixMarket matrix array real general\n";
0272   out << vec.size() << " "<< 1 << "\n";
0273   for (int i=0; i < vec.size(); i++){
0274     internal::putVectorElt(vec(i), out);
0275   }
0276   out.close();
0277   return true;
0278 }
0279
0280 } // end namespace Eigen
0281
0282 #endif // EIGEN_SPARSE_MARKET_IO_H