dune-istl  2.7.1
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 <algorithm>
7 #include <complex>
8 #include <cstddef>
9 #include <fstream>
10 #include <ios>
11 #include <iostream>
12 #include <istream>
13 #include <limits>
14 #include <ostream>
15 #include <set>
16 #include <sstream>
17 #include <string>
18 #include <tuple>
19 #include <type_traits>
20 #include <vector>
21 
22 #include <dune/common/exceptions.hh>
23 #include <dune/common/fmatrix.hh>
24 #include <dune/common/fvector.hh>
25 #include <dune/common/unused.hh>
26 #include <dune/common/hybridutilities.hh>
27 #include <dune/common/stdstreams.hh>
28 
29 #include <dune/istl/bcrsmatrix.hh>
30 #include <dune/istl/bvector.hh>
31 #include <dune/istl/matrixutils.hh> // countNonZeros()
33 
34 namespace Dune
35 {
36 
62  namespace MatrixMarketImpl
63  {
73  template<class T>
74  struct mm_numeric_type {
75  enum {
79  is_numeric=false
80  };
81  };
82 
83  template<>
84  struct mm_numeric_type<int>
85  {
86  enum {
90  is_numeric=true
91  };
92 
93  static std::string str()
94  {
95  return "integer";
96  }
97  };
98 
99  template<>
100  struct mm_numeric_type<double>
101  {
102  enum {
106  is_numeric=true
107  };
108 
109  static std::string str()
110  {
111  return "real";
112  }
113  };
114 
115  template<>
116  struct mm_numeric_type<float>
117  {
118  enum {
122  is_numeric=true
123  };
124 
125  static std::string str()
126  {
127  return "real";
128  }
129  };
130 
131  template<>
132  struct mm_numeric_type<std::complex<double> >
133  {
134  enum {
138  is_numeric=true
139  };
140 
141  static std::string str()
142  {
143  return "complex";
144  }
145  };
146 
147  template<>
148  struct mm_numeric_type<std::complex<float> >
149  {
150  enum {
154  is_numeric=true
155  };
156 
157  static std::string str()
158  {
159  return "complex";
160  }
161  };
162 
171  template<class M>
173 
174  template<typename T, typename A>
176  {
177  static void print(std::ostream& os)
178  {
179  os<<"%%MatrixMarket matrix coordinate ";
180  os<<mm_numeric_type<typename Imp::BlockTraits<T>::field_type>::str()<<" general"<<std::endl;
181  }
182  };
183 
184  template<typename B, typename A>
186  {
187  static void print(std::ostream& os)
188  {
189  os<<"%%MatrixMarket matrix array ";
190  os<<mm_numeric_type<typename Imp::BlockTraits<B>::field_type>::str()<<" general"<<std::endl;
191  }
192  };
193 
194  template<typename T, int j>
195  struct mm_header_printer<FieldVector<T,j> >
196  {
197  static void print(std::ostream& os)
198  {
199  os<<"%%MatrixMarket matrix array ";
200  os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
201  }
202  };
203 
204  template<typename T, int i, int j>
206  {
207  static void print(std::ostream& os)
208  {
209  os<<"%%MatrixMarket matrix array ";
210  os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
211  }
212  };
213 
222  template<class M>
224 
225  template<typename T, typename A>
227  {
229  static_assert(IsNumber<T>::value, "Only scalar entries are expected here!");
230 
231  static void print(std::ostream& os, const M&)
232  {
233  os<<"% ISTL_STRUCT blocked ";
234  os<<"1 1"<<std::endl;
235  }
236  };
237 
238  template<typename T, typename A, int i>
239  struct mm_block_structure_header<BlockVector<FieldVector<T,i>,A> >
240  {
242 
243  static void print(std::ostream& os, const M&)
244  {
245  os<<"% ISTL_STRUCT blocked ";
246  os<<i<<" "<<1<<std::endl;
247  }
248  };
249 
250  template<typename T, typename A>
252  {
254  static_assert(IsNumber<T>::value, "Only scalar entries are expected here!");
255 
256  static void print(std::ostream& os, const M&)
257  {
258  os<<"% ISTL_STRUCT blocked ";
259  os<<"1 1"<<std::endl;
260  }
261  };
262 
263  template<typename T, typename A, int i, int j>
265  {
267 
268  static void print(std::ostream& os, const M&)
269  {
270  os<<"% ISTL_STRUCT blocked ";
271  os<<i<<" "<<j<<std::endl;
272  }
273  };
274 
275 
276  template<typename T, int i, int j>
278  {
280 
281  static void print(std::ostream& os, const M& m)
282  {}
283  };
284 
285  template<typename T, int i>
286  struct mm_block_structure_header<FieldVector<T,i> >
287  {
288  typedef FieldVector<T,i> M;
289 
290  static void print(std::ostream& os, const M& m)
291  {}
292  };
293 
295  enum { MM_MAX_LINE_LENGTH=1025 };
296 
298 
300 
302 
303  struct MMHeader
304  {
307  {}
311  };
312 
313  inline bool lineFeed(std::istream& file)
314  {
315  char c;
316  if(!file.eof())
317  c=file.peek();
318  else
319  return false;
320  // ignore whitespace
321  while(c==' ')
322  {
323  file.get();
324  if(file.eof())
325  return false;
326  c=file.peek();
327  }
328 
329  if(c=='\n') {
330  /* eat the line feed */
331  file.get();
332  return true;
333  }
334  return false;
335  }
336 
337  inline void skipComments(std::istream& file)
338  {
339  lineFeed(file);
340  char c=file.peek();
341  // ignore comment lines
342  while(c=='%')
343  {
344  /* discard the rest of the line */
345  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
346  c=file.peek();
347  }
348  }
349 
350 
351  inline bool readMatrixMarketBanner(std::istream& file, MMHeader& mmHeader)
352  {
353  std::string buffer;
354  char c;
355  file >> buffer;
356  c=buffer[0];
357  mmHeader=MMHeader();
358  if(c!='%')
359  return false;
360  dverb<<buffer<<std::endl;
361  /* read the banner */
362  if(buffer!="%%MatrixMarket") {
363  /* discard the rest of the line */
364  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
365  return false;
366  }
367 
368  if(lineFeed(file))
369  /* premature end of line */
370  return false;
371 
372  /* read the matrix_type */
373  file >> buffer;
374 
375  if(buffer != "matrix")
376  {
377  /* discard the rest of the line */
378  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
379  return false;
380  }
381 
382  if(lineFeed(file))
383  /* premature end of line */
384  return false;
385 
386  /* The type of the matrix */
387  file >> buffer;
388 
389  if(buffer.empty())
390  return false;
391 
392  std::transform(buffer.begin(), buffer.end(), buffer.begin(),
393  ::tolower);
394 
395  switch(buffer[0])
396  {
397  case 'a' :
398  /* sanity check */
399  if(buffer != "array")
400  {
401  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
402  return false;
403  }
404  mmHeader.type=array_type;
405  break;
406  case 'c' :
407  /* sanity check */
408  if(buffer != "coordinate")
409  {
410  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
411  return false;
412  }
413  mmHeader.type=coordinate_type;
414  break;
415  default :
416  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
417  return false;
418  }
419 
420  if(lineFeed(file))
421  /* premature end of line */
422  return false;
423 
424  /* The numeric type used. */
425  file >> buffer;
426 
427  if(buffer.empty())
428  return false;
429 
430  std::transform(buffer.begin(), buffer.end(), buffer.begin(),
431  ::tolower);
432  switch(buffer[0])
433  {
434  case 'i' :
435  /* sanity check */
436  if(buffer != "integer")
437  {
438  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
439  return false;
440  }
441  mmHeader.ctype=integer_type;
442  break;
443  case 'r' :
444  /* sanity check */
445  if(buffer != "real")
446  {
447  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
448  return false;
449  }
450  mmHeader.ctype=double_type;
451  break;
452  case 'c' :
453  /* sanity check */
454  if(buffer != "complex")
455  {
456  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
457  return false;
458  }
459  mmHeader.ctype=complex_type;
460  break;
461  case 'p' :
462  /* sanity check */
463  if(buffer != "pattern")
464  {
465  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
466  return false;
467  }
468  mmHeader.ctype=pattern;
469  break;
470  default :
471  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
472  return false;
473  }
474 
475  if(lineFeed(file))
476  return false;
477 
478  file >> buffer;
479 
480  std::transform(buffer.begin(), buffer.end(), buffer.begin(),
481  ::tolower);
482  switch(buffer[0])
483  {
484  case 'g' :
485  /* sanity check */
486  if(buffer != "general")
487  {
488  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
489  return false;
490  }
491  mmHeader.structure=general;
492  break;
493  case 'h' :
494  /* sanity check */
495  if(buffer != "hermitian")
496  {
497  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
498  return false;
499  }
500  mmHeader.structure=hermitian;
501  break;
502  case 's' :
503  if(buffer.size()==1) {
504  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
505  return false;
506  }
507 
508  switch(buffer[1])
509  {
510  case 'y' :
511  /* sanity check */
512  if(buffer != "symmetric")
513  {
514  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
515  return false;
516  }
517  mmHeader.structure=symmetric;
518  break;
519  case 'k' :
520  /* sanity check */
521  if(buffer != "skew-symmetric")
522  {
523  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
524  return false;
525  }
526  mmHeader.structure=skew_symmetric;
527  break;
528  default :
529  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
530  return false;
531  }
532  break;
533  default :
534  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
535  return false;
536  }
537  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
538  c=file.peek();
539  return true;
540 
541  }
542 
543  template<std::size_t brows, std::size_t bcols>
544  std::tuple<std::size_t, std::size_t, std::size_t>
545  calculateNNZ(std::size_t rows, std::size_t cols, std::size_t entries, const MMHeader& header)
546  {
547  std::size_t blockrows=rows/brows;
548  std::size_t blockcols=cols/bcols;
549  std::size_t blocksize=brows*bcols;
550  std::size_t blockentries=0;
551 
552  switch(header.structure)
553  {
554  case general :
555  blockentries = entries/blocksize; break;
556  case skew_symmetric :
557  blockentries = 2*entries/blocksize; break;
558  case symmetric :
559  blockentries = (2*entries-rows)/blocksize; break;
560  case hermitian :
561  blockentries = (2*entries-rows)/blocksize; break;
562  default :
563  throw Dune::NotImplemented();
564  }
565  return std::make_tuple(blockrows, blockcols, blockentries);
566  }
567 
568  /*
569  * @brief Storage class for the column index and the numeric value.
570  *
571  * \tparam T Either a NumericWrapper of the numeric type or PatternDummy
572  * for MatrixMarket pattern case.
573  */
574  template<typename T>
575  struct IndexData : public T
576  {
577  std::size_t index;
578  };
579 
580 
591  template<typename T>
593  {
595  operator T&()
596  {
597  return number;
598  }
599  };
600 
605  {};
606 
607  template<>
609  {};
610 
611  template<typename T>
612  std::istream& operator>>(std::istream& is, NumericWrapper<T>& num)
613  {
614  return is>>num.number;
615  }
616 
617  inline std::istream& operator>>(std::istream& is, NumericWrapper<PatternDummy>& num)
618  {
619  DUNE_UNUSED_PARAMETER(num);
620  return is;
621  }
622 
628  template<typename T>
629  bool operator<(const IndexData<T>& i1, const IndexData<T>& i2)
630  {
631  return i1.index<i2.index;
632  }
633 
639  template<typename T>
640  std::istream& operator>>(std::istream& is, IndexData<T>& data)
641  {
642  is>>data.index;
643  /* MatrixMarket indices are one based. Decrement for C++ */
644  --data.index;
645  return is>>data.number;
646  }
647 
653  template<typename T>
654  std::istream& operator>>(std::istream& is, IndexData<NumericWrapper<std::complex<T>>>& data)
655  {
656  is>>data.index;
657  /* MatrixMarket indices are one based. Decrement for C++ */
658  --data.index;
659  // real and imaginary part needs to be read separately as
660  // complex numbers are not provided in pair form. (x,y)
661  NumericWrapper<T> real, imag;
662  is>>real;
663  is>>imag;
664  data.number = {real.number, imag.number};
665  return is;
666  }
667 
674  template<typename D, int brows, int bcols>
676  {
682  template<typename T>
683  void operator()(const std::vector<std::set<IndexData<D> > >& rows,
684  BCRSMatrix<T>& matrix)
685  {
686  static_assert(IsNumber<T>::value && brows==1 && bcols==1, "Only scalar entries are expected here!");
687  for (auto iter=matrix.begin(); iter!= matrix.end(); ++iter)
688  {
689  auto brow=iter.index();
690  for (auto siter=rows[brow].begin(); siter != rows[brow].end(); ++siter)
691  (*iter)[siter->index] = siter->number;
692  }
693  }
694 
700  template<typename T>
701  void operator()(const std::vector<std::set<IndexData<D> > >& rows,
703  {
704  for (auto iter=matrix.begin(); iter!= matrix.end(); ++iter)
705  {
706  for (auto brow=iter.index()*brows,
707  browend=iter.index()*brows+brows;
708  brow<browend; ++brow)
709  {
710  for (auto siter=rows[brow].begin(), send=rows[brow].end();
711  siter != send; ++siter)
712  (*iter)[siter->index/bcols][brow%brows][siter->index%bcols]=siter->number;
713  }
714  }
715  }
716  };
717 
718  template<int brows, int bcols>
719  struct MatrixValuesSetter<PatternDummy,brows,bcols>
720  {
721  template<typename M>
722  void operator()(const std::vector<std::set<IndexData<PatternDummy> > >& rows,
723  M& matrix)
724  {}
725  };
726 
727  template<class T> struct is_complex : std::false_type {};
728  template<class T> struct is_complex<std::complex<T>> : std::true_type {};
729 
730  // wrapper for std::conj. Returns T if T is not complex.
731  template<class T>
732  std::enable_if_t<!is_complex<T>::value, T> conj(const T& r){
733  return r;
734  }
735 
736  template<class T>
737  std::enable_if_t<is_complex<T>::value, T> conj(const T& r){
738  return std::conj(r);
739  }
740 
741  template<typename M>
743  {};
744 
745  template<typename B, typename A>
747  {
748  enum {
749  rows = 1,
750  cols = 1
751  };
752  };
753 
754  template<typename B, int i, int j, typename A>
756  {
757  enum {
758  rows = i,
759  cols = j
760  };
761  };
762 
763  template<typename T, typename A, typename D>
765  std::istream& file, std::size_t entries,
766  const MMHeader& mmHeader, const D&)
767  {
769 
770  // Number of rows and columns of T, if it is a matrix (1x1 otherwise)
771  constexpr int brows = mm_multipliers<Matrix>::rows;
772  constexpr int bcols = mm_multipliers<Matrix>::cols;
773 
774  // First path
775  // store entries together with column index in a separate
776  // data structure
777  std::vector<std::set<IndexData<D> > > rows(matrix.N()*brows);
778 
779  auto readloop = [&] (auto symmetryFixup) {
780  for(std::size_t i = 0; i < entries; ++i) {
781  std::size_t row;
782  IndexData<D> data;
783  skipComments(file);
784  file>>row;
785  --row; // Index was 1 based.
786  assert(row/bcols<matrix.N());
787  file>>data;
788  assert(data.index/bcols<matrix.M());
789  rows[row].insert(data);
790  if(row!=data.index)
791  symmetryFixup(row, data);
792  }
793  };
794 
795  switch(mmHeader.structure)
796  {
797  case general:
798  readloop([](auto...){});
799  break;
800  case symmetric :
801  readloop([&](auto row, auto data) {
802  IndexData<D> data_sym(data);
803  data_sym.index = row;
804  rows[data.index].insert(data_sym);
805  });
806  break;
807  case skew_symmetric :
808  readloop([&](auto row, auto data) {
809  IndexData<D> data_sym;
810  data_sym.number = -data.number;
811  data_sym.index = row;
812  rows[data.index].insert(data_sym);
813  });
814  break;
815  case hermitian :
816  readloop([&](auto row, auto data) {
817  IndexData<D> data_sym;
818  data_sym.number = conj(data.number);
819  data_sym.index = row;
820  rows[data.index].insert(data_sym);
821  });
822  break;
823  default:
824  DUNE_THROW(Dune::NotImplemented,
825  "Only general, symmetric, skew-symmetric and hermitian is supported right now!");
826  }
827 
828  // Setup the matrix sparsity pattern
829  int nnz=0;
830  for(typename Matrix::CreateIterator iter=matrix.createbegin();
831  iter!= matrix.createend(); ++iter)
832  {
833  for(std::size_t brow=iter.index()*brows, browend=iter.index()*brows+brows;
834  brow<browend; ++brow)
835  {
836  typedef typename std::set<IndexData<D> >::const_iterator Siter;
837  for(Siter siter=rows[brow].begin(), send=rows[brow].end();
838  siter != send; ++siter, ++nnz)
839  iter.insert(siter->index/bcols);
840  }
841  }
842 
843  //Set the matrix values
844  matrix=0;
845 
847 
848  Setter(rows, matrix);
849  }
850  } // end namespace MatrixMarketImpl
851 
852  class MatrixMarketFormatError : public Dune::Exception
853  {};
854 
855 
856  inline void mm_read_header(std::size_t& rows, std::size_t& cols,
857  MatrixMarketImpl::MMHeader& header, std::istream& istr,
858  bool isVector)
859  {
860  using namespace MatrixMarketImpl;
861 
862  if(!readMatrixMarketBanner(istr, header)) {
863  std::cerr << "First line was not a correct Matrix Market banner. Using default:\n"
864  << "%%MatrixMarket matrix coordinate real general"<<std::endl;
865  // Go to the beginning of the file
866  istr.clear() ;
867  istr.seekg(0, std::ios::beg);
868  if(isVector)
869  header.type=array_type;
870  }
871 
872  skipComments(istr);
873 
874  if(lineFeed(istr))
875  throw MatrixMarketFormatError();
876 
877  istr >> rows;
878 
879  if(lineFeed(istr))
880  throw MatrixMarketFormatError();
881  istr >> cols;
882  }
883 
884  template<typename T, typename A>
886  std::size_t size,
887  std::istream& istr)
888  {
889  for (int i=0; size>0; ++i, --size)
890  istr>>vector[i];
891  }
892 
893  template<typename T, typename A, int entries>
894  void mm_read_vector_entries(Dune::BlockVector<Dune::FieldVector<T,entries>,A>& vector,
895  std::size_t size,
896  std::istream& istr)
897  {
898  for(int i=0; size>0; ++i, --size) {
899  T val;
900  istr>>val;
901  vector[i/entries][i%entries]=val;
902  }
903  }
904 
905 
912  template<typename T, typename A>
914  std::istream& istr)
915  {
916  using namespace MatrixMarketImpl;
917 
918  MMHeader header;
919  std::size_t rows, cols;
920  mm_read_header(rows,cols,header,istr, true);
921  if(cols!=1)
922  DUNE_THROW(MatrixMarketFormatError, "cols!=1, therefore this is no vector!");
923 
924  if(header.type!=array_type)
925  DUNE_THROW(MatrixMarketFormatError, "Vectors have to be stored in array format!");
926 
927 
928  Dune::Hybrid::ifElse(Dune::IsNumber<T>(),
929  [&](auto id){
930  vector.resize(rows);
931  },
932  [&](auto id){ // Assume that T is a FieldVector
933  T dummy;
934  auto blocksize = id(dummy).size();
935  std::size_t size=rows/blocksize;
936  if(size*blocksize!=rows)
937  DUNE_THROW(MatrixMarketFormatError, "Block size of vector is not correct!");
938 
939  vector.resize(size);
940  });
941 
942  istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
943  mm_read_vector_entries(vector, rows, istr);
944  }
945 
952  template<typename T, typename A>
954  std::istream& istr)
955  {
956  using namespace MatrixMarketImpl;
958 
959  MMHeader header;
960  if(!readMatrixMarketBanner(istr, header)) {
961  std::cerr << "First line was not a correct Matrix Market banner. Using default:\n"
962  << "%%MatrixMarket matrix coordinate real general"<<std::endl;
963  // Go to the beginning of the file
964  istr.clear() ;
965  istr.seekg(0, std::ios::beg);
966  }
967  skipComments(istr);
968 
969  std::size_t rows, cols, entries;
970 
971  if(lineFeed(istr))
972  throw MatrixMarketFormatError();
973 
974  istr >> rows;
975 
976  if(lineFeed(istr))
977  throw MatrixMarketFormatError();
978  istr >> cols;
979 
980  if(lineFeed(istr))
981  throw MatrixMarketFormatError();
982 
983  istr >>entries;
984 
985  std::size_t nnz, blockrows, blockcols;
986 
987  // Number of rows and columns of T, if it is a matrix (1x1 otherwise)
988  constexpr int brows = mm_multipliers<Matrix>::rows;
989  constexpr int bcols = mm_multipliers<Matrix>::cols;
990 
991  std::tie(blockrows, blockcols, nnz) = calculateNNZ<brows, bcols>(rows, cols, entries, header);
992 
993  istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
994 
995 
996  matrix.setSize(blockrows, blockcols);
998 
999  if(header.type==array_type)
1000  DUNE_THROW(Dune::NotImplemented, "Array format currently not supported for matrices!");
1001 
1002  readSparseEntries(matrix, istr, entries, header, NumericWrapper<typename Matrix::field_type>());
1003  }
1004 
1005  // Print a scalar entry
1006  template<typename B>
1007  void mm_print_entry(const B& entry,
1008  std::size_t rowidx,
1009  std::size_t colidx,
1010  std::ostream& ostr)
1011  {
1012  Hybrid::ifElse(IsNumber<B>(),
1013  [&](auto id) {
1014  ostr << rowidx << " " << colidx << " " << entry << std::endl;
1015  },
1016  [&](auto id) {
1017  for (auto row=id(entry).begin(); row != id(entry).end(); ++row, ++rowidx) {
1018  int coli=colidx;
1019  for (auto col = row->begin(); col != row->end(); ++col, ++coli)
1020  ostr<< rowidx<<" "<<coli<<" "<<*col<<std::endl;
1021  }
1022  });
1023  }
1024 
1025  // Write a vector entry
1026  template<typename V>
1027  void mm_print_vector_entry(const V& entry, std::ostream& ostr,
1028  const std::integral_constant<int,1>&)
1029  {
1030  ostr<<entry<<std::endl;
1031  }
1032 
1033  // Write a vector
1034  template<typename V>
1035  void mm_print_vector_entry(const V& vector, std::ostream& ostr,
1036  const std::integral_constant<int,0>&)
1037  {
1038  using namespace MatrixMarketImpl;
1039 
1040  // Is the entry a supported numeric type?
1042  typedef typename V::const_iterator VIter;
1043 
1044  for(VIter i=vector.begin(); i != vector.end(); ++i)
1045 
1046  mm_print_vector_entry(*i, ostr,
1047  std::integral_constant<int,isnumeric>());
1048  }
1049 
1050  template<typename T, typename A>
1051  std::size_t countEntries(const BlockVector<T,A>& vector)
1052  {
1053  return vector.size();
1054  }
1055 
1056  template<typename T, typename A, int i>
1057  std::size_t countEntries(const BlockVector<FieldVector<T,i>,A>& vector)
1058  {
1059  return vector.size()*i;
1060  }
1061 
1062  // Version for writing vectors.
1063  template<typename V>
1064  void writeMatrixMarket(const V& vector, std::ostream& ostr,
1065  const std::integral_constant<int,0>&)
1066  {
1067  using namespace MatrixMarketImpl;
1068 
1069  ostr<<countEntries(vector)<<" "<<1<<std::endl;
1071  mm_print_vector_entry(vector,ostr, std::integral_constant<int,isnumeric>());
1072  }
1073 
1074  // Versions for writing matrices
1075  template<typename M>
1076  void writeMatrixMarket(const M& matrix,
1077  std::ostream& ostr,
1078  const std::integral_constant<int,1>&)
1079  {
1080  ostr<<matrix.N()*MatrixMarketImpl::mm_multipliers<M>::rows<<" "
1082  <<countNonZeros(matrix)<<std::endl;
1083 
1084  typedef typename M::const_iterator riterator;
1085  typedef typename M::ConstColIterator citerator;
1086  for(riterator row=matrix.begin(); row != matrix.end(); ++row)
1087  for(citerator col = row->begin(); col != row->end(); ++col)
1088  // Matrix Market indexing start with 1!
1091  }
1092 
1093 
1097  template<typename M>
1098  void writeMatrixMarket(const M& matrix,
1099  std::ostream& ostr)
1100  {
1101  using namespace MatrixMarketImpl;
1102 
1103  // Write header information
1104  mm_header_printer<M>::print(ostr);
1105  mm_block_structure_header<M>::print(ostr,matrix);
1106  // Choose the correct function for matrix and vector
1107  writeMatrixMarket(matrix,ostr,std::integral_constant<int,IsMatrix<M>::value>());
1108  }
1109 
1110 
1121  template<typename M>
1122  void storeMatrixMarket(const M& matrix,
1123  std::string filename)
1124  {
1125  std::ofstream file(filename.c_str());
1126  file.setf(std::ios::scientific,std::ios::floatfield);
1127  writeMatrixMarket(matrix, file);
1128  file.close();
1129  }
1130 
1131 #if HAVE_MPI
1145  template<typename M, typename G, typename L>
1146  void storeMatrixMarket(const M& matrix,
1147  std::string filename,
1149  bool storeIndices=true)
1150  {
1151  // Get our rank
1152  int rank = comm.communicator().rank();
1153  // Write the local matrix
1154  std::ostringstream rfilename;
1155  rfilename<<filename <<"_"<<rank<<".mm";
1156  dverb<< rfilename.str()<<std::endl;
1157  std::ofstream file(rfilename.str().c_str());
1158  file.setf(std::ios::scientific,std::ios::floatfield);
1159  writeMatrixMarket(matrix, file);
1160  file.close();
1161 
1162  if(!storeIndices)
1163  return;
1164 
1165  // Write the global to local index mapping
1166  rfilename.str("");
1167  rfilename<<filename<<"_"<<rank<<".idx";
1168  file.open(rfilename.str().c_str());
1169  file.setf(std::ios::scientific,std::ios::floatfield);
1170  typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet IndexSet;
1171  typedef typename IndexSet::const_iterator Iterator;
1172  for(Iterator iter = comm.indexSet().begin();
1173  iter != comm.indexSet().end(); ++iter) {
1174  file << iter->global()<<" "<<(std::size_t)iter->local()<<" "
1175  <<(int)iter->local().attribute()<<" "<<(int)iter->local().isPublic()<<std::endl;
1176  }
1177  // Store neighbour information for efficient remote indices setup.
1178  file<<"neighbours:";
1179  const std::set<int>& neighbours=comm.remoteIndices().getNeighbours();
1180  typedef std::set<int>::const_iterator SIter;
1181  for(SIter neighbour=neighbours.begin(); neighbour != neighbours.end(); ++neighbour) {
1182  file<<" "<< *neighbour;
1183  }
1184  file.close();
1185  }
1186 
1201  template<typename M, typename G, typename L>
1202  void loadMatrixMarket(M& matrix,
1203  const std::string& filename,
1205  bool readIndices=true)
1206  {
1207  using namespace MatrixMarketImpl;
1208 
1210  typedef typename LocalIndex::Attribute Attribute;
1211  // Get our rank
1212  int rank = comm.communicator().rank();
1213  // load local matrix
1214  std::ostringstream rfilename;
1215  rfilename<<filename <<"_"<<rank<<".mm";
1216  std::ifstream file;
1217  file.open(rfilename.str().c_str(), std::ios::in);
1218  if(!file)
1219  DUNE_THROW(IOError, "Could not open file: " << rfilename.str().c_str());
1220  //if(!file.is_open())
1221  //
1222  readMatrixMarket(matrix,file);
1223  file.close();
1224 
1225  if(!readIndices)
1226  return;
1227 
1228  // read indices
1229  typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet IndexSet;
1230  IndexSet& pis=comm.pis;
1231  rfilename.str("");
1232  rfilename<<filename<<"_"<<rank<<".idx";
1233  file.open(rfilename.str().c_str());
1234  if(pis.size()!=0)
1235  DUNE_THROW(InvalidIndexSetState, "Index set is not empty!");
1236 
1237  pis.beginResize();
1238  while(!file.eof() && file.peek()!='n') {
1239  G g;
1240  file >>g;
1241  std::size_t l;
1242  file >>l;
1243  int c;
1244  file >>c;
1245  bool b;
1246  file >> b;
1247  pis.add(g,LocalIndex(l,Attribute(c),b));
1248  lineFeed(file);
1249  }
1250  pis.endResize();
1251  if(!file.eof()) {
1252  // read neighbours
1253  std::string s;
1254  file>>s;
1255  if(s!="neighbours:")
1256  DUNE_THROW(MatrixMarketFormatError, "was expecting the string: \"neighbours:\"");
1257  std::set<int> nb;
1258  while(!file.eof()) {
1259  int i;
1260  file >> i;
1261  nb.insert(i);
1262  }
1263  file.close();
1264  comm.ri.setNeighbours(nb);
1265  }
1266  comm.ri.template rebuild<false>();
1267  }
1268 
1269  #endif
1270 
1281  template<typename M>
1282  void loadMatrixMarket(M& matrix,
1283  const std::string& filename)
1284  {
1285  std::ifstream file;
1286  file.open(filename.c_str(), std::ios::in);
1287  if(!file)
1288  DUNE_THROW(IOError, "Could not open file: " << filename);
1289  readMatrixMarket(matrix,file);
1290  file.close();
1291  }
1292 
1294 }
1295 #endif
Implementation of the BCRSMatrix class.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Some handy generic functions for ISTL matrices.
Classes providing communication interfaces for overlapping Schwarz methods.
Col col
Definition: matrixmatrix.hh:349
auto countNonZeros(const M &matrix, typename std::enable_if_t< Dune::IsNumber< M >::value > *sfinae=nullptr)
Get the number of nonzero fields in the matrix.
Definition: matrixutils.hh:119
void readMatrixMarket(Dune::BlockVector< T, A > &vector, std::istream &istr)
Reads a BlockVector from a matrix market file.
Definition: matrixmarket.hh:913
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:1202
std::size_t countEntries(const BlockVector< T, A > &vector)
Definition: matrixmarket.hh:1051
void writeMatrixMarket(const V &vector, std::ostream &ostr, const std::integral_constant< int, 0 > &)
Definition: matrixmarket.hh:1064
void mm_print_vector_entry(const V &entry, std::ostream &ostr, const std::integral_constant< int, 1 > &)
Definition: matrixmarket.hh:1027
void storeMatrixMarket(const M &matrix, std::string filename)
Stores a parallel matrix/vector in matrix market format in a file.
Definition: matrixmarket.hh:1122
void mm_read_vector_entries(Dune::BlockVector< T, A > &vector, std::size_t size, std::istream &istr)
Definition: matrixmarket.hh:885
void mm_read_header(std::size_t &rows, std::size_t &cols, MatrixMarketImpl::MMHeader &header, std::istream &istr, bool isVector)
Definition: matrixmarket.hh:856
void mm_print_entry(const B &entry, std::size_t rowidx, std::size_t colidx, std::ostream &ostr)
Definition: matrixmarket.hh:1007
Definition: allocator.hh:7
std::istream & operator>>(std::istream &is, NumericWrapper< T > &num)
Definition: matrixmarket.hh:612
bool operator<(const IndexData< T > &i1, const IndexData< T > &i2)
LessThan operator.
Definition: matrixmarket.hh:629
LineType
Definition: matrixmarket.hh:294
@ DATA
Definition: matrixmarket.hh:294
@ MM_HEADER
Definition: matrixmarket.hh:294
@ MM_ISTLSTRUCT
Definition: matrixmarket.hh:294
bool readMatrixMarketBanner(std::istream &file, MMHeader &mmHeader)
Definition: matrixmarket.hh:351
std::enable_if_t< is_complex< T >::value, T > conj(const T &r)
Definition: matrixmarket.hh:737
std::tuple< std::size_t, std::size_t, std::size_t > calculateNNZ(std::size_t rows, std::size_t cols, std::size_t entries, const MMHeader &header)
Definition: matrixmarket.hh:545
void readSparseEntries(Dune::BCRSMatrix< T, A > &matrix, std::istream &file, std::size_t entries, const MMHeader &mmHeader, const D &)
Definition: matrixmarket.hh:764
MM_TYPE
Definition: matrixmarket.hh:297
@ array_type
Definition: matrixmarket.hh:297
@ coordinate_type
Definition: matrixmarket.hh:297
@ unknown_type
Definition: matrixmarket.hh:297
std::enable_if_t<!is_complex< T >::value, T > conj(const T &r)
Definition: matrixmarket.hh:732
void skipComments(std::istream &file)
Definition: matrixmarket.hh:337
bool lineFeed(std::istream &file)
Definition: matrixmarket.hh:313
@ MM_MAX_LINE_LENGTH
Definition: matrixmarket.hh:295
MM_STRUCTURE
Definition: matrixmarket.hh:301
@ skew_symmetric
Definition: matrixmarket.hh:301
@ general
Definition: matrixmarket.hh:301
@ hermitian
Definition: matrixmarket.hh:301
@ unknown_structure
Definition: matrixmarket.hh:301
@ symmetric
Definition: matrixmarket.hh:301
MM_CTYPE
Definition: matrixmarket.hh:299
@ unknown_ctype
Definition: matrixmarket.hh:299
@ pattern
Definition: matrixmarket.hh:299
@ complex_type
Definition: matrixmarket.hh:299
@ double_type
Definition: matrixmarket.hh:299
@ integer_type
Definition: matrixmarket.hh:299
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:426
Iterator begin()
Get iterator to first row.
Definition: bcrsmatrix.hh:634
Iterator end()
Get iterator to one beyond last row.
Definition: bcrsmatrix.hh:640
CreateIterator createend()
get create iterator pointing to one after the last block
Definition: bcrsmatrix.hh:1062
size_type M() const
number of columns (counted in blocks)
Definition: bcrsmatrix.hh:1937
CreateIterator createbegin()
get initial create iterator
Definition: bcrsmatrix.hh:1056
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1931
void setBuildMode(BuildMode bm)
Sets the build mode of the matrix.
Definition: bcrsmatrix.hh:792
void setSize(size_type rows, size_type columns, size_type nnz=0)
Set the size of the matrix.
Definition: bcrsmatrix.hh:820
A vector of blocks with memory management.
Definition: bvector.hh:403
void resize(size_type size)
Resize the vector.
Definition: bvector.hh:536
A generic dynamic dense matrix.
Definition: matrix.hh:558
Helper metaprogram to get the matrix market string representation of the numeric type.
Definition: matrixmarket.hh:74
@ is_numeric
Whether T is a supported numeric type.
Definition: matrixmarket.hh:79
static std::string str()
Definition: matrixmarket.hh:93
static std::string str()
Definition: matrixmarket.hh:109
static std::string str()
Definition: matrixmarket.hh:125
static std::string str()
Definition: matrixmarket.hh:141
static std::string str()
Definition: matrixmarket.hh:157
Meta program to write the correct Matrix Market header.
Definition: matrixmarket.hh:172
static void print(std::ostream &os)
Definition: matrixmarket.hh:177
static void print(std::ostream &os)
Definition: matrixmarket.hh:187
static void print(std::ostream &os)
Definition: matrixmarket.hh:197
static void print(std::ostream &os)
Definition: matrixmarket.hh:207
Metaprogram for writing the ISTL block structure header.
Definition: matrixmarket.hh:223
static void print(std::ostream &os, const M &)
Definition: matrixmarket.hh:231
BlockVector< T, A > M
Definition: matrixmarket.hh:228
static void print(std::ostream &os, const M &)
Definition: matrixmarket.hh:243
BCRSMatrix< T, A > M
Definition: matrixmarket.hh:253
static void print(std::ostream &os, const M &)
Definition: matrixmarket.hh:256
static void print(std::ostream &os, const M &)
Definition: matrixmarket.hh:268
static void print(std::ostream &os, const M &m)
Definition: matrixmarket.hh:281
FieldMatrix< T, i, j > M
Definition: matrixmarket.hh:279
static void print(std::ostream &os, const M &m)
Definition: matrixmarket.hh:290
FieldVector< T, i > M
Definition: matrixmarket.hh:288
Definition: matrixmarket.hh:304
MM_STRUCTURE structure
Definition: matrixmarket.hh:310
MM_TYPE type
Definition: matrixmarket.hh:308
MMHeader()
Definition: matrixmarket.hh:305
MM_CTYPE ctype
Definition: matrixmarket.hh:309
Definition: matrixmarket.hh:576
std::size_t index
Definition: matrixmarket.hh:577
a wrapper class of numeric values.
Definition: matrixmarket.hh:593
T number
Definition: matrixmarket.hh:594
Utility class for marking the pattern type of the MatrixMarket matrices.
Definition: matrixmarket.hh:605
Functor to the data values of the matrix.
Definition: matrixmarket.hh:676
void operator()(const std::vector< std::set< IndexData< D > > > &rows, BCRSMatrix< T > &matrix)
Sets the matrix values.
Definition: matrixmarket.hh:683
void operator()(const std::vector< std::set< IndexData< D > > > &rows, BCRSMatrix< FieldMatrix< T, brows, bcols > > &matrix)
Sets the matrix values.
Definition: matrixmarket.hh:701
void operator()(const std::vector< std::set< IndexData< PatternDummy > > > &rows, M &matrix)
Definition: matrixmarket.hh:722
Definition: matrixmarket.hh:727
Definition: matrixmarket.hh:743
Definition: matrixmarket.hh:853
Definition: matrixutils.hh:26
Test whether a type is an ISTL Matrix.
Definition: matrixutils.hh:502
A class setting up standard communication for a two-valued attribute set with owner/overlap/copy sema...
Definition: owneroverlapcopy.hh:172
const CollectiveCommunication< MPI_Comm > & communicator() const
Definition: owneroverlapcopy.hh:311
const ParallelIndexSet & indexSet() const
Get the underlying parallel index set.
Definition: owneroverlapcopy.hh:472
const RemoteIndices & remoteIndices() const
Get the underlying remote indices.
Definition: owneroverlapcopy.hh:481
Dune::ParallelIndexSet< GlobalIdType, LI, 512 > ParallelIndexSet
The type of the parallel index set.
Definition: owneroverlapcopy.hh:459