Added routines to read and write parallel matrices in a matrix market format.

[[Imported from SVN: r1478]]

dune/istl/matrixmarket.hh

@@ -4,7 +4,12 @@
 #define DUNE_MATRIXMARKET_HH
 
 #include <ostream>
+#include <istream>
+#include <fstream>
+#include <sstream>
 #include <limits>
+#include <ios>
+#include "matrixutils.hh"
 #include "bcrsmatrix.hh"
 #include <dune/common/fmatrix.hh>
 #include <dune/common/tuples.hh>
@@ -155,11 +160,17 @@ namespace Dune
 
     bool lineFeed(std::istream& file)
     {
-      char c=file.peek();
+      char c;
+      if(!file.eof())
+        c=file.peek();
+      else
+        return false;
       // ignore whitespace
       while(c==' ')
       {
         file.get();
+        if(file.eof())
+          return false;
         c=file.peek();
       }
 
@@ -189,13 +200,13 @@ namespace Dune
     {
       std::string buffer;
       char c;
-      c=file.peek();
+      file >> buffer;
+      c=buffer[0];
       mmHeader=MMHeader();
       if(c!='%')
         return false;
-
+      std::cout<<buffer<<std::endl;
       /* read the banner */
-      file >> buffer;
       if(buffer!="%%MatrixMarket") {
         /* disgard the rest of the line */
         file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
@@ -652,12 +663,152 @@ namespace Dune
 
     readSparseEntries(matrix, istr, entries, header,d);
   }
+
+  template<typename M>
+  struct mm_multipliers
+  {};
+
+  template<typename B, int i, int j, typename A>
+  struct mm_multipliers<BCRSMatrix<FieldMatrix<B,i,j>,A> >
+  {
+    enum {
+      rows = i,
+      cols = j
+    };
+  };
+
+  template<typename B, int i, int j>
+  void mm_print_entry(const FieldMatrix<B,i,j>& entry,
+                      typename FieldMatrix<B,i,j>::size_type rowidx,
+                      typename FieldMatrix<B,i,j>::size_type colidx,
+                      std::ostream& ostr)
+  {
+    typedef typename FieldMatrix<B,i,j>::const_iterator riterator;
+    typedef typename FieldMatrix<B,i,j>::ConstColIterator citerator;
+    for(riterator row=entry.begin(); row != entry.end(); ++row, ++rowidx) {
+      int coli=colidx;
+      for(citerator col = row->begin(); col != row->end(); ++col, ++coli)
+        ostr<< rowidx<<" "<<coli<<" "<<*col<<std::endl;
+    }
+  }
+
+
+  /**
+   * @brief writes a ISTL block matrix to a stream in matrix market format.
+   */
   template<typename M>
-  void printMatrixMarket(M& matrix,
-                         std::ostream ostr)
+  void writeMatrixMarket(const M& matrix,
+                         std::ostream& ostr)
   {
-    mm_header_printer<M>::print(ostr, matrix);
+    mm_header_printer<M>::print(ostr);
     mm_block_structure_header<M>::print(ostr,matrix);
+    ostr<<matrix.M()*mm_multipliers<M>::rows<<" "
+        <<matrix.N()*mm_multipliers<M>::cols<<" "
+        <<countNonZeros(matrix)<<std::endl;
+
+    typedef typename M::const_iterator riterator;
+    typedef typename M::ConstColIterator citerator;
+    for(riterator row=matrix.begin(); row != matrix.end(); ++row)
+      for(citerator col = row->begin(); col != row->end(); ++col)
+        // Matrix Market indexing start with 1!
+        mm_print_entry(*col, row.index()*mm_multipliers<M>::rows+1,
+                       col.index()*mm_multipliers<M>::cols+1, ostr);
+  }
+
+  template<typename M, typename G, typename L>
+  void storeMatrixMarket(const M& matrix,
+                         std::string filename,
+                         const OwnerOverlapCopyCommunication<G,L>& comm)
+  {
+    // Get our rank
+    int rank = comm.communicator().rank();
+    // Write the local matrix
+    std::ostringstream rfilename;
+    rfilename<<filename <<"_"<<rank<<".mm";
+    std::cout<< rfilename.str()<<std::endl;
+    std::ofstream file(rfilename.str().c_str());
+    file.setf(std::ios::scientific,std::ios::floatfield);
+    writeMatrixMarket(matrix, file);
+    file.close();
+
+    // Write the global to local index mapping
+    rfilename.str("");
+    rfilename<<filename<<"_"<<rank<<".idx";
+    file.open(rfilename.str().c_str());
+    file.setf(std::ios::scientific,std::ios::floatfield);
+    typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet IndexSet;
+    typedef typename IndexSet::const_iterator Iterator;
+    for(Iterator iter = comm.indexSet().begin();
+        iter != comm.indexSet().end(); ++iter) {
+      file << iter->global()<<" "<<(std::size_t)iter->local()<<" "
+           <<iter->local().attribute()<<" "<<(int)iter->local().isPublic()<<std::endl;
+    }
+    // Store neighbour information for efficient remote indices setup.
+    file<<"neighbours:";
+    const std::set<int>& neighbours=comm.remoteIndices().getNeighbours();
+    typedef std::set<int>::const_iterator SIter;
+    for(SIter neighbour=neighbours.begin(); neighbour != neighbours.end(); ++neighbour) {
+      file<<" "<< *neighbour;
+    }
+    file.close();
+  }
+
+  template<typename M, typename G, typename L>
+  void loadMatrixMarket(M& matrix,
+                        const std::string& filename,
+                        OwnerOverlapCopyCommunication<G,L>& comm)
+  {
+    typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet::LocalIndex LocalIndex;
+    typedef typename LocalIndex::Attribute Attribute;
+    // Get our rank
+    int rank = comm.communicator().rank();
+    // load local matrix
+    std::ostringstream rfilename;
+    rfilename<<filename <<"_"<<rank<<".mm";
+    std::ifstream file;
+    file.open(rfilename.str().c_str(), std::ios::in);
+    if(!file)
+      DUNE_THROW(IOError, "Could not open file" << rfilename.str().c_str());
+    //if(!file.is_open())
+    //
+    readMatrixMarket(matrix,file);
+    file.close();
+    // read indices
+    typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet IndexSet;
+    IndexSet& pis=comm.pis;
+    rfilename.str("");
+    rfilename<<filename<<"_"<<rank<<".idx";
+    file.open(rfilename.str().c_str());
+    pis.beginResize();
+    while(!file.eof() && file.peek()!='n') {
+      G g;
+      file >>g;
+      std::size_t l;
+      file >>l;
+      char c;
+      file >>c;
+      bool b;
+      file >> b;
+      pis.add(g,LocalIndex(l,Attribute(c),b));
+      lineFeed(file);
+    }
+    pis.endResize();
+    if(!file.eof()) {
+      // read neighbours
+      std::string s;
+      file>>s;
+      if(s!="neighbours:")
+        DUNE_THROW(MatrixMarketFormatError, "was exspecting the string: \"neighbours:\"");
+      std::set<int> nb;
+      while(!file.eof()) {
+        int i;
+        file >> i;
+        nb.insert(i);
+      }
+      file.close();
+      comm.ri.setNeighbours(nb);
+    }
+    comm.ri.template rebuild<false>();
   }
 
 }

dune/istl/test/matrixmarkettest.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#include "config.h"
+#include "mpi.h"
+#include <dune/istl/io.hh>
+#include <dune/istl/bvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/istl/paamg/test/anisotropic.hh>
+#include <dune/common/timer.hh>
+#include <dune/istl/matrixmarket.hh>
+
+#include <iterator>
+
+int main(int argc, char** argv)
+{
+
+  MPI_Init(&argc, &argv);
+  int size;
+  MPI_Comm_size(MPI_COMM_WORLD, &size);
+  const int BS=1;
+  int N=100;
+
+  if(argc>1)
+    N = atoi(argv[1]);
+  std::cout<<"testing for N="<<N<<" BS="<<1<<std::endl;
+
+
+  typedef Dune::FieldMatrix<double,BS,BS> MatrixBlock;
+  typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
+  typedef Dune::FieldVector<double,BS> VectorBlock;
+  typedef Dune::BlockVector<VectorBlock> BVector;
+  typedef int GlobalId;
+  typedef Dune::OwnerOverlapCopyCommunication<GlobalId> Communication;
+  Communication comm(MPI_COMM_WORLD);
+  int n;
+
+  std::cout<<comm.communicator().rank()<<" "<<comm.communicator().size()<<
+  " "<<size<<std::endl;
+
+  BCRSMat mat = setupAnisotropic2d<BS,double>(N, comm.indexSet(), comm.communicator(), &n, .011);
+
+  storeMatrixMarket(mat, std::string("testmat"), comm);
+
+  BCRSMat mat1;
+  Communication comm1(MPI_COMM_WORLD);
+
+  loadMatrixMarket(mat1, std::string("testmat"), comm1);
+
+  int ret=0;
+  // if(mat!=mat1)
+  //   {
+  //     std::cerr<<"written and read matrix do not match"<<std::endl;
+  //     ++ret;
+  //   }
+  if(comm1.indexSet()!=comm.indexSet())
+  {
+    std::cerr<<"written and read idxset do not match"<<std::endl;
+    ++ret;
+  }
+  storeMatrixMarket(mat1, std::string("testmat1"), comm1);
+  if(ret!=0)
+    MPI_Abort(MPI_COMM_WORLD, ret);
+  MPI_Finalize();
+}