Home | History | Annotate | Download | only in SparseExtra
      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 (at) inria.fr>
      5 // Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam (at) 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>
     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 
    109 inline 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); 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 
    132 template<typename SparseMatrixType>
    133 bool 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 
    191 template<typename VectorType>
    192 bool 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); assert(in.good());
    204   } while (line[0] == '%');
    205   std::istringstream newline(line);
    206   newline  >> n >> col;
    207   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 
    223 template<typename SparseMatrixType>
    224 bool 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 
    249 template<typename VectorType>
    250 bool 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