cart-elc

Source code for CART-ELC
git clone git://git.laack.co/cart-elc.git
Log | Files | Refs | README | LICENSE

MarketIO.h (8416B)

      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 #include <vector>
     16 
     17 namespace Eigen { 
     18 
     19 namespace internal 
     20 {
     21   template <typename Scalar, typename StorageIndex>
     22   inline void GetMarketLine (const char* line, StorageIndex& i, StorageIndex& j, Scalar& value)
     23   {
     24     std::stringstream sline(line);
     25     sline >> i >> j >> value;
     26   }
     27 
     28   template<> inline void GetMarketLine (const char* line, int& i, int& j, float& value)
     29   { std::sscanf(line, "%d %d %g", &i, &j, &value); }
     30 
     31   template<> inline void GetMarketLine (const char* line, int& i, int& j, double& value)
     32   { std::sscanf(line, "%d %d %lg", &i, &j, &value); }
     33 
     34   template<> inline void GetMarketLine (const char* line, int& i, int& j, std::complex<float>& value)
     35   { std::sscanf(line, "%d %d %g %g", &i, &j, &numext::real_ref(value), &numext::imag_ref(value)); }
     36 
     37   template<> inline void GetMarketLine (const char* line, int& i, int& j, std::complex<double>& value)
     38   { std::sscanf(line, "%d %d %lg %lg", &i, &j, &numext::real_ref(value), &numext::imag_ref(value)); }
     39 
     40   template <typename Scalar, typename StorageIndex>
     41   inline void GetMarketLine (const char* line, StorageIndex& i, StorageIndex& j, std::complex<Scalar>& value)
     42   {
     43     std::stringstream sline(line);
     44     Scalar valR, valI;
     45     sline >> i >> j >> valR >> valI;
     46     value = std::complex<Scalar>(valR,valI);
     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, typename StorageIndex>
     85   inline void PutMatrixElt(Scalar value, StorageIndex row, StorageIndex col, std::ofstream& out)
     86   {
     87     out << row << " "<< col << " " << value << "\n";
     88   }
     89   template<typename Scalar, typename StorageIndex>
     90   inline void PutMatrixElt(std::complex<Scalar> value, StorageIndex row, StorageIndex 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 namespace internal
    108 
    109 inline bool getMarketHeader(const std::string& filename, int& sym, bool& iscomplex, bool& isvector)
    110 {
    111   sym = 0; 
    112   iscomplex = false;
    113   isvector = false;
    114   std::ifstream in(filename.c_str(),std::ios::in);
    115   if(!in)
    116     return false;
    117   
    118   std::string line; 
    119   // The matrix header is always the first line in the file 
    120   std::getline(in, line); eigen_assert(in.good());
    121   
    122   std::stringstream fmtline(line); 
    123   std::string substr[5];
    124   fmtline>> substr[0] >> substr[1] >> substr[2] >> substr[3] >> substr[4];
    125   if(substr[2].compare("array") == 0) isvector = true;
    126   if(substr[3].compare("complex") == 0) iscomplex = true;
    127   if(substr[4].compare("symmetric") == 0) sym = Symmetric;
    128   else if (substr[4].compare("Hermitian") == 0) sym = SelfAdjoint;
    129   
    130   return true;
    131 }
    132   
    133 template<typename SparseMatrixType>
    134 bool loadMarket(SparseMatrixType& mat, const std::string& filename)
    135 {
    136   typedef typename SparseMatrixType::Scalar Scalar;
    137   typedef typename SparseMatrixType::StorageIndex StorageIndex;
    138   std::ifstream input(filename.c_str(),std::ios::in);
    139   if(!input)
    140     return false;
    141 
    142   char rdbuffer[4096];
    143   input.rdbuf()->pubsetbuf(rdbuffer, 4096);
    144   
    145   const int maxBuffersize = 2048;
    146   char buffer[maxBuffersize];
    147   
    148   bool readsizes = false;
    149 
    150   typedef Triplet<Scalar,StorageIndex> T;
    151   std::vector<T> elements;
    152   
    153   Index M(-1), N(-1), NNZ(-1);
    154   Index count = 0;
    155   while(input.getline(buffer, maxBuffersize))
    156   {
    157     // skip comments   
    158     //NOTE An appropriate test should be done on the header to get the  symmetry
    159     if(buffer[0]=='%')
    160       continue;
    161 
    162     if(!readsizes)
    163     {
    164       std::stringstream line(buffer);
    165       line >> M >> N >> NNZ;
    166       if(M > 0 && N > 0)
    167       {
    168         readsizes = true;
    169         mat.resize(M,N);
    170         mat.reserve(NNZ);
    171       }
    172     }
    173     else
    174     { 
    175       StorageIndex i(-1), j(-1);
    176       Scalar value; 
    177       internal::GetMarketLine(buffer, i, j, value);
    178 
    179       i--;
    180       j--;
    181       if(i>=0 && j>=0 && i<M && j<N)
    182       {
    183         ++count;
    184         elements.push_back(T(i,j,value));
    185       }
    186       else
    187         std::cerr << "Invalid read: " << i << "," << j << "\n";        
    188     }
    189   }
    190 
    191   mat.setFromTriplets(elements.begin(), elements.end());
    192   if(count!=NNZ)
    193     std::cerr << count << "!=" << NNZ << "\n";
    194   
    195   input.close();
    196   return true;
    197 }
    198 
    199 template<typename VectorType>
    200 bool loadMarketVector(VectorType& vec, const std::string& filename)
    201 {
    202    typedef typename VectorType::Scalar Scalar;
    203   std::ifstream in(filename.c_str(), std::ios::in);
    204   if(!in)
    205     return false;
    206   
    207   std::string line; 
    208   int n(0), col(0); 
    209   do 
    210   { // Skip comments
    211     std::getline(in, line); eigen_assert(in.good());
    212   } while (line[0] == '%');
    213   std::istringstream newline(line);
    214   newline  >> n >> col; 
    215   eigen_assert(n>0 && col>0);
    216   vec.resize(n);
    217   int i = 0; 
    218   Scalar value; 
    219   while ( std::getline(in, line) && (i < n) ){
    220     internal::GetVectorElt(line, value); 
    221     vec(i++) = value; 
    222   }
    223   in.close();
    224   if (i!=n){
    225     std::cerr<< "Unable to read all elements from file " << filename << "\n";
    226     return false;
    227   }
    228   return true;
    229 }
    230 
    231 template<typename SparseMatrixType>
    232 bool saveMarket(const SparseMatrixType& mat, const std::string& filename, int sym = 0)
    233 {
    234   typedef typename SparseMatrixType::Scalar Scalar;
    235   typedef typename SparseMatrixType::RealScalar RealScalar;
    236   std::ofstream out(filename.c_str(),std::ios::out);
    237   if(!out)
    238     return false;
    239   
    240   out.flags(std::ios_base::scientific);
    241   out.precision(std::numeric_limits<RealScalar>::digits10 + 2);
    242   std::string header; 
    243   internal::putMarketHeader<Scalar>(header, sym); 
    244   out << header << std::endl; 
    245   out << mat.rows() << " " << mat.cols() << " " << mat.nonZeros() << "\n";
    246   int count = 0;
    247   for(int j=0; j<mat.outerSize(); ++j)
    248     for(typename SparseMatrixType::InnerIterator it(mat,j); it; ++it)
    249     {
    250       ++ count;
    251       internal::PutMatrixElt(it.value(), it.row()+1, it.col()+1, out);
    252     }
    253   out.close();
    254   return true;
    255 }
    256 
    257 template<typename VectorType>
    258 bool saveMarketVector (const VectorType& vec, const std::string& filename)
    259 {
    260  typedef typename VectorType::Scalar Scalar;
    261  typedef typename VectorType::RealScalar RealScalar;
    262  std::ofstream out(filename.c_str(),std::ios::out);
    263   if(!out)
    264     return false;
    265   
    266   out.flags(std::ios_base::scientific);
    267   out.precision(std::numeric_limits<RealScalar>::digits10 + 2);
    268   if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value)
    269       out << "%%MatrixMarket matrix array complex general\n"; 
    270   else
    271     out << "%%MatrixMarket matrix array real general\n"; 
    272   out << vec.size() << " "<< 1 << "\n";
    273   for (int i=0; i < vec.size(); i++){
    274     internal::putVectorElt(vec(i), out); 
    275   }
    276   out.close();
    277   return true; 
    278 }
    279 
    280 } // end namespace Eigen
    281 
    282 #endif // EIGEN_SPARSE_MARKET_IO_H