dune-istl  2.2.1
matrixmarket.hh
Go to the documentation of this file.
1 #ifndef DUNE_MATRIXMARKET_HH
2 #define DUNE_MATRIXMARKET_HH
3 
4 #include<ostream>
5 #include<istream>
6 #include<fstream>
7 #include<sstream>
8 #include<limits>
9 #include<ios>
10 #include"matrixutils.hh"
11 #include "bcrsmatrix.hh"
12 #include"owneroverlapcopy.hh"
13 #include<dune/common/fmatrix.hh>
14 #include<dune/common/tuples.hh>
15 #include<dune/common/misc.hh>
16 
17 namespace Dune
18 {
19 
45  namespace
46  {
56  template<class T>
57  struct mm_numeric_type{
58  enum{
62  is_numeric=false
63  };
64  };
65 
66  template<>
67  struct mm_numeric_type<int>
68  {
69  enum{
73  is_numeric=true
74  };
75 
76  static std::string str()
77  {
78  return "integer";
79  }
80  };
81 
82  template<>
83  struct mm_numeric_type<double>
84  {
85  enum{
89  is_numeric=true
90  };
91 
92  static std::string str()
93  {
94  return "real";
95  }
96  };
97 
98  template<>
99  struct mm_numeric_type<float>
100  {
101  enum{
105  is_numeric=true
106  };
107 
108  static std::string str()
109  {
110  return "real";
111  }
112  };
113 
114  template<>
115  struct mm_numeric_type<std::complex<double> >
116  {
117  enum{
121  is_numeric=true
122  };
123 
124  static std::string str()
125  {
126  return "complex";
127  }
128  };
129 
130  template<>
131  struct mm_numeric_type<std::complex<float> >
132  {
133  enum{
137  is_numeric=true
138  };
139 
140  static std::string str()
141  {
142  return "complex";
143  }
144  };
145 
154  template<class M>
155  struct mm_header_printer;
156 
157  template<typename T, typename A, int i, int j>
158  struct mm_header_printer<BCRSMatrix<FieldMatrix<T,i,j>,A> >
159  {
160  static void print(std::ostream& os)
161  {
162  os<<"%%MatrixMarket matrix coordinate ";
163  os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
164  }
165  };
166 
167  template<typename B, typename A>
168  struct mm_header_printer<BlockVector<B,A> >
169  {
170  static void print(std::ostream& os)
171  {
172  os<<"%%MatrixMarket matrix array ";
173  os<<mm_numeric_type<typename B::field_type>::str()<<" general"<<std::endl;
174  }
175  };
176 
177  template<typename T, int j>
178  struct mm_header_printer<FieldVector<T,j> >
179  {
180  static void print(std::ostream& os)
181  {
182  os<<"%%MatrixMarket matrix array ";
183  os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
184  }
185  };
186 
187  template<typename T, int i, int j>
188  struct mm_header_printer<FieldMatrix<T,i,j> >
189  {
190  static void print(std::ostream& os)
191  {
192  os<<"%%MatrixMarket matrix array ";
193  os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
194  }
195  };
196 
205  template<class M>
206  struct mm_block_structure_header;
207 
208  template<typename T, typename A, int i>
209  struct mm_block_structure_header<BlockVector<FieldVector<T,i>,A> >
210  {
211  typedef BlockVector<FieldVector<T,i>,A> M;
212 
213  static void print(std::ostream& os, const M& m)
214  {
215  os<<"% ISTL_STRUCT blocked ";
216  os<<i<<" "<<1<<std::endl;
217  }
218  };
219 
220  template<typename T, typename A, int i, int j>
221  struct mm_block_structure_header<BCRSMatrix<FieldMatrix<T,i,j>,A> >
222  {
223  typedef BCRSMatrix<FieldMatrix<T,i,j>,A> M;
224 
225  static void print(std::ostream& os, const M& m)
226  {
227  os<<"% ISTL_STRUCT blocked ";
228  os<<i<<" "<<j<<std::endl;
229  }
230  };
231 
232 
233  template<typename T, int i, int j>
234  struct mm_block_structure_header<FieldMatrix<T,i,j> >
235  {
236  typedef FieldMatrix<T,i,j> M;
237 
238  static void print(std::ostream& os, const M& m)
239  {}
240  };
241 
242  template<typename T, int i>
243  struct mm_block_structure_header<FieldVector<T,i> >
244  {
245  typedef FieldVector<T,i> M;
246 
247  static void print(std::ostream& os, const M& m)
248  {}
249  };
250 
251  enum LineType{ MM_HEADER, MM_ISTLSTRUCT, DATA };
252  enum{ MM_MAX_LINE_LENGTH=1025 };
253 
254  enum MM_TYPE { coordinate_type, array_type, unknown_type };
255 
256  enum MM_CTYPE { integer_type, double_type, complex_type, pattern, unknown_ctype };
257 
258  enum MM_STRUCTURE { general, symmetric, skew_symmetric, hermitian, unknown_structure };
259 
260  struct MMHeader
261  {
262  MMHeader()
263  : type(coordinate_type), ctype(double_type), structure(general)
264  {}
265  MM_TYPE type;
266  MM_CTYPE ctype;
267  MM_STRUCTURE structure;
268  };
269 
270  bool lineFeed(std::istream& file)
271  {
272  char c;
273  if(!file.eof())
274  c=file.peek();
275  else
276  return false;
277  // ignore whitespace
278  while(c==' ')
279  {
280  file.get();
281  if(file.eof())
282  return false;
283  c=file.peek();
284  }
285 
286  if(c=='\n'){
287  /* eat the line feed */
288  file.get();
289  return true;
290  }
291  return false;
292  }
293 
294  void skipComments(std::istream& file)
295  {
296  lineFeed(file);
297  char c=file.peek();
298  // ignore comment lines
299  while(c=='%')
300  {
301  /* disgard the rest of the line */
302  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
303  c=file.peek();
304  }
305  }
306 
307 
308  bool readMatrixMarketBanner(std::istream& file, MMHeader& mmHeader)
309  {
310  std::string buffer;
311  char c;
312  file >> buffer;
313  c=buffer[0];
314  mmHeader=MMHeader();
315  if(c!='%')
316  return false;
317  std::cout<<buffer<<std::endl;
318  /* read the banner */
319  if(buffer!="%%MatrixMarket"){
320  /* disgard the rest of the line */
321  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
322  return false;
323  }
324 
325  if(lineFeed(file))
326  /* premature end of line */
327  return false;
328 
329  /* read the matrix_type */
330  file >> buffer;
331 
332  if(buffer != "matrix")
333  {
334  /* disgard the rest of the line */
335  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
336  return false;
337  }
338 
339  if(lineFeed(file))
340  /* premature end of line */
341  return false;
342 
343  /* The type of the matrix */
344  file >> buffer;
345 
346  if(buffer.empty())
347  return false;
348 
349  std::transform(buffer.begin(), buffer.end(), buffer.begin(),
350  tolower);
351 
352  switch(buffer[0])
353  {
354  case 'a':
355  /* sanity check */
356  if(buffer != "array")
357  {
358  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
359  return false;
360  }
361  mmHeader.type=array_type;
362  break;
363  case 'c':
364  /* sanity check */
365  if(buffer != "coordinate")
366  {
367  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
368  return false;
369  }
370  mmHeader.type=coordinate_type;
371  break;
372  default:
373  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
374  return false;
375  }
376 
377  if(lineFeed(file))
378  /* premature end of line */
379  return false;
380 
381  /* The numeric type used. */
382  file >> buffer;
383 
384  if(buffer.empty())
385  return false;
386 
387  std::transform(buffer.begin(), buffer.end(), buffer.begin(),
388  tolower);
389  switch(buffer[0])
390  {
391  case 'i':
392  /* sanity check */
393  if(buffer != "integer")
394  {
395  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
396  return false;
397  }
398  mmHeader.ctype=integer_type;
399  break;
400  case 'r':
401  /* sanity check */
402  if(buffer != "real")
403  {
404  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
405  return false;
406  }
407  mmHeader.ctype=double_type;
408  break;
409  case 'c':
410  /* sanity check */
411  if(buffer != "complex")
412  {
413  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
414  return false;
415  }
416  mmHeader.ctype=complex_type;
417  break;
418  case 'p':
419  /* sanity check */
420  if(buffer != "pattern")
421  {
422  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
423  return false;
424  }
425  mmHeader.ctype=pattern;
426  break;
427  default:
428  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
429  return false;
430  }
431 
432  if(lineFeed(file))
433  return false;
434 
435  file >> buffer;
436 
437  std::transform(buffer.begin(), buffer.end(), buffer.begin(),
438  tolower);
439  switch(buffer[0])
440  {
441  case 'g':
442  /* sanity check */
443  if(buffer != "general")
444  {
445  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
446  return false;
447  }
448  mmHeader.structure=general;
449  break;
450  case 'h':
451  /* sanity check */
452  if(buffer != "hermitian")
453  {
454  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
455  return false;
456  }
457  mmHeader.structure=hermitian;
458  break;
459  case 's':
460  if(buffer.size()==1){
461  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
462  return false;
463  }
464 
465  switch(buffer[1])
466  {
467  case 'y':
468  /* sanity check */
469  if(buffer != "symmetric")
470  {
471  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
472  return false;
473  }
474  mmHeader.structure=symmetric;
475  break;
476  case 'k':
477  /* sanity check */
478  if(buffer != "skew-symmetric")
479  {
480  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
481  return false;
482  }
483  mmHeader.structure=skew_symmetric;
484  break;
485  default:
486  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
487  return false;
488  }
489  default:
490  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
491  return false;
492  }
493  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
494  c=file.peek();
495  return true;
496 
497  }
498 
499  void readNextLine(std::istream& file, std::ostringstream& line, LineType& type)
500  {
501  char c;
502  std::size_t index=0;
503 
504  //empty lines will be disgarded and we will simply read the next line
505  while(index==0&&!file.eof())
506  {
507  // strip spaces
508  while(!file.eof() && (c=file.get())==' ');
509 
510  //read the rest of the line until comment
511  while(!file.eof() && (c=file.get())=='\n'){
512  switch(c)
513  {
514  case '%':
515  /* disgard the rest of the line */
516  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
517  }
518  }
519  }
520 
521  // buffer[index]='\0';
522  }
523 
524  template<std::size_t brows, std::size_t bcols>
525  Dune::tuple<std::size_t, std::size_t, std::size_t>
526  calculateNNZ(std::size_t rows, std::size_t cols, std::size_t entries, const MMHeader& header)
527  {
528  std::size_t blockrows=rows/brows;
529  std::size_t blockcols=cols/bcols;
530  std::size_t blocksize=brows*bcols;
531  std::size_t blockentries=0;
532 
533  switch(header.structure)
534  {
535  case general:
536  blockentries = entries/blocksize; break;
537  case skew_symmetric:
538  blockentries = 2*entries/blocksize; break;
539  case symmetric:
540  blockentries = (2*entries-rows)/blocksize; break;
541  case hermitian:
542  blockentries = (2*entries-rows)/blocksize; break;
543  default:
544  throw Dune::NotImplemented();
545  }
546  return Dune::make_tuple(blockrows, blockcols, blockentries);
547  }
548 
549  /*
550  * @brief Storage class for the column index and the numeric value.
551  *
552  * \tparam T Either a NumericWrapper of the numeric type or PatternDummy
553  * for MatrixMarket pattern case.
554  */
555  template<typename T>
556  struct IndexData : public T
557  {
558  std::size_t index;
559  };
560 
561 
572  template<typename T>
573  struct NumericWrapper
574  {
576  operator T&()
577  {
578  return number;
579  }
580  };
581 
585  struct PatternDummy
586  {};
587 
588  template<>
589  struct NumericWrapper<PatternDummy>
590  {};
591 
592  template<typename T>
593  std::istream& operator>>(std::istream& is, NumericWrapper<T>& num)
594  {
595  return is>>num.number;
596  }
597 
598  std::istream& operator>>(std::istream& is, NumericWrapper<PatternDummy>& num)
599  {
600  return is;
601  }
602 
608  template<typename T>
609  bool operator<(const IndexData<T>& i1, const IndexData<T>& i2)
610  {
611  return i1.index<i2.index;
612  }
613 
619  template<typename T>
620  std::istream& operator>>(std::istream& is, IndexData<T>& data)
621  {
622  is>>data.index;
623  /* MatrixMarket indices are one based. Decrement for C++ */
624  --data.index;
625  return is>>data.number;
626  }
627 
634  template<typename D, int brows, int bcols>
635  struct MatrixValuesSetter
636  {
642  template<typename M>
643  void operator()(const std::vector<std::set<IndexData<D> > >& rows,
644  M& matrix)
645  {
646  for(typename M::RowIterator iter=matrix.begin();
647  iter!= matrix.end(); ++iter)
648  {
649  for(typename M::size_type brow=iter.index()*brows,
650  browend=iter.index()*brows+brows;
651  brow<browend; ++brow)
652  {
653  typedef typename std::set<IndexData<D> >::const_iterator Siter;
654  for(Siter siter=rows[brow].begin(), send=rows[brow].end();
655  siter != send; ++siter)
656  (*iter)[siter->index/bcols][brow%brows][siter->index%bcols]=siter->number;
657  }
658  }
659  }
660  };
661 
662  template<int brows, int bcols>
663  struct MatrixValuesSetter<PatternDummy,brows,bcols>
664  {
665  template<typename M>
666  void operator()(const std::vector<std::set<IndexData<PatternDummy> > >& rows,
667  M& matrix)
668  {}
669  };
670 
671  template<typename T, typename A, int brows, int bcols, typename D>
672  void readSparseEntries(Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A>& matrix,
673  std::istream& file, std::size_t entries,
674  const MMHeader& mmHeader, const D&)
675  {
677  // First path
678  // store entries together with column index in a speparate
679  // data structure
680  std::vector<std::set<IndexData<D> > > rows(matrix.N()*brows);
681 
682  for(; entries>0;--entries){
683  std::size_t row;
684  IndexData<D> data;
685  skipComments(file);
686  file>>row;
687  --row; // Index was 1 based.
688  assert(row/bcols<matrix.N());
689  file>>data;
690  assert(data.index/bcols<matrix.M());
691  rows[row].insert(data);
692  }
693 
694  // TODO extend to capture the nongeneral cases.
695  if(mmHeader.structure!= general)
696  DUNE_THROW(Dune::NotImplemented, "Only general is supported right now!");
697 
698  // Setup the matrix sparsity pattern
699  int nnz=0;
700  for(typename Matrix::CreateIterator iter=matrix.createbegin();
701  iter!= matrix.createend(); ++iter)
702  {
703  for(std::size_t brow=iter.index()*brows, browend=iter.index()*brows+brows;
704  brow<browend; ++brow)
705  {
706  typedef typename std::set<IndexData<D> >::const_iterator Siter;
707  for(Siter siter=rows[brow].begin(), send=rows[brow].end();
708  siter != send; ++siter, ++nnz)
709  iter.insert(siter->index/bcols);
710  }
711  }
712 
713  //Set the matrix values
714  matrix=0;
715 
716  MatrixValuesSetter<D,brows,bcols> Setter;
717 
718  Setter(rows, matrix);
719  }
720  } // end anonymous namespace
721 
722  class MatrixMarketFormatError : public Dune::Exception
723  {};
724 
725 
726  void mm_read_header(std::size_t& rows, std::size_t& cols, MMHeader& header, std::istream& istr,
727  bool isVector)
728  {
729  if(!readMatrixMarketBanner(istr, header)){
730  std::cerr << "First line was not a correct Matrix Market banner. Using default:\n"
731  << "%%MatrixMarket matrix coordinate real general"<<std::endl;
732  // Go to the beginning of the file
733  istr.clear() ;
734  istr.seekg(0, std::ios::beg);
735  if(isVector)
736  header.type=array_type;
737  }
738 
739  skipComments(istr);
740 
741  if(lineFeed(istr))
742  throw MatrixMarketFormatError();
743 
744  istr >> rows;
745 
746  if(lineFeed(istr))
747  throw MatrixMarketFormatError();
748  istr >> cols;
749  }
750 
751  template<typename T, typename A, int entries>
752  void mm_read_vector_entries(Dune::BlockVector<Dune::FieldVector<T,entries>,A>& vector,
753  std::size_t size,
754  std::istream& istr)
755  {
756  for(int i=0;size>0; ++i, --size){
757  T val;
758  istr>>val;
759  vector[i/entries][i%entries]=val;
760  }
761  }
762 
763 
770  template<typename T, typename A, int entries>
771  void readMatrixMarket(Dune::BlockVector<Dune::FieldVector<T,entries>,A>& vector,
772  std::istream& istr)
773  {
774  MMHeader header;
775  std::size_t rows, cols;
776  mm_read_header(rows,cols,header,istr, true);
777  if(cols!=1)
778  DUNE_THROW(MatrixMarketFormatError, "cols!=1, therefore this is no vector!");
779 
780  if(header.type!=array_type)
781  DUNE_THROW(MatrixMarketFormatError, "Vectors have to be stored in array format!");
782 
783  std::size_t size=rows/entries;
784  if(size*entries!=rows)
785  DUNE_THROW(MatrixMarketFormatError, "Block size of vector is not correct1");
786 
787  vector.resize(size);
788 
789  istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
790 
791  mm_read_vector_entries(vector, rows, istr);
792  }
793 
794 
801  template<typename T, typename A, int brows, int bcols>
803  std::istream& istr)
804  {
805 
807 
808  MMHeader header;
809  if(!readMatrixMarketBanner(istr, header)){
810  std::cerr << "First line was not a correct Matrix Market banner. Using default:\n"
811  << "%%MatrixMarket matrix coordinate real general"<<std::endl;
812  // Go to the beginning of the file
813  istr.clear() ;
814  istr.seekg(0, std::ios::beg);
815  }
816  skipComments(istr);
817 
818  std::size_t rows, cols, entries;
819 
820  if(lineFeed(istr))
821  throw MatrixMarketFormatError();
822 
823  istr >> rows;
824 
825  if(lineFeed(istr))
826  throw MatrixMarketFormatError();
827  istr >> cols;
828 
829  if(lineFeed(istr))
830  throw MatrixMarketFormatError();
831 
832  istr >>entries;
833 
834  std::size_t nnz, blockrows, blockcols;
835 
836  Dune::tie(blockrows, blockcols, nnz) = calculateNNZ<brows, bcols>(rows, cols, entries, header);
837 
838  istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
839 
840 
841  matrix.setSize(blockrows, blockcols);
842  matrix.setBuildMode(Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A>::row_wise);
843 
844  if(header.type==array_type)
845  DUNE_THROW(Dune::NotImplemented, "Array format currently not supported for matrices!");
846 
847  readSparseEntries(matrix, istr, entries, header, NumericWrapper<double>());
848  }
849 
850  template<typename M>
852  {};
853 
854  template<typename B, int i, int j, typename A>
856  {
857  enum{
858  rows = i,
859  cols = j
860  };
861  };
862 
863  template<typename B, int i, int j>
864  void mm_print_entry(const FieldMatrix<B,i,j>& entry,
865  typename FieldMatrix<B,i,j>::size_type rowidx,
866  typename FieldMatrix<B,i,j>::size_type colidx,
867  std::ostream& ostr)
868  {
869  typedef typename FieldMatrix<B,i,j>::const_iterator riterator;
870  typedef typename FieldMatrix<B,i,j>::ConstColIterator citerator;
871  for(riterator row=entry.begin(); row != entry.end(); ++row, ++rowidx){
872  int coli=colidx;
873  for(citerator col = row->begin(); col != row->end(); ++col, ++coli)
874  ostr<< rowidx<<" "<<coli<<" "<<*col<<std::endl;
875  }
876  }
877 
878  // Write a vector entry
879  template<typename V>
880  void mm_print_vector_entry(const V& entry, std::ostream& ostr,
881  const integral_constant<int,1>&)
882  {
883  ostr<<entry<<std::endl;
884  }
885 
886  // Write a vector
887  template<typename V>
888  void mm_print_vector_entry(const V& vector, std::ostream& ostr,
889  const integral_constant<int,0>&)
890  {
891  // Is the entry a supported numeric type?
892  const int isnumeric = mm_numeric_type<typename V::block_type>::is_numeric;
893  typedef typename V::const_iterator VIter;
894 
895  for(VIter i=vector.begin(); i != vector.end(); ++i)
896 
897  mm_print_vector_entry(*i, ostr,
898  integral_constant<int,isnumeric>());
899  }
900 
901  template<typename T, typename A, int i>
902  std::size_t countEntries(const BlockVector<FieldVector<T,i>,A>& vector)
903  {
904  return vector.size()*i;
905  }
906 
907  // Version for writing vectors.
908  template<typename V>
909  void writeMatrixMarket(const V& vector, std::ostream& ostr,
910  const integral_constant<int,0>&)
911  {
912  ostr<<countEntries(vector)<<" "<<1<<std::endl;
913  const int isnumeric = mm_numeric_type<typename V::block_type>::is_numeric;
914  mm_print_vector_entry(vector,ostr, integral_constant<int,isnumeric>());
915  }
916 
917  // Versions for writing matrices
918  template<typename M>
919  void writeMatrixMarket(const M& matrix,
920  std::ostream& ostr,
921  const integral_constant<int,1>&)
922  {
923  ostr<<matrix.M()*mm_multipliers<M>::rows<<" "
924  <<matrix.N()*mm_multipliers<M>::cols<<" "
925  <<countNonZeros(matrix)<<std::endl;
926 
927  typedef typename M::const_iterator riterator;
928  typedef typename M::ConstColIterator citerator;
929  for(riterator row=matrix.begin(); row != matrix.end(); ++row)
930  for(citerator col = row->begin(); col != row->end(); ++col)
931  // Matrix Market indexing start with 1!
932  mm_print_entry(*col, row.index()*mm_multipliers<M>::rows+1,
933  col.index()*mm_multipliers<M>::cols+1, ostr);
934  }
935 
936 
940  template<typename M>
941  void writeMatrixMarket(const M& matrix,
942  std::ostream& ostr)
943  {
944  // Write header information
945  mm_header_printer<M>::print(ostr);
946  mm_block_structure_header<M>::print(ostr,matrix);
947  // Choose the correct function for matrix and vector
948  writeMatrixMarket(matrix,ostr,integral_constant<int,IsMatrix<M>::value>());
949  }
950 
951 
962  template<typename M>
963  void storeMatrixMarket(const M& matrix,
964  std::string filename)
965  {
966  std::ofstream file(filename.c_str());
967  file.setf(std::ios::scientific,std::ios::floatfield);
968  writeMatrixMarket(matrix, file);
969  file.close();
970  }
971 
972 #if HAVE_MPI
973 
986  template<typename M, typename G, typename L>
987  void storeMatrixMarket(const M& matrix,
988  std::string filename,
990  bool storeIndices=true)
991  {
992  // Get our rank
993  int rank = comm.communicator().rank();
994  // Write the local matrix
995  std::ostringstream rfilename;
996  rfilename<<filename <<"_"<<rank<<".mm";
997  std::cout<< rfilename.str()<<std::endl;
998  std::ofstream file(rfilename.str().c_str());
999  file.setf(std::ios::scientific,std::ios::floatfield);
1000  writeMatrixMarket(matrix, file);
1001  file.close();
1002 
1003  if(!storeIndices)
1004  return;
1005 
1006  // Write the global to local index mapping
1007  rfilename.str("");
1008  rfilename<<filename<<"_"<<rank<<".idx";
1009  file.open(rfilename.str().c_str());
1010  file.setf(std::ios::scientific,std::ios::floatfield);
1011  typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet IndexSet;
1012  typedef typename IndexSet::const_iterator Iterator;
1013  for(Iterator iter = comm.indexSet().begin();
1014  iter != comm.indexSet().end(); ++iter){
1015  file << iter->global()<<" "<<(std::size_t)iter->local()<<" "
1016  <<(int)iter->local().attribute()<<" "<<(int)iter->local().isPublic()<<std::endl;
1017  }
1018  // Store neighbour information for efficient remote indices setup.
1019  file<<"neighbours:";
1020  const std::set<int>& neighbours=comm.remoteIndices().getNeighbours();
1021  typedef std::set<int>::const_iterator SIter;
1022  for(SIter neighbour=neighbours.begin(); neighbour != neighbours.end(); ++neighbour){
1023  file<<" "<< *neighbour;
1024  }
1025  file.close();
1026  }
1027 
1042  template<typename M, typename G, typename L>
1043  void loadMatrixMarket(M& matrix,
1044  const std::string& filename,
1046  bool readIndices=true)
1047  {
1049  typedef typename LocalIndex::Attribute Attribute;
1050  // Get our rank
1051  int rank = comm.communicator().rank();
1052  // load local matrix
1053  std::ostringstream rfilename;
1054  rfilename<<filename <<"_"<<rank<<".mm";
1055  std::ifstream file;
1056  file.open(rfilename.str().c_str(), std::ios::in);
1057  if(!file)
1058  DUNE_THROW(IOError, "Could not open file" << rfilename.str().c_str());
1059  //if(!file.is_open())
1060  //
1061  readMatrixMarket(matrix,file);
1062  file.close();
1063 
1064  if(!readIndices)
1065  return;
1066 
1067  // read indices
1068  typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet IndexSet;
1069  IndexSet& pis=comm.pis;
1070  rfilename.str("");
1071  rfilename<<filename<<"_"<<rank<<".idx";
1072  file.open(rfilename.str().c_str());
1073  if(pis.size()!=0)
1074  DUNE_THROW(InvalidIndexSetState, "Index set is not empty!");
1075 
1076  pis.beginResize();
1077  while(!file.eof() && file.peek()!='n'){
1078  G g;
1079  file >>g;
1080  std::size_t l;
1081  file >>l;
1082  int c;
1083  file >>c;
1084  bool b;
1085  file >> b;
1086  pis.add(g,LocalIndex(l,Attribute(c),b));
1087  lineFeed(file);
1088  }
1089  pis.endResize();
1090  if(!file.eof()){
1091  // read neighbours
1092  std::string s;
1093  file>>s;
1094  if(s!="neighbours:")
1095  DUNE_THROW(MatrixMarketFormatError, "was exspecting the string: \"neighbours:\"");
1096  std::set<int> nb;
1097  while(!file.eof()){
1098  int i;
1099  file >> i;
1100  nb.insert(i);
1101  }
1102  file.close();
1103  comm.ri.setNeighbours(nb);
1104  }
1105  comm.ri.template rebuild<false>();
1106  }
1107 
1108  #endif
1109 
1120  template<typename M>
1121  void loadMatrixMarket(M& matrix,
1122  const std::string& filename)
1123  {
1124  std::ifstream file;
1125  file.open(filename.c_str(), std::ios::in);
1126  if(!file)
1127  DUNE_THROW(IOError, "Could not open file" << filename);
1128  readMatrixMarket(matrix,file);
1129  file.close();
1130  }
1131 
1133 }
1134 #endif