Dune Core Modules (2.6.0)

matrixmarket.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ISTL_MATRIXMARKET_HH
4 #define DUNE_ISTL_MATRIXMARKET_HH
5 
6 #include <ostream>
7 #include <istream>
8 #include <fstream>
9 #include <sstream>
10 #include <limits>
11 #include <ios>
12 #include <tuple>
13 #include "matrixutils.hh"
14 #include "bcrsmatrix.hh"
15 #include "owneroverlapcopy.hh"
16 #include <dune/common/fmatrix.hh>
17 #include <dune/common/unused.hh>
18 
19 namespace Dune
20 {
21 
47  namespace MatrixMarketImpl
48  {
58  template<class T>
59  struct mm_numeric_type {
60  enum {
64  is_numeric=false
65  };
66  };
67 
68  template<>
69  struct mm_numeric_type<int>
70  {
71  enum {
75  is_numeric=true
76  };
77 
78  static std::string str()
79  {
80  return "integer";
81  }
82  };
83 
84  template<>
85  struct mm_numeric_type<double>
86  {
87  enum {
91  is_numeric=true
92  };
93 
94  static std::string str()
95  {
96  return "real";
97  }
98  };
99 
100  template<>
101  struct mm_numeric_type<float>
102  {
103  enum {
107  is_numeric=true
108  };
109 
110  static std::string str()
111  {
112  return "real";
113  }
114  };
115 
116  template<>
117  struct mm_numeric_type<std::complex<double> >
118  {
119  enum {
123  is_numeric=true
124  };
125 
126  static std::string str()
127  {
128  return "complex";
129  }
130  };
131 
132  template<>
133  struct mm_numeric_type<std::complex<float> >
134  {
135  enum {
139  is_numeric=true
140  };
141 
142  static std::string str()
143  {
144  return "complex";
145  }
146  };
147 
156  template<class M>
158 
159  template<typename T, typename A, int i, int j>
160  struct mm_header_printer<BCRSMatrix<FieldMatrix<T,i,j>,A> >
161  {
162  static void print(std::ostream& os)
163  {
164  os<<"%%MatrixMarket matrix coordinate ";
165  os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
166  }
167  };
168 
169  template<typename B, typename A>
170  struct mm_header_printer<BlockVector<B,A> >
171  {
172  static void print(std::ostream& os)
173  {
174  os<<"%%MatrixMarket matrix array ";
175  os<<mm_numeric_type<typename B::field_type>::str()<<" general"<<std::endl;
176  }
177  };
178 
179  template<typename T, int j>
180  struct mm_header_printer<FieldVector<T,j> >
181  {
182  static void print(std::ostream& os)
183  {
184  os<<"%%MatrixMarket matrix array ";
185  os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
186  }
187  };
188 
189  template<typename T, int i, int j>
190  struct mm_header_printer<FieldMatrix<T,i,j> >
191  {
192  static void print(std::ostream& os)
193  {
194  os<<"%%MatrixMarket matrix array ";
195  os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
196  }
197  };
198 
207  template<class M>
209 
210  template<typename T, typename A, int i>
212  {
213  typedef BlockVector<FieldVector<T,i>,A> M;
214 
215  static void print(std::ostream& os, const M&)
216  {
217  os<<"% ISTL_STRUCT blocked ";
218  os<<i<<" "<<1<<std::endl;
219  }
220  };
221 
222  template<typename T, typename A, int i, int j>
223  struct mm_block_structure_header<BCRSMatrix<FieldMatrix<T,i,j>,A> >
224  {
225  typedef BCRSMatrix<FieldMatrix<T,i,j>,A> M;
226 
227  static void print(std::ostream& os, const M&)
228  {
229  os<<"% ISTL_STRUCT blocked ";
230  os<<i<<" "<<j<<std::endl;
231  }
232  };
233 
234 
235  template<typename T, int i, int j>
236  struct mm_block_structure_header<FieldMatrix<T,i,j> >
237  {
238  typedef FieldMatrix<T,i,j> M;
239 
240  static void print(std::ostream& os, const M& m)
241  {}
242  };
243 
244  template<typename T, int i>
245  struct mm_block_structure_header<FieldVector<T,i> >
246  {
247  typedef FieldVector<T,i> M;
248 
249  static void print(std::ostream& os, const M& m)
250  {}
251  };
252 
253  enum LineType { MM_HEADER, MM_ISTLSTRUCT, DATA };
254  enum { MM_MAX_LINE_LENGTH=1025 };
255 
256  enum MM_TYPE { coordinate_type, array_type, unknown_type };
257 
258  enum MM_CTYPE { integer_type, double_type, complex_type, pattern, unknown_ctype };
259 
260  enum MM_STRUCTURE { general, symmetric, skew_symmetric, hermitian, unknown_structure };
261 
262  struct MMHeader
263  {
264  MMHeader()
265  : type(coordinate_type), ctype(double_type), structure(general)
266  {}
267  MM_TYPE type;
268  MM_CTYPE ctype;
269  MM_STRUCTURE structure;
270  };
271 
272  inline bool lineFeed(std::istream& file)
273  {
274  char c;
275  if(!file.eof())
276  c=file.peek();
277  else
278  return false;
279  // ignore whitespace
280  while(c==' ')
281  {
282  file.get();
283  if(file.eof())
284  return false;
285  c=file.peek();
286  }
287 
288  if(c=='\n') {
289  /* eat the line feed */
290  file.get();
291  return true;
292  }
293  return false;
294  }
295 
296  inline void skipComments(std::istream& file)
297  {
298  lineFeed(file);
299  char c=file.peek();
300  // ignore comment lines
301  while(c=='%')
302  {
303  /* discard the rest of the line */
304  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
305  c=file.peek();
306  }
307  }
308 
309 
310  inline bool readMatrixMarketBanner(std::istream& file, MMHeader& mmHeader)
311  {
312  std::string buffer;
313  char c;
314  file >> buffer;
315  c=buffer[0];
316  mmHeader=MMHeader();
317  if(c!='%')
318  return false;
319  std::cout<<buffer<<std::endl;
320  /* read the banner */
321  if(buffer!="%%MatrixMarket") {
322  /* discard the rest of the line */
323  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
324  return false;
325  }
326 
327  if(lineFeed(file))
328  /* premature end of line */
329  return false;
330 
331  /* read the matrix_type */
332  file >> buffer;
333 
334  if(buffer != "matrix")
335  {
336  /* discard the rest of the line */
337  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
338  return false;
339  }
340 
341  if(lineFeed(file))
342  /* premature end of line */
343  return false;
344 
345  /* The type of the matrix */
346  file >> buffer;
347 
348  if(buffer.empty())
349  return false;
350 
351  std::transform(buffer.begin(), buffer.end(), buffer.begin(),
352  ::tolower);
353 
354  switch(buffer[0])
355  {
356  case 'a' :
357  /* sanity check */
358  if(buffer != "array")
359  {
360  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
361  return false;
362  }
363  mmHeader.type=array_type;
364  break;
365  case 'c' :
366  /* sanity check */
367  if(buffer != "coordinate")
368  {
369  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
370  return false;
371  }
372  mmHeader.type=coordinate_type;
373  break;
374  default :
375  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
376  return false;
377  }
378 
379  if(lineFeed(file))
380  /* premature end of line */
381  return false;
382 
383  /* The numeric type used. */
384  file >> buffer;
385 
386  if(buffer.empty())
387  return false;
388 
389  std::transform(buffer.begin(), buffer.end(), buffer.begin(),
390  ::tolower);
391  switch(buffer[0])
392  {
393  case 'i' :
394  /* sanity check */
395  if(buffer != "integer")
396  {
397  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
398  return false;
399  }
400  mmHeader.ctype=integer_type;
401  break;
402  case 'r' :
403  /* sanity check */
404  if(buffer != "real")
405  {
406  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
407  return false;
408  }
409  mmHeader.ctype=double_type;
410  break;
411  case 'c' :
412  /* sanity check */
413  if(buffer != "complex")
414  {
415  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
416  return false;
417  }
418  mmHeader.ctype=complex_type;
419  break;
420  case 'p' :
421  /* sanity check */
422  if(buffer != "pattern")
423  {
424  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
425  return false;
426  }
427  mmHeader.ctype=pattern;
428  break;
429  default :
430  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
431  return false;
432  }
433 
434  if(lineFeed(file))
435  return false;
436 
437  file >> buffer;
438 
439  std::transform(buffer.begin(), buffer.end(), buffer.begin(),
440  ::tolower);
441  switch(buffer[0])
442  {
443  case 'g' :
444  /* sanity check */
445  if(buffer != "general")
446  {
447  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
448  return false;
449  }
450  mmHeader.structure=general;
451  break;
452  case 'h' :
453  /* sanity check */
454  if(buffer != "hermitian")
455  {
456  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
457  return false;
458  }
459  mmHeader.structure=hermitian;
460  break;
461  case 's' :
462  if(buffer.size()==1) {
463  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
464  return false;
465  }
466 
467  switch(buffer[1])
468  {
469  case 'y' :
470  /* sanity check */
471  if(buffer != "symmetric")
472  {
473  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
474  return false;
475  }
476  mmHeader.structure=symmetric;
477  break;
478  case 'k' :
479  /* sanity check */
480  if(buffer != "skew-symmetric")
481  {
482  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
483  return false;
484  }
485  mmHeader.structure=skew_symmetric;
486  break;
487  default :
488  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
489  return false;
490  }
491  default :
492  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
493  return false;
494  }
495  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
496  c=file.peek();
497  return true;
498 
499  }
500 
501  template<std::size_t brows, std::size_t bcols>
502  std::tuple<std::size_t, std::size_t, std::size_t>
503  calculateNNZ(std::size_t rows, std::size_t cols, std::size_t entries, const MMHeader& header)
504  {
505  std::size_t blockrows=rows/brows;
506  std::size_t blockcols=cols/bcols;
507  std::size_t blocksize=brows*bcols;
508  std::size_t blockentries=0;
509 
510  switch(header.structure)
511  {
512  case general :
513  blockentries = entries/blocksize; break;
514  case skew_symmetric :
515  blockentries = 2*entries/blocksize; break;
516  case symmetric :
517  blockentries = (2*entries-rows)/blocksize; break;
518  case hermitian :
519  blockentries = (2*entries-rows)/blocksize; break;
520  default :
521  throw Dune::NotImplemented();
522  }
523  return std::make_tuple(blockrows, blockcols, blockentries);
524  }
525 
526  /*
527  * @brief Storage class for the column index and the numeric value.
528  *
529  * \tparam T Either a NumericWrapper of the numeric type or PatternDummy
530  * for MatrixMarket pattern case.
531  */
532  template<typename T>
533  struct IndexData : public T
534  {
535  std::size_t index;
536  };
537 
538 
549  template<typename T>
551  {
552  T number;
553  operator T&()
554  {
555  return number;
556  }
557  };
558 
563  {};
564 
565  template<>
567  {};
568 
569  template<typename T>
570  std::istream& operator>>(std::istream& is, NumericWrapper<T>& num)
571  {
572  return is>>num.number;
573  }
574 
575  inline std::istream& operator>>(std::istream& is, NumericWrapper<PatternDummy>& num)
576  {
578  return is;
579  }
580 
586  template<typename T>
587  bool operator<(const IndexData<T>& i1, const IndexData<T>& i2)
588  {
589  return i1.index<i2.index;
590  }
591 
597  template<typename T>
598  std::istream& operator>>(std::istream& is, IndexData<T>& data)
599  {
600  is>>data.index;
601  /* MatrixMarket indices are one based. Decrement for C++ */
602  --data.index;
603  return is>>data.number;
604  }
605 
612  template<typename D, int brows, int bcols>
614  {
620  template<typename M>
621  void operator()(const std::vector<std::set<IndexData<D> > >& rows,
622  M& matrix)
623  {
624  for(typename M::RowIterator iter=matrix.begin();
625  iter!= matrix.end(); ++iter)
626  {
627  for(typename M::size_type brow=iter.index()*brows,
628  browend=iter.index()*brows+brows;
629  brow<browend; ++brow)
630  {
631  typedef typename std::set<IndexData<D> >::const_iterator Siter;
632  for(Siter siter=rows[brow].begin(), send=rows[brow].end();
633  siter != send; ++siter)
634  (*iter)[siter->index/bcols][brow%brows][siter->index%bcols]=siter->number;
635  }
636  }
637  }
638  };
639 
640  template<int brows, int bcols>
641  struct MatrixValuesSetter<PatternDummy,brows,bcols>
642  {
643  template<typename M>
644  void operator()(const std::vector<std::set<IndexData<PatternDummy> > >& rows,
645  M& matrix)
646  {}
647  };
648 
649  template<typename T, typename A, int brows, int bcols, typename D>
650  void readSparseEntries(Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A>& matrix,
651  std::istream& file, std::size_t entries,
652  const MMHeader& mmHeader, const D&)
653  {
655  // First path
656  // store entries together with column index in a speparate
657  // data structure
658  std::vector<std::set<IndexData<D> > > rows(matrix.N()*brows);
659 
660  for(; entries>0; --entries) {
661  std::size_t row;
662  IndexData<D> data;
663  skipComments(file);
664  file>>row;
665  --row; // Index was 1 based.
666  assert(row/bcols<matrix.N());
667  file>>data;
668  assert(data.index/bcols<matrix.M());
669  rows[row].insert(data);
670  }
671 
672  // TODO extend to capture the nongeneral cases.
673  if(mmHeader.structure!= general)
674  DUNE_THROW(Dune::NotImplemented, "Only general is supported right now!");
675 
676  // Setup the matrix sparsity pattern
677  int nnz=0;
678  for(typename Matrix::CreateIterator iter=matrix.createbegin();
679  iter!= matrix.createend(); ++iter)
680  {
681  for(std::size_t brow=iter.index()*brows, browend=iter.index()*brows+brows;
682  brow<browend; ++brow)
683  {
684  typedef typename std::set<IndexData<D> >::const_iterator Siter;
685  for(Siter siter=rows[brow].begin(), send=rows[brow].end();
686  siter != send; ++siter, ++nnz)
687  iter.insert(siter->index/bcols);
688  }
689  }
690 
691  //Set the matrix values
692  matrix=0;
693 
694  MatrixValuesSetter<D,brows,bcols> Setter;
695 
696  Setter(rows, matrix);
697  }
698  } // end namespace MatrixMarketImpl
699 
700  class MatrixMarketFormatError : public Dune::Exception
701  {};
702 
703 
704  inline void mm_read_header(std::size_t& rows, std::size_t& cols,
705  MatrixMarketImpl::MMHeader& header, std::istream& istr,
706  bool isVector)
707  {
708  using namespace MatrixMarketImpl;
709 
710  if(!readMatrixMarketBanner(istr, header)) {
711  std::cerr << "First line was not a correct Matrix Market banner. Using default:\n"
712  << "%%MatrixMarket matrix coordinate real general"<<std::endl;
713  // Go to the beginning of the file
714  istr.clear() ;
715  istr.seekg(0, std::ios::beg);
716  if(isVector)
717  header.type=array_type;
718  }
719 
720  skipComments(istr);
721 
722  if(lineFeed(istr))
723  throw MatrixMarketFormatError();
724 
725  istr >> rows;
726 
727  if(lineFeed(istr))
728  throw MatrixMarketFormatError();
729  istr >> cols;
730  }
731 
732  template<typename T, typename A, int entries>
733  void mm_read_vector_entries(Dune::BlockVector<Dune::FieldVector<T,entries>,A>& vector,
734  std::size_t size,
735  std::istream& istr)
736  {
737  for(int i=0; size>0; ++i, --size) {
738  T val;
739  istr>>val;
740  vector[i/entries][i%entries]=val;
741  }
742  }
743 
744 
751  template<typename T, typename A, int entries>
753  std::istream& istr)
754  {
755  using namespace MatrixMarketImpl;
756 
757  MMHeader header;
758  std::size_t rows, cols;
759  mm_read_header(rows,cols,header,istr, true);
760  if(cols!=1)
761  DUNE_THROW(MatrixMarketFormatError, "cols!=1, therefore this is no vector!");
762 
763  if(header.type!=array_type)
764  DUNE_THROW(MatrixMarketFormatError, "Vectors have to be stored in array format!");
765 
766  std::size_t size=rows/entries;
767  if(size*entries!=rows)
768  DUNE_THROW(MatrixMarketFormatError, "Block size of vector is not correct!");
769 
770  vector.resize(size);
771 
772  istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
773 
774  mm_read_vector_entries(vector, rows, istr);
775  }
776 
777 
784  template<typename T, typename A, int brows, int bcols>
786  std::istream& istr)
787  {
788  using namespace MatrixMarketImpl;
789 
790  MMHeader header;
791  if(!readMatrixMarketBanner(istr, header)) {
792  std::cerr << "First line was not a correct Matrix Market banner. Using default:\n"
793  << "%%MatrixMarket matrix coordinate real general"<<std::endl;
794  // Go to the beginning of the file
795  istr.clear() ;
796  istr.seekg(0, std::ios::beg);
797  }
798  skipComments(istr);
799 
800  std::size_t rows, cols, entries;
801 
802  if(lineFeed(istr))
803  throw MatrixMarketFormatError();
804 
805  istr >> rows;
806 
807  if(lineFeed(istr))
808  throw MatrixMarketFormatError();
809  istr >> cols;
810 
811  if(lineFeed(istr))
812  throw MatrixMarketFormatError();
813 
814  istr >>entries;
815 
816  std::size_t nnz, blockrows, blockcols;
817 
818  std::tie(blockrows, blockcols, nnz) = calculateNNZ<brows, bcols>(rows, cols, entries, header);
819 
820  istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
821 
822 
823  matrix.setSize(blockrows, blockcols);
824  matrix.setBuildMode(Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A>::row_wise);
825 
826  if(header.type==array_type)
827  DUNE_THROW(Dune::NotImplemented, "Array format currently not supported for matrices!");
828 
829  readSparseEntries(matrix, istr, entries, header, NumericWrapper<T>());
830  }
831 
832  template<typename M>
833  struct mm_multipliers
834  {};
835 
836  template<typename B, int i, int j, typename A>
837  struct mm_multipliers<BCRSMatrix<FieldMatrix<B,i,j>,A> >
838  {
839  enum {
840  rows = i,
841  cols = j
842  };
843  };
844 
845  template<typename B, int i, int j>
846  void mm_print_entry(const FieldMatrix<B,i,j>& entry,
847  typename FieldMatrix<B,i,j>::size_type rowidx,
848  typename FieldMatrix<B,i,j>::size_type colidx,
849  std::ostream& ostr)
850  {
851  typedef typename FieldMatrix<B,i,j>::const_iterator riterator;
852  typedef typename FieldMatrix<B,i,j>::ConstColIterator citerator;
853  for(riterator row=entry.begin(); row != entry.end(); ++row, ++rowidx) {
854  int coli=colidx;
855  for(citerator col = row->begin(); col != row->end(); ++col, ++coli)
856  ostr<< rowidx<<" "<<coli<<" "<<*col<<std::endl;
857  }
858  }
859 
860  // Write a vector entry
861  template<typename V>
862  void mm_print_vector_entry(const V& entry, std::ostream& ostr,
863  const std::integral_constant<int,1>&)
864  {
865  ostr<<entry<<std::endl;
866  }
867 
868  // Write a vector
869  template<typename V>
870  void mm_print_vector_entry(const V& vector, std::ostream& ostr,
871  const std::integral_constant<int,0>&)
872  {
873  using namespace MatrixMarketImpl;
874 
875  // Is the entry a supported numeric type?
877  typedef typename V::const_iterator VIter;
878 
879  for(VIter i=vector.begin(); i != vector.end(); ++i)
880 
881  mm_print_vector_entry(*i, ostr,
882  std::integral_constant<int,isnumeric>());
883  }
884 
885  template<typename T, typename A, int i>
886  std::size_t countEntries(const BlockVector<FieldVector<T,i>,A>& vector)
887  {
888  return vector.size()*i;
889  }
890 
891  // Version for writing vectors.
892  template<typename V>
893  void writeMatrixMarket(const V& vector, std::ostream& ostr,
894  const std::integral_constant<int,0>&)
895  {
896  using namespace MatrixMarketImpl;
897 
898  ostr<<countEntries(vector)<<" "<<1<<std::endl;
900  mm_print_vector_entry(vector,ostr, std::integral_constant<int,isnumeric>());
901  }
902 
903  // Versions for writing matrices
904  template<typename M>
905  void writeMatrixMarket(const M& matrix,
906  std::ostream& ostr,
907  const std::integral_constant<int,1>&)
908  {
909  ostr<<matrix.N()*mm_multipliers<M>::rows<<" "
910  <<matrix.M()*mm_multipliers<M>::cols<<" "
911  <<countNonZeros(matrix)<<std::endl;
912 
913  typedef typename M::const_iterator riterator;
914  typedef typename M::ConstColIterator citerator;
915  for(riterator row=matrix.begin(); row != matrix.end(); ++row)
916  for(citerator col = row->begin(); col != row->end(); ++col)
917  // Matrix Market indexing start with 1!
918  mm_print_entry(*col, row.index()*mm_multipliers<M>::rows+1,
919  col.index()*mm_multipliers<M>::cols+1, ostr);
920  }
921 
922 
926  template<typename M>
927  void writeMatrixMarket(const M& matrix,
928  std::ostream& ostr)
929  {
930  using namespace MatrixMarketImpl;
931 
932  // Write header information
933  mm_header_printer<M>::print(ostr);
934  mm_block_structure_header<M>::print(ostr,matrix);
935  // Choose the correct function for matrix and vector
936  writeMatrixMarket(matrix,ostr,std::integral_constant<int,IsMatrix<M>::value>());
937  }
938 
939 
950  template<typename M>
951  void storeMatrixMarket(const M& matrix,
952  std::string filename)
953  {
954  std::ofstream file(filename.c_str());
955  file.setf(std::ios::scientific,std::ios::floatfield);
956  writeMatrixMarket(matrix, file);
957  file.close();
958  }
959 
960 #if HAVE_MPI
974  template<typename M, typename G, typename L>
975  void storeMatrixMarket(const M& matrix,
976  std::string filename,
978  bool storeIndices=true)
979  {
980  // Get our rank
981  int rank = comm.communicator().rank();
982  // Write the local matrix
983  std::ostringstream rfilename;
984  rfilename<<filename <<"_"<<rank<<".mm";
985  std::cout<< rfilename.str()<<std::endl;
986  std::ofstream file(rfilename.str().c_str());
987  file.setf(std::ios::scientific,std::ios::floatfield);
988  writeMatrixMarket(matrix, file);
989  file.close();
990 
991  if(!storeIndices)
992  return;
993 
994  // Write the global to local index mapping
995  rfilename.str("");
996  rfilename<<filename<<"_"<<rank<<".idx";
997  file.open(rfilename.str().c_str());
998  file.setf(std::ios::scientific,std::ios::floatfield);
1000  typedef typename IndexSet::const_iterator Iterator;
1001  for(Iterator iter = comm.indexSet().begin();
1002  iter != comm.indexSet().end(); ++iter) {
1003  file << iter->global()<<" "<<(std::size_t)iter->local()<<" "
1004  <<(int)iter->local().attribute()<<" "<<(int)iter->local().isPublic()<<std::endl;
1005  }
1006  // Store neighbour information for efficient remote indices setup.
1007  file<<"neighbours:";
1008  const std::set<int>& neighbours=comm.remoteIndices().getNeighbours();
1009  typedef std::set<int>::const_iterator SIter;
1010  for(SIter neighbour=neighbours.begin(); neighbour != neighbours.end(); ++neighbour) {
1011  file<<" "<< *neighbour;
1012  }
1013  file.close();
1014  }
1015 
1030  template<typename M, typename G, typename L>
1031  void loadMatrixMarket(M& matrix,
1032  const std::string& filename,
1034  bool readIndices=true)
1035  {
1036  using namespace MatrixMarketImpl;
1037 
1039  typedef typename LocalIndex::Attribute Attribute;
1040  // Get our rank
1041  int rank = comm.communicator().rank();
1042  // load local matrix
1043  std::ostringstream rfilename;
1044  rfilename<<filename <<"_"<<rank<<".mm";
1045  std::ifstream file;
1046  file.open(rfilename.str().c_str(), std::ios::in);
1047  if(!file)
1048  DUNE_THROW(IOError, "Could not open file: " << rfilename.str().c_str());
1049  //if(!file.is_open())
1050  //
1051  readMatrixMarket(matrix,file);
1052  file.close();
1053 
1054  if(!readIndices)
1055  return;
1056 
1057  // read indices
1059  IndexSet& pis=comm.pis;
1060  rfilename.str("");
1061  rfilename<<filename<<"_"<<rank<<".idx";
1062  file.open(rfilename.str().c_str());
1063  if(pis.size()!=0)
1064  DUNE_THROW(InvalidIndexSetState, "Index set is not empty!");
1065 
1066  pis.beginResize();
1067  while(!file.eof() && file.peek()!='n') {
1068  G g;
1069  file >>g;
1070  std::size_t l;
1071  file >>l;
1072  int c;
1073  file >>c;
1074  bool b;
1075  file >> b;
1076  pis.add(g,LocalIndex(l,Attribute(c),b));
1077  lineFeed(file);
1078  }
1079  pis.endResize();
1080  if(!file.eof()) {
1081  // read neighbours
1082  std::string s;
1083  file>>s;
1084  if(s!="neighbours:")
1085  DUNE_THROW(MatrixMarketFormatError, "was expecting the string: \"neighbours:\"");
1086  std::set<int> nb;
1087  while(!file.eof()) {
1088  int i;
1089  file >> i;
1090  nb.insert(i);
1091  }
1092  file.close();
1093  comm.ri.setNeighbours(nb);
1094  }
1095  comm.ri.template rebuild<false>();
1096  }
1097 
1098  #endif
1099 
1110  template<typename M>
1111  void loadMatrixMarket(M& matrix,
1112  const std::string& filename)
1113  {
1114  std::ifstream file;
1115  file.open(filename.c_str(), std::ios::in);
1116  if(!file)
1117  DUNE_THROW(IOError, "Could not open file: " << filename);
1118  readMatrixMarket(matrix,file);
1119  file.close();
1120  }
1121 
1123 }
1124 #endif
Implementation of the BCRSMatrix class.
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:423
A vector of blocks with memory management.
Definition: bvector.hh:317
int rank() const
Return rank, is between 0 and size()-1.
Definition: mpicollectivecommunication.hh:163
ConstIterator const_iterator
typedef for stl compliant access
Definition: densematrix.hh:267
std::remove_reference< const_row_reference >::type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition: densematrix.hh:271
Base class for Dune-Exceptions.
Definition: exceptions.hh:94
A dense n x m matrix.
Definition: fmatrix.hh:68
vector space out of a tensor product of fields.
Definition: fvector.hh:93
Default exception class for I/O errors.
Definition: exceptions.hh:229
Index Set Interface base class.
Definition: indexidset.hh:76
IndexType size(GeometryType type) const
Return total number of entities of given geometry type in entity set .
Definition: indexidset.hh:220
Exception indicating that the index set is not in the expected state.
Definition: indexset.hh:204
An index present on the local process.
Definition: localindex.hh:33
Default exception for dummy implementations.
Definition: exceptions.hh:261
A class setting up standard communication for a two-valued attribute set with owner/overlap/copy sema...
Definition: owneroverlapcopy.hh:172
const ParallelIndexSet & indexSet() const
Get the underlying parallel index set.
Definition: owneroverlapcopy.hh:463
const RemoteIndices & remoteIndices() const
Get the underlying remote indices.
Definition: owneroverlapcopy.hh:472
Manager class for the mapping between local indices and globally unique indices.
Definition: indexset.hh:217
Implements a matrix constructed from a given type representing a field and compile-time given number ...
iterator begin()
Get an iterator over the indices positioned at the first index.
iterator end()
Get an iterator over the indices positioned after the last index.
TL LocalIndex
The type of the local index, e.g. ParallelLocalIndex.
Definition: indexset.hh:238
Stream & operator>>(Stream &stream, std::tuple< Ts... > &t)
Read a std::tuple.
Definition: streamoperators.hh:41
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:25
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
void readMatrixMarket(Dune::BlockVector< Dune::FieldVector< T, entries >, A > &vector, std::istream &istr)
Reads a BlockVector from a matrix market file.
Definition: matrixmarket.hh:752
void loadMatrixMarket(M &matrix, const std::string &filename, OwnerOverlapCopyCommunication< G, L > &comm, bool readIndices=true)
Load a parallel matrix/vector stored in matrix market format.
Definition: matrixmarket.hh:1031
void storeMatrixMarket(const M &matrix, std::string filename)
Stores a parallel matrix/vector in matrix market format in a file.
Definition: matrixmarket.hh:951
int countNonZeros(const M &matrix)
Get the number of nonzero fields in the matrix.
Definition: matrixutils.hh:153
EnableIfInterOperable< T1, T2, bool >::type operator<(const RandomAccessIteratorFacade< T1, V1, R1, D > &lhs, const RandomAccessIteratorFacade< T2, V2, R2, D > &rhs)
Comparison operator.
Definition: iteratorfacades.hh:629
Some handy generic functions for ISTL matrices.
Dune namespace.
Definition: alignedallocator.hh:10
Classes providing communication interfaces for overlapping Schwarz methods.
Test whether a type is an ISTL Matrix.
Definition: matrixutils.hh:478
Functor to the data values of the matrix.
Definition: matrixmarket.hh:614
void operator()(const std::vector< std::set< IndexData< D > > > &rows, M &matrix)
Sets the matrixvalues.
Definition: matrixmarket.hh:621
a wrapper class of numeric values.
Definition: matrixmarket.hh:551
Utility class for marking the pattern type of the MatrixMarket matrices.
Definition: matrixmarket.hh:563
Metaprogram for writing the ISTL block structure header.
Definition: matrixmarket.hh:208
Meta program to write the correct Matrix Market header.
Definition: matrixmarket.hh:157
Helper metaprogram to get the matrix market string representation of the numeric type.
Definition: matrixmarket.hh:59
@ is_numeric
Whether T is a supported numeric type.
Definition: matrixmarket.hh:64
Definition of the DUNE_UNUSED macro for the case that config.h is not available.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.80.0 (Mar 28, 23:30, 2024)