MarketIO.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam@inria.fr>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_SPARSE_MARKET_IO_H
00012 #define EIGEN_SPARSE_MARKET_IO_H
00013 
00014 #include <iostream>
00015 
00016 namespace Eigen { 
00017 
00018 namespace internal 
00019 {
00020   template <typename Scalar>
00021   inline bool GetMarketLine (std::stringstream& line, Index& M, Index& N, Index& i, Index& j, Scalar& value)
00022   {
00023     line >> i >> j >> value;
00024     i--;
00025     j--;
00026     if(i>=0 && j>=0 && i<M && j<N)
00027     {
00028       return true; 
00029     }
00030     else
00031       return false;
00032   }
00033   template <typename Scalar>
00034   inline bool GetMarketLine (std::stringstream& line, Index& M, Index& N, Index& i, Index& j, std::complex<Scalar>& value)
00035   {
00036     Scalar valR, valI;
00037     line >> i >> j >> valR >> valI;
00038     i--;
00039     j--;
00040     if(i>=0 && j>=0 && i<M && j<N)
00041     {
00042       value = std::complex<Scalar>(valR, valI);
00043       return true; 
00044     }
00045     else
00046       return false;
00047   }
00048 
00049   template <typename RealScalar>
00050   inline void  GetVectorElt (const std::string& line, RealScalar& val)
00051   {
00052     std::istringstream newline(line);
00053     newline >> val;  
00054   }
00055 
00056   template <typename RealScalar>
00057   inline void GetVectorElt (const std::string& line, std::complex<RealScalar>& val)
00058   {
00059     RealScalar valR, valI; 
00060     std::istringstream newline(line);
00061     newline >> valR >> valI; 
00062     val = std::complex<RealScalar>(valR, valI);
00063   }
00064   
00065   template<typename Scalar>
00066   inline void putMarketHeader(std::string& header,int sym)
00067   {
00068     header= "%%MatrixMarket matrix coordinate ";
00069     if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value)
00070     {
00071       header += " complex"; 
00072       if(sym == Symmetric) header += " symmetric";
00073       else if (sym == SelfAdjoint) header += " Hermitian";
00074       else header += " general";
00075     }
00076     else
00077     {
00078       header += " real"; 
00079       if(sym == Symmetric) header += " symmetric";
00080       else header += " general";
00081     }
00082   }
00083 
00084   template<typename Scalar>
00085   inline void PutMatrixElt(Scalar value, int row, int col, std::ofstream& out)
00086   {
00087     out << row << " "<< col << " " << value << "\n";
00088   }
00089   template<typename Scalar>
00090   inline void PutMatrixElt(std::complex<Scalar> value, int row, int col, std::ofstream& out)
00091   {
00092     out << row << " " << col << " " << value.real() << " " << value.imag() << "\n";
00093   }
00094 
00095 
00096   template<typename Scalar>
00097   inline void putVectorElt(Scalar value, std::ofstream& out)
00098   {
00099     out << value << "\n"; 
00100   }
00101   template<typename Scalar>
00102   inline void putVectorElt(std::complex<Scalar> value, std::ofstream& out)
00103   {
00104     out << value.real << " " << value.imag()<< "\n"; 
00105   }
00106 
00107 } // end namepsace internal
00108 
00109 inline bool getMarketHeader(const std::string& filename, int& sym, bool& iscomplex, bool& isvector)
00110 {
00111   sym = 0; 
00112   isvector = false;
00113   std::ifstream in(filename.c_str(),std::ios::in);
00114   if(!in)
00115     return false;
00116   
00117   std::string line; 
00118   // The matrix header is always the first line in the file 
00119   std::getline(in, line); eigen_assert(in.good());
00120   
00121   std::stringstream fmtline(line); 
00122   std::string substr[5];
00123   fmtline>> substr[0] >> substr[1] >> substr[2] >> substr[3] >> substr[4];
00124   if(substr[2].compare("array") == 0) isvector = true;
00125   if(substr[3].compare("complex") == 0) iscomplex = true;
00126   if(substr[4].compare("symmetric") == 0) sym = Symmetric;
00127   else if (substr[4].compare("Hermitian") == 0) sym = SelfAdjoint;
00128   
00129   return true;
00130 }
00131   
00132 template<typename SparseMatrixType>
00133 bool loadMarket(SparseMatrixType& mat, const std::string& filename)
00134 {
00135   typedef typename SparseMatrixType::Scalar Scalar;
00136   typedef typename SparseMatrixType::Index Index;
00137   std::ifstream input(filename.c_str(),std::ios::in);
00138   if(!input)
00139     return false;
00140   
00141   const int maxBuffersize = 2048;
00142   char buffer[maxBuffersize];
00143   
00144   bool readsizes = false;
00145 
00146   typedef Triplet<Scalar,Index> T;
00147   std::vector<T> elements;
00148   
00149   Index M(-1), N(-1), NNZ(-1);
00150   Index count = 0;
00151   while(input.getline(buffer, maxBuffersize))
00152   {
00153     // skip comments   
00154     //NOTE An appropriate test should be done on the header to get the  symmetry
00155     if(buffer[0]=='%')
00156       continue;
00157     
00158     std::stringstream line(buffer);
00159     
00160     if(!readsizes)
00161     {
00162       line >> M >> N >> NNZ;
00163       if(M > 0 && N > 0 && NNZ > 0) 
00164       {
00165         readsizes = true;
00166         //std::cout << "sizes: " << M << "," << N << "," << NNZ << "\n";
00167         mat.resize(M,N);
00168         mat.reserve(NNZ);
00169       }
00170     }
00171     else
00172     { 
00173       Index i(-1), j(-1);
00174       Scalar value; 
00175       if( internal::GetMarketLine(line, M, N, i, j, value) ) 
00176       {
00177         ++ count;
00178         elements.push_back(T(i,j,value));
00179       }
00180       else 
00181         std::cerr << "Invalid read: " << i << "," << j << "\n";        
00182     }
00183   }
00184   mat.setFromTriplets(elements.begin(), elements.end());
00185   if(count!=NNZ)
00186     std::cerr << count << "!=" << NNZ << "\n";
00187   
00188   input.close();
00189   return true;
00190 }
00191 
00192 template<typename VectorType>
00193 bool loadMarketVector(VectorType& vec, const std::string& filename)
00194 {
00195    typedef typename VectorType::Scalar Scalar;
00196   std::ifstream in(filename.c_str(), std::ios::in);
00197   if(!in)
00198     return false;
00199   
00200   std::string line; 
00201   int n(0), col(0); 
00202   do 
00203   { // Skip comments
00204     std::getline(in, line); eigen_assert(in.good());
00205   } while (line[0] == '%');
00206   std::istringstream newline(line);
00207   newline  >> n >> col; 
00208   eigen_assert(n>0 && col>0);
00209   vec.resize(n);
00210   int i = 0; 
00211   Scalar value; 
00212   while ( std::getline(in, line) && (i < n) ){
00213     internal::GetVectorElt(line, value); 
00214     vec(i++) = value; 
00215   }
00216   in.close();
00217   if (i!=n){
00218     std::cerr<< "Unable to read all elements from file " << filename << "\n";
00219     return false;
00220   }
00221   return true;
00222 }
00223 
00224 template<typename SparseMatrixType>
00225 bool saveMarket(const SparseMatrixType& mat, const std::string& filename, int sym = 0)
00226 {
00227   typedef typename SparseMatrixType::Scalar Scalar;
00228   std::ofstream out(filename.c_str(),std::ios::out);
00229   if(!out)
00230     return false;
00231   
00232   out.flags(std::ios_base::scientific);
00233   out.precision(64);
00234   std::string header; 
00235   internal::putMarketHeader<Scalar>(header, sym); 
00236   out << header << std::endl; 
00237   out << mat.rows() << " " << mat.cols() << " " << mat.nonZeros() << "\n";
00238   int count = 0;
00239   for(int j=0; j<mat.outerSize(); ++j)
00240     for(typename SparseMatrixType::InnerIterator it(mat,j); it; ++it)
00241     {
00242       ++ count;
00243       internal::PutMatrixElt(it.value(), it.row()+1, it.col()+1, out);
00244       // out << it.row()+1 << " " << it.col()+1 << " " << it.value() << "\n";
00245     }
00246   out.close();
00247   return true;
00248 }
00249 
00250 template<typename VectorType>
00251 bool saveMarketVector (const VectorType& vec, const std::string& filename)
00252 {
00253  typedef typename VectorType::Scalar Scalar; 
00254  std::ofstream out(filename.c_str(),std::ios::out);
00255   if(!out)
00256     return false;
00257   
00258   out.flags(std::ios_base::scientific);
00259   out.precision(64);
00260   if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value)
00261       out << "%%MatrixMarket matrix array complex general\n"; 
00262   else
00263     out << "%%MatrixMarket matrix array real general\n"; 
00264   out << vec.size() << " "<< 1 << "\n";
00265   for (int i=0; i < vec.size(); i++){
00266     internal::putVectorElt(vec(i), out); 
00267   }
00268   out.close();
00269   return true; 
00270 }
00271 
00272 } // end namespace Eigen
00273 
00274 #endif // EIGEN_SPARSE_MARKET_IO_H
 All Classes Functions Variables Typedefs Enumerator