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 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00025 
00026 #ifndef EIGEN_SPARSE_MARKET_IO_H
00027 #define EIGEN_SPARSE_MARKET_IO_H
00028 
00029 #include <iostream>
00030 
00031 namespace Eigen { 
00032 
00033 namespace internal 
00034 {
00035   template <typename Scalar>
00036   inline bool GetMarketLine (std::stringstream& line, int& M, int& N, int& i, int& j, Scalar& value)
00037   {
00038     line >> i >> j >> value;
00039     i--;
00040     j--;
00041     if(i>=0 && j>=0 && i<M && j<N)
00042     {
00043       return true; 
00044     }
00045     else
00046       return false;
00047   }
00048   template <typename Scalar>
00049   inline bool GetMarketLine (std::stringstream& line, int& M, int& N, int& i, int& j, std::complex<Scalar>& value)
00050   {
00051     Scalar valR, valI;
00052     line >> i >> j >> valR >> valI;
00053     i--;
00054     j--;
00055     if(i>=0 && j>=0 && i<M && j<N)
00056     {
00057       value = std::complex<Scalar>(valR, valI);
00058       return true; 
00059     }
00060     else
00061       return false;
00062   }
00063 
00064   template <typename RealScalar>
00065   inline void  GetVectorElt (const std::string& line, RealScalar& val)
00066   {
00067     std::istringstream newline(line);
00068     newline >> val;  
00069   }
00070 
00071   template <typename RealScalar>
00072   inline void GetVectorElt (const std::string& line, std::complex<RealScalar>& val)
00073   {
00074     RealScalar valR, valI; 
00075     std::istringstream newline(line);
00076     newline >> valR >> valI; 
00077     val = std::complex<RealScalar>(valR, valI);
00078   }
00079   
00080   template<typename Scalar>
00081   inline void putMarketHeader(std::string& header,int sym)
00082   {
00083     header= "%%MatrixMarket matrix coordinate ";
00084     if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value)
00085     {
00086       header += " complex"; 
00087       if(sym == Symmetric) header += " symmetric";
00088       else if (sym == SelfAdjoint) header += " Hermitian";
00089       else header += " general";
00090     }
00091     else
00092     {
00093       header += " real"; 
00094       if(sym == Symmetric) header += " symmetric";
00095       else header += " general";
00096     }
00097   }
00098 
00099   template<typename Scalar>
00100   inline void PutMatrixElt(Scalar value, int row, int col, std::ofstream& out)
00101   {
00102     out << row << " "<< col << " " << value << "\n";
00103   }
00104   template<typename Scalar>
00105   inline void PutMatrixElt(std::complex<Scalar> value, int row, int col, std::ofstream& out)
00106   {
00107     out << row << " " << col << " " << value.real() << " " << value.imag() << "\n";
00108   }
00109 
00110 
00111   template<typename Scalar>
00112   inline void putVectorElt(Scalar value, std::ofstream& out)
00113   {
00114     out << value << "\n"; 
00115   }
00116   template<typename Scalar>
00117   inline void putVectorElt(std::complex<Scalar> value, std::ofstream& out)
00118   {
00119     out << value.real << " " << value.imag()<< "\n"; 
00120   }
00121 
00122 } // end namepsace internal
00123 
00124 inline bool getMarketHeader(const std::string& filename, int& sym, bool& iscomplex, bool& isvector)
00125 {
00126   sym = 0; 
00127   isvector = false;
00128   std::ifstream in(filename.c_str(),std::ios::in);
00129   if(!in)
00130     return false;
00131   
00132   std::string line; 
00133   // The matrix header is always the first line in the file 
00134   std::getline(in, line); assert(in.good());
00135   
00136   std::stringstream fmtline(line); 
00137   std::string substr[5];
00138   fmtline>> substr[0] >> substr[1] >> substr[2] >> substr[3] >> substr[4];
00139   if(substr[2].compare("array") == 0) isvector = true;
00140   if(substr[3].compare("complex") == 0) iscomplex = true;
00141   if(substr[4].compare("symmetric") == 0) sym = Symmetric;
00142   else if (substr[4].compare("Hermitian") == 0) sym = SelfAdjoint;
00143   
00144   return true;
00145 }
00146   
00147 template<typename SparseMatrixType>
00148 bool loadMarket(SparseMatrixType& mat, const std::string& filename)
00149 {
00150   typedef typename SparseMatrixType::Scalar Scalar;
00151   std::ifstream input(filename.c_str(),std::ios::in);
00152   if(!input)
00153     return false;
00154   
00155   const int maxBuffersize = 2048;
00156   char buffer[maxBuffersize];
00157   
00158   bool readsizes = false;
00159   
00160   int M(-1), N(-1), NNZ(-1);
00161   int count = 0;
00162   while(input.getline(buffer, maxBuffersize))
00163   {
00164     // skip comments   
00165     //NOTE An appropriate test should be done on the header to get the  symmetry
00166     if(buffer[0]=='%')
00167       continue;
00168     
00169     std::stringstream line(buffer);
00170     
00171     if(!readsizes)
00172     {
00173       line >> M >> N >> NNZ;
00174       if(M > 0 && N > 0 && NNZ > 0) 
00175       {
00176         readsizes = true;
00177         std::cout << "sizes: " << M << "," << N << "," << NNZ << "\n";
00178         mat.resize(M,N);
00179         mat.reserve(NNZ);
00180       }
00181     }
00182     else
00183     { 
00184       int i(-1), j(-1);
00185       Scalar value; 
00186       if( internal::GetMarketLine(line, M, N, i, j, value) ) 
00187       {
00188         ++ count;
00189         mat.insert(i,j) = value;
00190       }
00191       else 
00192         std::cerr << "Invalid read: " << i << "," << j << "\n";        
00193     }
00194   }
00195   mat.makeCompressed();
00196   if(count!=NNZ)
00197     std::cerr << count << "!=" << NNZ << "\n";
00198   
00199   input.close();
00200   return true;
00201 }
00202 
00203 template<typename VectorType>
00204 bool loadMarketVector(VectorType& vec, const std::string& filename)
00205 {
00206    typedef typename VectorType::Scalar Scalar;
00207   std::ifstream in(filename.c_str(), std::ios::in);
00208   if(!in)
00209     return false;
00210   
00211   std::string line; 
00212   int n(0), col(0); 
00213   do 
00214   { // Skip comments
00215     std::getline(in, line); assert(in.good());
00216   } while (line[0] == '%');
00217   std::istringstream newline(line);
00218   newline  >> n >> col; 
00219   assert(n>0 && col>0);
00220   vec.resize(n);
00221   int i = 0; 
00222   Scalar value; 
00223   while ( std::getline(in, line) && (i < n) ){
00224     internal::GetVectorElt(line, value); 
00225     vec(i++) = value; 
00226   }
00227   in.close();
00228   if (i!=n){
00229     std::cerr<< "Unable to read all elements from file " << filename << "\n";
00230     return false;
00231   }
00232   return true;
00233 }
00234 
00235 template<typename SparseMatrixType>
00236 bool saveMarket(const SparseMatrixType& mat, const std::string& filename, int sym = 0)
00237 {
00238   typedef typename SparseMatrixType::Scalar Scalar;
00239   std::ofstream out(filename.c_str(),std::ios::out);
00240   if(!out)
00241     return false;
00242   
00243   out.flags(std::ios_base::scientific);
00244   out.precision(64);
00245   std::string header; 
00246   internal::putMarketHeader<Scalar>(header, sym); 
00247   out << header << std::endl; 
00248   out << mat.rows() << " " << mat.cols() << " " << mat.nonZeros() << "\n";
00249   int count = 0;
00250   for(int j=0; j<mat.outerSize(); ++j)
00251     for(typename SparseMatrixType::InnerIterator it(mat,j); it; ++it)
00252     {
00253         ++ count;
00254         internal::PutMatrixElt(it.value(), it.row()+1, it.col()+1, out);
00255         // out << it.row()+1 << " " << it.col()+1 << " " << it.value() << "\n";
00256     }
00257   out.close();
00258   return true;
00259 }
00260 
00261 template<typename VectorType>
00262 bool saveMarketVector (const VectorType& vec, const std::string& filename)
00263 {
00264  typedef typename VectorType::Scalar Scalar; 
00265  std::ofstream out(filename.c_str(),std::ios::out);
00266   if(!out)
00267     return false;
00268   
00269   out.flags(std::ios_base::scientific);
00270   out.precision(64);
00271   if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value)
00272       out << "%%MatrixMarket matrix array complex general\n"; 
00273   else
00274     out << "%%MatrixMarket matrix array real general\n"; 
00275   out << vec.size() << " "<< 1 << "\n";
00276   for (int i=0; i < vec.size(); i++){
00277     internal::putVectorElt(vec(i), out); 
00278   }
00279   out.close();
00280   return true; 
00281 }
00282 
00283 } // end namespace Eigen
00284 
00285 #endif // EIGEN_SPARSE_MARKET_IO_H