dune-istl  2.8.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 <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/hybridutilities.hh>
26 #include <dune/common/stdstreams.hh>
27 #include <dune/common/simd/simd.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<Simd::Scalar<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<Simd::Scalar<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, [[maybe_unused]] NumericWrapper<PatternDummy>& num)
618  {
619  return is;
620  }
621 
627  template<typename T>
628  bool operator<(const IndexData<T>& i1, const IndexData<T>& i2)
629  {
630  return i1.index<i2.index;
631  }
632 
638  template<typename T>
639  std::istream& operator>>(std::istream& is, IndexData<T>& data)
640  {
641  is>>data.index;
642  /* MatrixMarket indices are one based. Decrement for C++ */
643  --data.index;
644  return is>>data.number;
645  }
646 
652  template<typename T>
653  std::istream& operator>>(std::istream& is, IndexData<NumericWrapper<std::complex<T>>>& data)
654  {
655  is>>data.index;
656  /* MatrixMarket indices are one based. Decrement for C++ */
657  --data.index;
658  // real and imaginary part needs to be read separately as
659  // complex numbers are not provided in pair form. (x,y)
660  NumericWrapper<T> real, imag;
661  is>>real;
662  is>>imag;
663  data.number = {real.number, imag.number};
664  return is;
665  }
666 
673  template<typename D, int brows, int bcols>
675  {
681  template<typename T>
682  void operator()(const std::vector<std::set<IndexData<D> > >& rows,
683  BCRSMatrix<T>& matrix)
684  {
685  static_assert(IsNumber<T>::value && brows==1 && bcols==1, "Only scalar entries are expected here!");
686  for (auto iter=matrix.begin(); iter!= matrix.end(); ++iter)
687  {
688  auto brow=iter.index();
689  for (auto siter=rows[brow].begin(); siter != rows[brow].end(); ++siter)
690  (*iter)[siter->index] = siter->number;
691  }
692  }
693 
699  template<typename T>
700  void operator()(const std::vector<std::set<IndexData<D> > >& rows,
702  {
703  for (auto iter=matrix.begin(); iter!= matrix.end(); ++iter)
704  {
705  for (auto brow=iter.index()*brows,
706  browend=iter.index()*brows+brows;
707  brow<browend; ++brow)
708  {
709  for (auto siter=rows[brow].begin(), send=rows[brow].end();
710  siter != send; ++siter)
711  (*iter)[siter->index/bcols][brow%brows][siter->index%bcols]=siter->number;
712  }
713  }
714  }
715  };
716 
717  template<int brows, int bcols>
718  struct MatrixValuesSetter<PatternDummy,brows,bcols>
719  {
720  template<typename M>
721  void operator()(const std::vector<std::set<IndexData<PatternDummy> > >& rows,
722  M& matrix)
723  {}
724  };
725 
726  template<class T> struct is_complex : std::false_type {};
727  template<class T> struct is_complex<std::complex<T>> : std::true_type {};
728 
729  // wrapper for std::conj. Returns T if T is not complex.
730  template<class T>
731  std::enable_if_t<!is_complex<T>::value, T> conj(const T& r){
732  return r;
733  }
734 
735  template<class T>
736  std::enable_if_t<is_complex<T>::value, T> conj(const T& r){
737  return std::conj(r);
738  }
739 
740  template<typename M>
742  {};
743 
744  template<typename B, typename A>
746  {
747  enum {
748  rows = 1,
749  cols = 1
750  };
751  };
752 
753  template<typename B, int i, int j, typename A>
755  {
756  enum {
757  rows = i,
758  cols = j
759  };
760  };
761 
762  template<typename T, typename A, typename D>
764  std::istream& file, std::size_t entries,
765  const MMHeader& mmHeader, const D&)
766  {
768 
769  // Number of rows and columns of T, if it is a matrix (1x1 otherwise)
770  constexpr int brows = mm_multipliers<Matrix>::rows;
771  constexpr int bcols = mm_multipliers<Matrix>::cols;
772 
773  // First path
774  // store entries together with column index in a separate
775  // data structure
776  std::vector<std::set<IndexData<D> > > rows(matrix.N()*brows);
777 
778  auto readloop = [&] (auto symmetryFixup) {
779  for(std::size_t i = 0; i < entries; ++i) {
780  std::size_t row;
781  IndexData<D> data;
782  skipComments(file);
783  file>>row;
784  --row; // Index was 1 based.
785  assert(row/bcols<matrix.N());
786  file>>data;
787  assert(data.index/bcols<matrix.M());
788  rows[row].insert(data);
789  if(row!=data.index)
790  symmetryFixup(row, data);
791  }
792  };
793 
794  switch(mmHeader.structure)
795  {
796  case general:
797  readloop([](auto...){});
798  break;
799  case symmetric :
800  readloop([&](auto row, auto data) {
801  IndexData<D> data_sym(data);
802  data_sym.index = row;
803  rows[data.index].insert(data_sym);
804  });
805  break;
806  case skew_symmetric :
807  readloop([&](auto row, auto data) {
808  IndexData<D> data_sym;
809  data_sym.number = -data.number;
810  data_sym.index = row;
811  rows[data.index].insert(data_sym);
812  });
813  break;
814  case hermitian :
815  readloop([&](auto row, auto data) {
816  IndexData<D> data_sym;
817  data_sym.number = conj(data.number);
818  data_sym.index = row;
819  rows[data.index].insert(data_sym);
820  });
821  break;
822  default:
823  DUNE_THROW(Dune::NotImplemented,
824  "Only general, symmetric, skew-symmetric and hermitian is supported right now!");
825  }
826 
827  // Setup the matrix sparsity pattern
828  int nnz=0;
829  for(typename Matrix::CreateIterator iter=matrix.createbegin();
830  iter!= matrix.createend(); ++iter)
831  {
832  for(std::size_t brow=iter.index()*brows, browend=iter.index()*brows+brows;
833  brow<browend; ++brow)
834  {
835  typedef typename std::set<IndexData<D> >::const_iterator Siter;
836  for(Siter siter=rows[brow].begin(), send=rows[brow].end();
837  siter != send; ++siter, ++nnz)
838  iter.insert(siter->index/bcols);
839  }
840  }
841 
842  //Set the matrix values
843  matrix=0;
844 
846 
847  Setter(rows, matrix);
848  }
849  } // end namespace MatrixMarketImpl
850 
851  class MatrixMarketFormatError : public Dune::Exception
852  {};
853 
854 
855  inline void mm_read_header(std::size_t& rows, std::size_t& cols,
856  MatrixMarketImpl::MMHeader& header, std::istream& istr,
857  bool isVector)
858  {
859  using namespace MatrixMarketImpl;
860 
861  if(!readMatrixMarketBanner(istr, header)) {
862  std::cerr << "First line was not a correct Matrix Market banner. Using default:\n"
863  << "%%MatrixMarket matrix coordinate real general"<<std::endl;
864  // Go to the beginning of the file
865  istr.clear() ;
866  istr.seekg(0, std::ios::beg);
867  if(isVector)
868  header.type=array_type;
869  }
870 
871  skipComments(istr);
872 
873  if(lineFeed(istr))
874  throw MatrixMarketFormatError();
875 
876  istr >> rows;
877 
878  if(lineFeed(istr))
879  throw MatrixMarketFormatError();
880  istr >> cols;
881  }
882 
883  template<typename T, typename A>
885  std::size_t size,
886  std::istream& istr,
887  size_t lane)
888  {
889  for (int i=0; size>0; ++i, --size)
890  istr>>Simd::lane(lane,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  size_t lane)
898  {
899  for(int i=0; size>0; ++i, --size) {
900  Simd::Scalar<T> val;
901  istr>>val;
902  Simd::lane(lane, vector[i/entries][i%entries])=val;
903  }
904  }
905 
906 
913  template<typename T, typename A>
915  std::istream& istr)
916  {
917  typedef typename Dune::BlockVector<T,A>::field_type field_type;
918  using namespace MatrixMarketImpl;
919 
920  MMHeader header;
921  std::size_t rows, cols;
922  mm_read_header(rows,cols,header,istr, true);
923  if(cols!=Simd::lanes<field_type>()) {
924  if(Simd::lanes<field_type>() == 1)
925  DUNE_THROW(MatrixMarketFormatError, "cols!=1, therefore this is no vector!");
926  else
927  DUNE_THROW(MatrixMarketFormatError, "cols does not match the number of lanes in the field_type!");
928  }
929 
930  if(header.type!=array_type)
931  DUNE_THROW(MatrixMarketFormatError, "Vectors have to be stored in array format!");
932 
933 
934  if constexpr (Dune::IsNumber<T>())
935  vector.resize(rows);
936  else
937  {
938  T dummy;
939  auto blocksize = dummy.size();
940  std::size_t size=rows/blocksize;
941  if(size*blocksize!=rows)
942  DUNE_THROW(MatrixMarketFormatError, "Block size of vector is not correct!");
943 
944  vector.resize(size);
945  }
946 
947  istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
948  for(size_t l=0;l<Simd::lanes<field_type>();++l){
949  mm_read_vector_entries(vector, rows, istr, l);
950  }
951  }
952 
959  template<typename T, typename A>
961  std::istream& istr)
962  {
963  using namespace MatrixMarketImpl;
965 
966  MMHeader header;
967  if(!readMatrixMarketBanner(istr, header)) {
968  std::cerr << "First line was not a correct Matrix Market banner. Using default:\n"
969  << "%%MatrixMarket matrix coordinate real general"<<std::endl;
970  // Go to the beginning of the file
971  istr.clear() ;
972  istr.seekg(0, std::ios::beg);
973  }
974  skipComments(istr);
975 
976  std::size_t rows, cols, entries;
977 
978  if(lineFeed(istr))
979  throw MatrixMarketFormatError();
980 
981  istr >> rows;
982 
983  if(lineFeed(istr))
984  throw MatrixMarketFormatError();
985  istr >> cols;
986 
987  if(lineFeed(istr))
988  throw MatrixMarketFormatError();
989 
990  istr >>entries;
991 
992  std::size_t nnz, blockrows, blockcols;
993 
994  // Number of rows and columns of T, if it is a matrix (1x1 otherwise)
995  constexpr int brows = mm_multipliers<Matrix>::rows;
996  constexpr int bcols = mm_multipliers<Matrix>::cols;
997 
998  std::tie(blockrows, blockcols, nnz) = calculateNNZ<brows, bcols>(rows, cols, entries, header);
999 
1000  istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
1001 
1002 
1003  matrix.setSize(blockrows, blockcols);
1005 
1006  if(header.type==array_type)
1007  DUNE_THROW(Dune::NotImplemented, "Array format currently not supported for matrices!");
1008 
1009  readSparseEntries(matrix, istr, entries, header, NumericWrapper<typename Matrix::field_type>());
1010  }
1011 
1012  // Print a scalar entry
1013  template<typename B>
1014  void mm_print_entry(const B& entry,
1015  std::size_t rowidx,
1016  std::size_t colidx,
1017  std::ostream& ostr)
1018  {
1019  if constexpr (IsNumber<B>())
1020  ostr << rowidx << " " << colidx << " " << entry << std::endl;
1021  else
1022  {
1023  for (auto row=entry.begin(); row != entry.end(); ++row, ++rowidx) {
1024  int coli=colidx;
1025  for (auto col = row->begin(); col != row->end(); ++col, ++coli)
1026  ostr<< rowidx<<" "<<coli<<" "<<*col<<std::endl;
1027  }
1028  }
1029  }
1030 
1031  // Write a vector entry
1032  template<typename V>
1033  void mm_print_vector_entry(const V& entry, std::ostream& ostr,
1034  const std::integral_constant<int,1>&,
1035  size_t lane)
1036  {
1037  ostr<<Simd::lane(lane,entry)<<std::endl;
1038  }
1039 
1040  // Write a vector
1041  template<typename V>
1042  void mm_print_vector_entry(const V& vector, std::ostream& ostr,
1043  const std::integral_constant<int,0>&,
1044  size_t lane)
1045  {
1046  using namespace MatrixMarketImpl;
1047 
1048  // Is the entry a supported numeric type?
1049  const int isnumeric = mm_numeric_type<Simd::Scalar<typename V::block_type>>::is_numeric;
1050  typedef typename V::const_iterator VIter;
1051 
1052  for(VIter i=vector.begin(); i != vector.end(); ++i)
1053 
1054  mm_print_vector_entry(*i, ostr,
1055  std::integral_constant<int,isnumeric>(),
1056  lane);
1057  }
1058 
1059  template<typename T, typename A>
1060  std::size_t countEntries(const BlockVector<T,A>& vector)
1061  {
1062  return vector.size();
1063  }
1064 
1065  template<typename T, typename A, int i>
1066  std::size_t countEntries(const BlockVector<FieldVector<T,i>,A>& vector)
1067  {
1068  return vector.size()*i;
1069  }
1070 
1071  // Version for writing vectors.
1072  template<typename V>
1073  void writeMatrixMarket(const V& vector, std::ostream& ostr,
1074  const std::integral_constant<int,0>&)
1075  {
1076  using namespace MatrixMarketImpl;
1077  typedef typename V::field_type field_type;
1078 
1079  ostr<<countEntries(vector)<<" "<<Simd::lanes<field_type>()<<std::endl;
1080  const int isnumeric = mm_numeric_type<Simd::Scalar<V>>::is_numeric;
1081  for(size_t l=0;l<Simd::lanes<field_type>(); ++l){
1082  mm_print_vector_entry(vector,ostr, std::integral_constant<int,isnumeric>(), l);
1083  }
1084  }
1085 
1086  // Versions for writing matrices
1087  template<typename M>
1088  void writeMatrixMarket(const M& matrix,
1089  std::ostream& ostr,
1090  const std::integral_constant<int,1>&)
1091  {
1092  ostr<<matrix.N()*MatrixMarketImpl::mm_multipliers<M>::rows<<" "
1094  <<countNonZeros(matrix)<<std::endl;
1095 
1096  typedef typename M::const_iterator riterator;
1097  typedef typename M::ConstColIterator citerator;
1098  for(riterator row=matrix.begin(); row != matrix.end(); ++row)
1099  for(citerator col = row->begin(); col != row->end(); ++col)
1100  // Matrix Market indexing start with 1!
1103  }
1104 
1105 
1109  template<typename M>
1110  void writeMatrixMarket(const M& matrix,
1111  std::ostream& ostr)
1112  {
1113  using namespace MatrixMarketImpl;
1114 
1115  // Write header information
1116  mm_header_printer<M>::print(ostr);
1117  mm_block_structure_header<M>::print(ostr,matrix);
1118  // Choose the correct function for matrix and vector
1119  writeMatrixMarket(matrix,ostr,std::integral_constant<int,IsMatrix<M>::value>());
1120  }
1121 
1122  static const int default_precision = -1;
1134  template<typename M>
1135  void storeMatrixMarket(const M& matrix,
1136  std::string filename,
1137  int prec=default_precision)
1138  {
1139  std::ofstream file(filename.c_str());
1140  file.setf(std::ios::scientific,std::ios::floatfield);
1141  if(prec>0)
1142  file.precision(prec);
1143  writeMatrixMarket(matrix, file);
1144  file.close();
1145  }
1146 
1147 #if HAVE_MPI
1162  template<typename M, typename G, typename L>
1163  void storeMatrixMarket(const M& matrix,
1164  std::string filename,
1166  bool storeIndices=true,
1167  int prec=default_precision)
1168  {
1169  // Get our rank
1170  int rank = comm.communicator().rank();
1171  // Write the local matrix
1172  std::ostringstream rfilename;
1173  rfilename<<filename <<"_"<<rank<<".mm";
1174  dverb<< rfilename.str()<<std::endl;
1175  std::ofstream file(rfilename.str().c_str());
1176  file.setf(std::ios::scientific,std::ios::floatfield);
1177  if(prec>0)
1178  file.precision(prec);
1179  writeMatrixMarket(matrix, file);
1180  file.close();
1181 
1182  if(!storeIndices)
1183  return;
1184 
1185  // Write the global to local index mapping
1186  rfilename.str("");
1187  rfilename<<filename<<"_"<<rank<<".idx";
1188  file.open(rfilename.str().c_str());
1189  file.setf(std::ios::scientific,std::ios::floatfield);
1190  typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet IndexSet;
1191  typedef typename IndexSet::const_iterator Iterator;
1192  for(Iterator iter = comm.indexSet().begin();
1193  iter != comm.indexSet().end(); ++iter) {
1194  file << iter->global()<<" "<<(std::size_t)iter->local()<<" "
1195  <<(int)iter->local().attribute()<<" "<<(int)iter->local().isPublic()<<std::endl;
1196  }
1197  // Store neighbour information for efficient remote indices setup.
1198  file<<"neighbours:";
1199  const std::set<int>& neighbours=comm.remoteIndices().getNeighbours();
1200  typedef std::set<int>::const_iterator SIter;
1201  for(SIter neighbour=neighbours.begin(); neighbour != neighbours.end(); ++neighbour) {
1202  file<<" "<< *neighbour;
1203  }
1204  file.close();
1205  }
1206 
1221  template<typename M, typename G, typename L>
1222  void loadMatrixMarket(M& matrix,
1223  const std::string& filename,
1225  bool readIndices=true)
1226  {
1227  using namespace MatrixMarketImpl;
1228 
1230  typedef typename LocalIndexT::Attribute Attribute;
1231  // Get our rank
1232  int rank = comm.communicator().rank();
1233  // load local matrix
1234  std::ostringstream rfilename;
1235  rfilename<<filename <<"_"<<rank<<".mm";
1236  std::ifstream file;
1237  file.open(rfilename.str().c_str(), std::ios::in);
1238  if(!file)
1239  DUNE_THROW(IOError, "Could not open file: " << rfilename.str().c_str());
1240  //if(!file.is_open())
1241  //
1242  readMatrixMarket(matrix,file);
1243  file.close();
1244 
1245  if(!readIndices)
1246  return;
1247 
1248  // read indices
1249  typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet IndexSet;
1250  IndexSet& pis=comm.pis;
1251  rfilename.str("");
1252  rfilename<<filename<<"_"<<rank<<".idx";
1253  file.open(rfilename.str().c_str());
1254  if(pis.size()!=0)
1255  DUNE_THROW(InvalidIndexSetState, "Index set is not empty!");
1256 
1257  pis.beginResize();
1258  while(!file.eof() && file.peek()!='n') {
1259  G g;
1260  file >>g;
1261  std::size_t l;
1262  file >>l;
1263  int c;
1264  file >>c;
1265  bool b;
1266  file >> b;
1267  pis.add(g,LocalIndexT(l,Attribute(c),b));
1268  lineFeed(file);
1269  }
1270  pis.endResize();
1271  if(!file.eof()) {
1272  // read neighbours
1273  std::string s;
1274  file>>s;
1275  if(s!="neighbours:")
1276  DUNE_THROW(MatrixMarketFormatError, "was expecting the string: \"neighbours:\"");
1277  std::set<int> nb;
1278  while(!file.eof()) {
1279  int i;
1280  file >> i;
1281  nb.insert(i);
1282  }
1283  file.close();
1284  comm.ri.setNeighbours(nb);
1285  }
1286  comm.ri.template rebuild<false>();
1287  }
1288 
1289  #endif
1290 
1301  template<typename M>
1302  void loadMatrixMarket(M& matrix,
1303  const std::string& filename)
1304  {
1305  std::ifstream file;
1306  file.open(filename.c_str(), std::ios::in);
1307  if(!file)
1308  DUNE_THROW(IOError, "Could not open file: " << filename);
1309  readMatrixMarket(matrix,file);
1310  file.close();
1311  }
1312 
1314 }
1315 #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.
auto countNonZeros(const M &, [[maybe_unused]] typename std::enable_if_t< Dune::IsNumber< M >::value > *sfinae=nullptr)
Get the number of nonzero fields in the matrix.
Definition: matrixutils.hh:117
Col col
Definition: matrixmatrix.hh:349
void readMatrixMarket(Dune::BlockVector< T, A > &vector, std::istream &istr)
Reads a BlockVector from a matrix market file.
Definition: matrixmarket.hh:914
void storeMatrixMarket(const M &matrix, std::string filename, int prec=default_precision)
Stores a parallel matrix/vector in matrix market format in a file.
Definition: matrixmarket.hh:1135
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:1222
std::size_t countEntries(const BlockVector< T, A > &vector)
Definition: matrixmarket.hh:1060
void writeMatrixMarket(const V &vector, std::ostream &ostr, const std::integral_constant< int, 0 > &)
Definition: matrixmarket.hh:1073
void mm_print_vector_entry(const V &entry, std::ostream &ostr, const std::integral_constant< int, 1 > &, size_t lane)
Definition: matrixmarket.hh:1033
static const int default_precision
Definition: matrixmarket.hh:1122
void mm_read_vector_entries(Dune::BlockVector< T, A > &vector, std::size_t size, std::istream &istr, size_t lane)
Definition: matrixmarket.hh:884
void mm_read_header(std::size_t &rows, std::size_t &cols, MatrixMarketImpl::MMHeader &header, std::istream &istr, bool isVector)
Definition: matrixmarket.hh:855
void mm_print_entry(const B &entry, std::size_t rowidx, std::size_t colidx, std::ostream &ostr)
Definition: matrixmarket.hh:1014
Definition: allocator.hh:9
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:628
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:736
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:763
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:731
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:464
Iterator begin()
Get iterator to first row.
Definition: bcrsmatrix.hh:673
Iterator end()
Get iterator to one beyond last row.
Definition: bcrsmatrix.hh:679
CreateIterator createend()
get create iterator pointing to one after the last block
Definition: bcrsmatrix.hh:1101
size_type M() const
number of columns (counted in blocks)
Definition: bcrsmatrix.hh:1976
CreateIterator createbegin()
get initial create iterator
Definition: bcrsmatrix.hh:1095
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1970
void setBuildMode(BuildMode bm)
Sets the build mode of the matrix.
Definition: bcrsmatrix.hh:831
void setSize(size_type rows, size_type columns, size_type nnz=0)
Set the size of the matrix.
Definition: bcrsmatrix.hh:859
A vector of blocks with memory management.
Definition: bvector.hh:393
void resize(size_type size)
Resize the vector.
Definition: bvector.hh:501
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: bvector.hh:399
A generic dynamic dense matrix.
Definition: matrix.hh:559
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:675
void operator()(const std::vector< std::set< IndexData< D > > > &rows, BCRSMatrix< T > &matrix)
Sets the matrix values.
Definition: matrixmarket.hh:682
void operator()(const std::vector< std::set< IndexData< D > > > &rows, BCRSMatrix< FieldMatrix< T, brows, bcols > > &matrix)
Sets the matrix values.
Definition: matrixmarket.hh:700
void operator()(const std::vector< std::set< IndexData< PatternDummy > > > &rows, M &matrix)
Definition: matrixmarket.hh:721
Definition: matrixmarket.hh:726
Definition: matrixmarket.hh:742
Definition: matrixmarket.hh:852
Definition: matrixutils.hh:25
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:297
const ParallelIndexSet & indexSet() const
Get the underlying parallel index set.
Definition: owneroverlapcopy.hh:460
const RemoteIndices & remoteIndices() const
Get the underlying remote indices.
Definition: owneroverlapcopy.hh:469
Dune::ParallelIndexSet< GlobalIdType, LI, 512 > ParallelIndexSet
The type of the parallel index set.
Definition: owneroverlapcopy.hh:447