dune-istl  2.9.0
matrixmarket.hh
Go to the documentation of this file.
1 // SPDX-FileCopyrightText: Copyright (C) DUNE Project contributors, see file LICENSE.md in module root
2 // SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4 // vi: set et ts=4 sw=2 sts=2:
5 #ifndef DUNE_ISTL_MATRIXMARKET_HH
6 #define DUNE_ISTL_MATRIXMARKET_HH
7 
8 #include <algorithm>
9 #include <complex>
10 #include <cstddef>
11 #include <fstream>
12 #include <ios>
13 #include <iostream>
14 #include <istream>
15 #include <limits>
16 #include <ostream>
17 #include <set>
18 #include <sstream>
19 #include <string>
20 #include <tuple>
21 #include <type_traits>
22 #include <vector>
23 
24 #include <dune/common/exceptions.hh>
25 #include <dune/common/fmatrix.hh>
26 #include <dune/common/fvector.hh>
27 #include <dune/common/hybridutilities.hh>
28 #include <dune/common/stdstreams.hh>
29 #include <dune/common/simd/simd.hh>
30 
31 #include <dune/istl/bcrsmatrix.hh>
32 #include <dune/istl/bvector.hh>
33 #include <dune/istl/matrixutils.hh> // countNonZeros()
35 
36 namespace Dune
37 {
38 
64  namespace MatrixMarketImpl
65  {
75  template<class T>
76  struct mm_numeric_type {
77  enum {
81  is_numeric=false
82  };
83  };
84 
85  template<>
86  struct mm_numeric_type<int>
87  {
88  enum {
92  is_numeric=true
93  };
94 
95  static std::string str()
96  {
97  return "integer";
98  }
99  };
100 
101  template<>
102  struct mm_numeric_type<double>
103  {
104  enum {
108  is_numeric=true
109  };
110 
111  static std::string str()
112  {
113  return "real";
114  }
115  };
116 
117  template<>
118  struct mm_numeric_type<float>
119  {
120  enum {
124  is_numeric=true
125  };
126 
127  static std::string str()
128  {
129  return "real";
130  }
131  };
132 
133  template<>
134  struct mm_numeric_type<std::complex<double> >
135  {
136  enum {
140  is_numeric=true
141  };
142 
143  static std::string str()
144  {
145  return "complex";
146  }
147  };
148 
149  template<>
150  struct mm_numeric_type<std::complex<float> >
151  {
152  enum {
156  is_numeric=true
157  };
158 
159  static std::string str()
160  {
161  return "complex";
162  }
163  };
164 
173  template<class M>
175 
176  template<typename T, typename A>
178  {
179  static void print(std::ostream& os)
180  {
181  os<<"%%MatrixMarket matrix coordinate ";
182  os<<mm_numeric_type<Simd::Scalar<typename Imp::BlockTraits<T>::field_type>>::str()<<" general"<<std::endl;
183  }
184  };
185 
186  template<typename B, typename A>
188  {
189  static void print(std::ostream& os)
190  {
191  os<<"%%MatrixMarket matrix array ";
192  os<<mm_numeric_type<Simd::Scalar<typename Imp::BlockTraits<B>::field_type>>::str()<<" general"<<std::endl;
193  }
194  };
195 
196  template<typename T, int j>
197  struct mm_header_printer<FieldVector<T,j> >
198  {
199  static void print(std::ostream& os)
200  {
201  os<<"%%MatrixMarket matrix array ";
202  os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
203  }
204  };
205 
206  template<typename T, int i, int j>
208  {
209  static void print(std::ostream& os)
210  {
211  os<<"%%MatrixMarket matrix array ";
212  os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
213  }
214  };
215 
224  template<class M>
226 
227  template<typename T, typename A>
229  {
231  static_assert(IsNumber<T>::value, "Only scalar entries are expected here!");
232 
233  static void print(std::ostream& os, const M&)
234  {
235  os<<"% ISTL_STRUCT blocked ";
236  os<<"1 1"<<std::endl;
237  }
238  };
239 
240  template<typename T, typename A, int i>
241  struct mm_block_structure_header<BlockVector<FieldVector<T,i>,A> >
242  {
244 
245  static void print(std::ostream& os, const M&)
246  {
247  os<<"% ISTL_STRUCT blocked ";
248  os<<i<<" "<<1<<std::endl;
249  }
250  };
251 
252  template<typename T, typename A>
254  {
256  static_assert(IsNumber<T>::value, "Only scalar entries are expected here!");
257 
258  static void print(std::ostream& os, const M&)
259  {
260  os<<"% ISTL_STRUCT blocked ";
261  os<<"1 1"<<std::endl;
262  }
263  };
264 
265  template<typename T, typename A, int i, int j>
267  {
269 
270  static void print(std::ostream& os, const M&)
271  {
272  os<<"% ISTL_STRUCT blocked ";
273  os<<i<<" "<<j<<std::endl;
274  }
275  };
276 
277 
278  template<typename T, int i, int j>
280  {
282 
283  static void print(std::ostream& os, const M& m)
284  {}
285  };
286 
287  template<typename T, int i>
288  struct mm_block_structure_header<FieldVector<T,i> >
289  {
290  typedef FieldVector<T,i> M;
291 
292  static void print(std::ostream& os, const M& m)
293  {}
294  };
295 
297  enum { MM_MAX_LINE_LENGTH=1025 };
298 
300 
302 
304 
305  struct MMHeader
306  {
309  {}
313  };
314 
315  inline bool lineFeed(std::istream& file)
316  {
317  char c;
318  if(!file.eof())
319  c=file.peek();
320  else
321  return false;
322  // ignore whitespace
323  while(c==' ')
324  {
325  file.get();
326  if(file.eof())
327  return false;
328  c=file.peek();
329  }
330 
331  if(c=='\n') {
332  /* eat the line feed */
333  file.get();
334  return true;
335  }
336  return false;
337  }
338 
339  inline void skipComments(std::istream& file)
340  {
341  lineFeed(file);
342  char c=file.peek();
343  // ignore comment lines
344  while(c=='%')
345  {
346  /* discard the rest of the line */
347  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
348  c=file.peek();
349  }
350  }
351 
352 
353  inline bool readMatrixMarketBanner(std::istream& file, MMHeader& mmHeader)
354  {
355  std::string buffer;
356  char c;
357  file >> buffer;
358  c=buffer[0];
359  mmHeader=MMHeader();
360  if(c!='%')
361  return false;
362  dverb<<buffer<<std::endl;
363  /* read the banner */
364  if(buffer!="%%MatrixMarket") {
365  /* discard the rest of the line */
366  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
367  return false;
368  }
369 
370  if(lineFeed(file))
371  /* premature end of line */
372  return false;
373 
374  /* read the matrix_type */
375  file >> buffer;
376 
377  if(buffer != "matrix")
378  {
379  /* discard the rest of the line */
380  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
381  return false;
382  }
383 
384  if(lineFeed(file))
385  /* premature end of line */
386  return false;
387 
388  /* The type of the matrix */
389  file >> buffer;
390 
391  if(buffer.empty())
392  return false;
393 
394  std::transform(buffer.begin(), buffer.end(), buffer.begin(),
395  ::tolower);
396 
397  switch(buffer[0])
398  {
399  case 'a' :
400  /* sanity check */
401  if(buffer != "array")
402  {
403  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
404  return false;
405  }
406  mmHeader.type=array_type;
407  break;
408  case 'c' :
409  /* sanity check */
410  if(buffer != "coordinate")
411  {
412  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
413  return false;
414  }
415  mmHeader.type=coordinate_type;
416  break;
417  default :
418  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
419  return false;
420  }
421 
422  if(lineFeed(file))
423  /* premature end of line */
424  return false;
425 
426  /* The numeric type used. */
427  file >> buffer;
428 
429  if(buffer.empty())
430  return false;
431 
432  std::transform(buffer.begin(), buffer.end(), buffer.begin(),
433  ::tolower);
434  switch(buffer[0])
435  {
436  case 'i' :
437  /* sanity check */
438  if(buffer != "integer")
439  {
440  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
441  return false;
442  }
443  mmHeader.ctype=integer_type;
444  break;
445  case 'r' :
446  /* sanity check */
447  if(buffer != "real")
448  {
449  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
450  return false;
451  }
452  mmHeader.ctype=double_type;
453  break;
454  case 'c' :
455  /* sanity check */
456  if(buffer != "complex")
457  {
458  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
459  return false;
460  }
461  mmHeader.ctype=complex_type;
462  break;
463  case 'p' :
464  /* sanity check */
465  if(buffer != "pattern")
466  {
467  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
468  return false;
469  }
470  mmHeader.ctype=pattern;
471  break;
472  default :
473  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
474  return false;
475  }
476 
477  if(lineFeed(file))
478  return false;
479 
480  file >> buffer;
481 
482  std::transform(buffer.begin(), buffer.end(), buffer.begin(),
483  ::tolower);
484  switch(buffer[0])
485  {
486  case 'g' :
487  /* sanity check */
488  if(buffer != "general")
489  {
490  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
491  return false;
492  }
493  mmHeader.structure=general;
494  break;
495  case 'h' :
496  /* sanity check */
497  if(buffer != "hermitian")
498  {
499  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
500  return false;
501  }
502  mmHeader.structure=hermitian;
503  break;
504  case 's' :
505  if(buffer.size()==1) {
506  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
507  return false;
508  }
509 
510  switch(buffer[1])
511  {
512  case 'y' :
513  /* sanity check */
514  if(buffer != "symmetric")
515  {
516  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
517  return false;
518  }
519  mmHeader.structure=symmetric;
520  break;
521  case 'k' :
522  /* sanity check */
523  if(buffer != "skew-symmetric")
524  {
525  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
526  return false;
527  }
528  mmHeader.structure=skew_symmetric;
529  break;
530  default :
531  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
532  return false;
533  }
534  break;
535  default :
536  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
537  return false;
538  }
539  file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
540  c=file.peek();
541  return true;
542 
543  }
544 
545  template<std::size_t brows, std::size_t bcols>
546  std::tuple<std::size_t, std::size_t, std::size_t>
547  calculateNNZ(std::size_t rows, std::size_t cols, std::size_t entries, const MMHeader& header)
548  {
549  std::size_t blockrows=rows/brows;
550  std::size_t blockcols=cols/bcols;
551  std::size_t blocksize=brows*bcols;
552  std::size_t blockentries=0;
553 
554  switch(header.structure)
555  {
556  case general :
557  blockentries = entries/blocksize; break;
558  case skew_symmetric :
559  blockentries = 2*entries/blocksize; break;
560  case symmetric :
561  blockentries = (2*entries-rows)/blocksize; break;
562  case hermitian :
563  blockentries = (2*entries-rows)/blocksize; break;
564  default :
565  throw Dune::NotImplemented();
566  }
567  return std::make_tuple(blockrows, blockcols, blockentries);
568  }
569 
570  /*
571  * @brief Storage class for the column index and the numeric value.
572  *
573  * \tparam T Either a NumericWrapper of the numeric type or PatternDummy
574  * for MatrixMarket pattern case.
575  */
576  template<typename T>
577  struct IndexData : public T
578  {
579  std::size_t index = {};
580  };
581 
582 
593  template<typename T>
595  {
596  T number = {};
597  operator T&()
598  {
599  return number;
600  }
601  };
602 
607  {};
608 
609  template<>
611  {};
612 
613  template<typename T>
614  std::istream& operator>>(std::istream& is, NumericWrapper<T>& num)
615  {
616  return is>>num.number;
617  }
618 
619  inline std::istream& operator>>(std::istream& is, [[maybe_unused]] NumericWrapper<PatternDummy>& num)
620  {
621  return is;
622  }
623 
629  template<typename T>
630  bool operator<(const IndexData<T>& i1, const IndexData<T>& i2)
631  {
632  return i1.index<i2.index;
633  }
634 
640  template<typename T>
641  std::istream& operator>>(std::istream& is, IndexData<T>& data)
642  {
643  is>>data.index;
644  /* MatrixMarket indices are one based. Decrement for C++ */
645  --data.index;
646  return is>>data.number;
647  }
648 
654  template<typename T>
655  std::istream& operator>>(std::istream& is, IndexData<NumericWrapper<std::complex<T>>>& data)
656  {
657  is>>data.index;
658  /* MatrixMarket indices are one based. Decrement for C++ */
659  --data.index;
660  // real and imaginary part needs to be read separately as
661  // complex numbers are not provided in pair form. (x,y)
662  NumericWrapper<T> real, imag;
663  is>>real;
664  is>>imag;
665  data.number = {real.number, imag.number};
666  return is;
667  }
668 
675  template<typename D, int brows, int bcols>
677  {
683  template<typename T>
684  void operator()(const std::vector<std::set<IndexData<D> > >& rows,
685  BCRSMatrix<T>& matrix)
686  {
687  static_assert(IsNumber<T>::value && brows==1 && bcols==1, "Only scalar entries are expected here!");
688  for (auto iter=matrix.begin(); iter!= matrix.end(); ++iter)
689  {
690  auto brow=iter.index();
691  for (auto siter=rows[brow].begin(); siter != rows[brow].end(); ++siter)
692  (*iter)[siter->index] = siter->number;
693  }
694  }
695 
701  template<typename T>
702  void operator()(const std::vector<std::set<IndexData<D> > >& rows,
704  {
705  for (auto iter=matrix.begin(); iter!= matrix.end(); ++iter)
706  {
707  for (auto brow=iter.index()*brows,
708  browend=iter.index()*brows+brows;
709  brow<browend; ++brow)
710  {
711  for (auto siter=rows[brow].begin(), send=rows[brow].end();
712  siter != send; ++siter)
713  (*iter)[siter->index/bcols][brow%brows][siter->index%bcols]=siter->number;
714  }
715  }
716  }
717  };
718 
719  template<int brows, int bcols>
720  struct MatrixValuesSetter<PatternDummy,brows,bcols>
721  {
722  template<typename M>
723  void operator()(const std::vector<std::set<IndexData<PatternDummy> > >& rows,
724  M& matrix)
725  {}
726  };
727 
728  template<class T> struct is_complex : std::false_type {};
729  template<class T> struct is_complex<std::complex<T>> : std::true_type {};
730 
731  // wrapper for std::conj. Returns T if T is not complex.
732  template<class T>
733  std::enable_if_t<!is_complex<T>::value, T> conj(const T& r){
734  return r;
735  }
736 
737  template<class T>
738  std::enable_if_t<is_complex<T>::value, T> conj(const T& r){
739  return std::conj(r);
740  }
741 
742  template<typename M>
744  {};
745 
746  template<typename B, typename A>
748  {
749  enum {
750  rows = 1,
751  cols = 1
752  };
753  };
754 
755  template<typename B, int i, int j, typename A>
757  {
758  enum {
759  rows = i,
760  cols = j
761  };
762  };
763 
764  template<typename T, typename A, typename D>
766  std::istream& file, std::size_t entries,
767  const MMHeader& mmHeader, const D&)
768  {
770 
771  // Number of rows and columns of T, if it is a matrix (1x1 otherwise)
772  constexpr int brows = mm_multipliers<Matrix>::rows;
773  constexpr int bcols = mm_multipliers<Matrix>::cols;
774 
775  // First path
776  // store entries together with column index in a separate
777  // data structure
778  std::vector<std::set<IndexData<D> > > rows(matrix.N()*brows);
779 
780  auto readloop = [&] (auto symmetryFixup) {
781  for(std::size_t i = 0; i < entries; ++i) {
782  std::size_t row;
783  IndexData<D> data;
784  skipComments(file);
785  file>>row;
786  --row; // Index was 1 based.
787  assert(row/bcols<matrix.N());
788  file>>data;
789  assert(data.index/bcols<matrix.M());
790  rows[row].insert(data);
791  if(row!=data.index)
792  symmetryFixup(row, data);
793  }
794  };
795 
796  switch(mmHeader.structure)
797  {
798  case general:
799  readloop([](auto...){});
800  break;
801  case symmetric :
802  readloop([&](auto row, auto data) {
803  IndexData<D> data_sym(data);
804  data_sym.index = row;
805  rows[data.index].insert(data_sym);
806  });
807  break;
808  case skew_symmetric :
809  readloop([&](auto row, auto data) {
810  IndexData<D> data_sym;
811  data_sym.number = -data.number;
812  data_sym.index = row;
813  rows[data.index].insert(data_sym);
814  });
815  break;
816  case hermitian :
817  readloop([&](auto row, auto data) {
818  IndexData<D> data_sym;
819  data_sym.number = conj(data.number);
820  data_sym.index = row;
821  rows[data.index].insert(data_sym);
822  });
823  break;
824  default:
825  DUNE_THROW(Dune::NotImplemented,
826  "Only general, symmetric, skew-symmetric and hermitian is supported right now!");
827  }
828 
829  // Setup the matrix sparsity pattern
830  int nnz=0;
831  for(typename Matrix::CreateIterator iter=matrix.createbegin();
832  iter!= matrix.createend(); ++iter)
833  {
834  for(std::size_t brow=iter.index()*brows, browend=iter.index()*brows+brows;
835  brow<browend; ++brow)
836  {
837  typedef typename std::set<IndexData<D> >::const_iterator Siter;
838  for(Siter siter=rows[brow].begin(), send=rows[brow].end();
839  siter != send; ++siter, ++nnz)
840  iter.insert(siter->index/bcols);
841  }
842  }
843 
844  //Set the matrix values
845  matrix=0;
846 
848 
849  Setter(rows, matrix);
850  }
851 
852  inline std::tuple<std::string, std::string> splitFilename(const std::string& filename) {
853  std::size_t lastdot = filename.find_last_of(".");
854  if(lastdot == std::string::npos)
855  return std::make_tuple(filename, "");
856  else {
857  std::string potentialFileExtension = filename.substr(lastdot);
858  if (potentialFileExtension == ".mm" || potentialFileExtension == ".mtx")
859  return std::make_tuple(filename.substr(0, lastdot), potentialFileExtension);
860  else
861  return std::make_tuple(filename, "");
862  }
863  }
864 
865  } // end namespace MatrixMarketImpl
866 
867  class MatrixMarketFormatError : public Dune::Exception
868  {};
869 
870 
871  inline void mm_read_header(std::size_t& rows, std::size_t& cols,
872  MatrixMarketImpl::MMHeader& header, std::istream& istr,
873  bool isVector)
874  {
875  using namespace MatrixMarketImpl;
876 
877  if(!readMatrixMarketBanner(istr, header)) {
878  std::cerr << "First line was not a correct Matrix Market banner. Using default:\n"
879  << "%%MatrixMarket matrix coordinate real general"<<std::endl;
880  // Go to the beginning of the file
881  istr.clear() ;
882  istr.seekg(0, std::ios::beg);
883  if(isVector)
884  header.type=array_type;
885  }
886 
887  skipComments(istr);
888 
889  if(lineFeed(istr))
890  throw MatrixMarketFormatError();
891 
892  istr >> rows;
893 
894  if(lineFeed(istr))
895  throw MatrixMarketFormatError();
896  istr >> cols;
897  }
898 
899  template<typename T, typename A>
901  std::size_t size,
902  std::istream& istr,
903  size_t lane)
904  {
905  for (int i=0; size>0; ++i, --size)
906  istr>>Simd::lane(lane,vector[i]);
907  }
908 
909  template<typename T, typename A, int entries>
910  void mm_read_vector_entries(Dune::BlockVector<Dune::FieldVector<T,entries>,A>& vector,
911  std::size_t size,
912  std::istream& istr,
913  size_t lane)
914  {
915  for(int i=0; size>0; ++i, --size) {
916  Simd::Scalar<T> val;
917  istr>>val;
918  Simd::lane(lane, vector[i/entries][i%entries])=val;
919  }
920  }
921 
922 
929  template<typename T, typename A>
931  std::istream& istr)
932  {
933  typedef typename Dune::BlockVector<T,A>::field_type field_type;
934  using namespace MatrixMarketImpl;
935 
936  MMHeader header;
937  std::size_t rows, cols;
938  mm_read_header(rows,cols,header,istr, true);
939  if(cols!=Simd::lanes<field_type>()) {
940  if(Simd::lanes<field_type>() == 1)
941  DUNE_THROW(MatrixMarketFormatError, "cols!=1, therefore this is no vector!");
942  else
943  DUNE_THROW(MatrixMarketFormatError, "cols does not match the number of lanes in the field_type!");
944  }
945 
946  if(header.type!=array_type)
947  DUNE_THROW(MatrixMarketFormatError, "Vectors have to be stored in array format!");
948 
949 
950  if constexpr (Dune::IsNumber<T>())
951  vector.resize(rows);
952  else
953  {
954  T dummy;
955  auto blocksize = dummy.size();
956  std::size_t size=rows/blocksize;
957  if(size*blocksize!=rows)
958  DUNE_THROW(MatrixMarketFormatError, "Block size of vector is not correct!");
959 
960  vector.resize(size);
961  }
962 
963  istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
964  for(size_t l=0;l<Simd::lanes<field_type>();++l){
965  mm_read_vector_entries(vector, rows, istr, l);
966  }
967  }
968 
975  template<typename T, typename A>
977  std::istream& istr)
978  {
979  using namespace MatrixMarketImpl;
981 
982  MMHeader header;
983  if(!readMatrixMarketBanner(istr, header)) {
984  std::cerr << "First line was not a correct Matrix Market banner. Using default:\n"
985  << "%%MatrixMarket matrix coordinate real general"<<std::endl;
986  // Go to the beginning of the file
987  istr.clear() ;
988  istr.seekg(0, std::ios::beg);
989  }
990  skipComments(istr);
991 
992  std::size_t rows, cols, entries;
993 
994  if(lineFeed(istr))
995  throw MatrixMarketFormatError();
996 
997  istr >> rows;
998 
999  if(lineFeed(istr))
1000  throw MatrixMarketFormatError();
1001  istr >> cols;
1002 
1003  if(lineFeed(istr))
1004  throw MatrixMarketFormatError();
1005 
1006  istr >>entries;
1007 
1008  std::size_t nnz, blockrows, blockcols;
1009 
1010  // Number of rows and columns of T, if it is a matrix (1x1 otherwise)
1011  constexpr int brows = mm_multipliers<Matrix>::rows;
1012  constexpr int bcols = mm_multipliers<Matrix>::cols;
1013 
1014  std::tie(blockrows, blockcols, nnz) = calculateNNZ<brows, bcols>(rows, cols, entries, header);
1015 
1016  istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
1017 
1018 
1019  matrix.setSize(blockrows, blockcols, nnz);
1021 
1022  if(header.type==array_type)
1023  DUNE_THROW(Dune::NotImplemented, "Array format currently not supported for matrices!");
1024 
1025  readSparseEntries(matrix, istr, entries, header, NumericWrapper<typename Matrix::field_type>());
1026  }
1027 
1028  // Print a scalar entry
1029  template<typename B>
1030  void mm_print_entry(const B& entry,
1031  std::size_t rowidx,
1032  std::size_t colidx,
1033  std::ostream& ostr)
1034  {
1035  if constexpr (IsNumber<B>())
1036  ostr << rowidx << " " << colidx << " " << entry << std::endl;
1037  else
1038  {
1039  for (auto row=entry.begin(); row != entry.end(); ++row, ++rowidx) {
1040  int coli=colidx;
1041  for (auto col = row->begin(); col != row->end(); ++col, ++coli)
1042  ostr<< rowidx<<" "<<coli<<" "<<*col<<std::endl;
1043  }
1044  }
1045  }
1046 
1047  // Write a vector entry
1048  template<typename V>
1049  void mm_print_vector_entry(const V& entry, std::ostream& ostr,
1050  const std::integral_constant<int,1>&,
1051  size_t lane)
1052  {
1053  ostr<<Simd::lane(lane,entry)<<std::endl;
1054  }
1055 
1056  // Write a vector
1057  template<typename V>
1058  void mm_print_vector_entry(const V& vector, std::ostream& ostr,
1059  const std::integral_constant<int,0>&,
1060  size_t lane)
1061  {
1062  using namespace MatrixMarketImpl;
1063 
1064  // Is the entry a supported numeric type?
1065  const int isnumeric = mm_numeric_type<Simd::Scalar<typename V::block_type>>::is_numeric;
1066  typedef typename V::const_iterator VIter;
1067 
1068  for(VIter i=vector.begin(); i != vector.end(); ++i)
1069 
1070  mm_print_vector_entry(*i, ostr,
1071  std::integral_constant<int,isnumeric>(),
1072  lane);
1073  }
1074 
1075  template<typename T, typename A>
1076  std::size_t countEntries(const BlockVector<T,A>& vector)
1077  {
1078  return vector.size();
1079  }
1080 
1081  template<typename T, typename A, int i>
1082  std::size_t countEntries(const BlockVector<FieldVector<T,i>,A>& vector)
1083  {
1084  return vector.size()*i;
1085  }
1086 
1087  // Version for writing vectors.
1088  template<typename V>
1089  void writeMatrixMarket(const V& vector, std::ostream& ostr,
1090  const std::integral_constant<int,0>&)
1091  {
1092  using namespace MatrixMarketImpl;
1093  typedef typename V::field_type field_type;
1094 
1095  ostr<<countEntries(vector)<<" "<<Simd::lanes<field_type>()<<std::endl;
1096  const int isnumeric = mm_numeric_type<Simd::Scalar<V>>::is_numeric;
1097  for(size_t l=0;l<Simd::lanes<field_type>(); ++l){
1098  mm_print_vector_entry(vector,ostr, std::integral_constant<int,isnumeric>(), l);
1099  }
1100  }
1101 
1102  // Versions for writing matrices
1103  template<typename M>
1104  void writeMatrixMarket(const M& matrix,
1105  std::ostream& ostr,
1106  const std::integral_constant<int,1>&)
1107  {
1108  ostr<<matrix.N()*MatrixMarketImpl::mm_multipliers<M>::rows<<" "
1110  <<countNonZeros(matrix)<<std::endl;
1111 
1112  typedef typename M::const_iterator riterator;
1113  typedef typename M::ConstColIterator citerator;
1114  for(riterator row=matrix.begin(); row != matrix.end(); ++row)
1115  for(citerator col = row->begin(); col != row->end(); ++col)
1116  // Matrix Market indexing start with 1!
1119  }
1120 
1121 
1125  template<typename M>
1126  void writeMatrixMarket(const M& matrix,
1127  std::ostream& ostr)
1128  {
1129  using namespace MatrixMarketImpl;
1130 
1131  // Write header information
1132  mm_header_printer<M>::print(ostr);
1133  mm_block_structure_header<M>::print(ostr,matrix);
1134  // Choose the correct function for matrix and vector
1135  writeMatrixMarket(matrix,ostr,std::integral_constant<int,IsMatrix<M>::value>());
1136  }
1137 
1138  static const int default_precision = -1;
1150  template<typename M>
1151  void storeMatrixMarket(const M& matrix,
1152  std::string filename,
1153  int prec=default_precision)
1154  {
1155  auto [pureFilename, extension] = MatrixMarketImpl::splitFilename(filename);
1156  std::string rfilename;
1157  std::ofstream file;
1158  if (extension != "") {
1159  rfilename = pureFilename + extension;
1160  file.open(rfilename.c_str());
1161  if(!file)
1162  DUNE_THROW(IOError, "Could not open file for storage: " << rfilename.c_str());
1163  }
1164  else {
1165  // only try .mm so we do not ignore potential errors
1166  rfilename = pureFilename + ".mm";
1167  file.open(rfilename.c_str());
1168  if(!file)
1169  DUNE_THROW(IOError, "Could not open file for storage: " << rfilename.c_str());
1170  }
1171 
1172  file.setf(std::ios::scientific,std::ios::floatfield);
1173  if(prec>0)
1174  file.precision(prec);
1175  writeMatrixMarket(matrix, file);
1176  file.close();
1177  }
1178 
1179 #if HAVE_MPI
1194  template<typename M, typename G, typename L>
1195  void storeMatrixMarket(const M& matrix,
1196  std::string filename,
1198  bool storeIndices=true,
1199  int prec=default_precision)
1200  {
1201  // Get our rank
1202  int rank = comm.communicator().rank();
1203  // Write the local matrix
1204  auto [pureFilename, extension] = MatrixMarketImpl::splitFilename(filename);
1205  std::string rfilename;
1206  std::ofstream file;
1207  if (extension != "") {
1208  rfilename = pureFilename + "_" + std::to_string(rank) + extension;
1209  file.open(rfilename.c_str());
1210  dverb<< rfilename <<std::endl;
1211  if(!file)
1212  DUNE_THROW(IOError, "Could not open file for storage: " << rfilename.c_str());
1213  }
1214  else {
1215  // only try .mm so we do not ignore potential errors
1216  rfilename = pureFilename + "_" + std::to_string(rank) + ".mm";
1217  file.open(rfilename.c_str());
1218  dverb<< rfilename <<std::endl;
1219  if(!file)
1220  DUNE_THROW(IOError, "Could not open file for storage: " << rfilename.c_str());
1221  }
1222  file.setf(std::ios::scientific,std::ios::floatfield);
1223  if(prec>0)
1224  file.precision(prec);
1225  writeMatrixMarket(matrix, file);
1226  file.close();
1227 
1228  if(!storeIndices)
1229  return;
1230 
1231  // Write the global to local index mapping
1232  rfilename = pureFilename + "_" + std::to_string(rank) + ".idx";
1233  file.open(rfilename.c_str());
1234  if(!file)
1235  DUNE_THROW(IOError, "Could not open file for storage: " << rfilename.c_str());
1236  file.setf(std::ios::scientific,std::ios::floatfield);
1237  typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet IndexSet;
1238  typedef typename IndexSet::const_iterator Iterator;
1239  for(Iterator iter = comm.indexSet().begin();
1240  iter != comm.indexSet().end(); ++iter) {
1241  file << iter->global()<<" "<<(std::size_t)iter->local()<<" "
1242  <<(int)iter->local().attribute()<<" "<<(int)iter->local().isPublic()<<std::endl;
1243  }
1244  // Store neighbour information for efficient remote indices setup.
1245  file<<"neighbours:";
1246  const std::set<int>& neighbours=comm.remoteIndices().getNeighbours();
1247  typedef std::set<int>::const_iterator SIter;
1248  for(SIter neighbour=neighbours.begin(); neighbour != neighbours.end(); ++neighbour) {
1249  file<<" "<< *neighbour;
1250  }
1251  file.close();
1252  }
1253 
1268  template<typename M, typename G, typename L>
1269  void loadMatrixMarket(M& matrix,
1270  const std::string& filename,
1272  bool readIndices=true)
1273  {
1274  using namespace MatrixMarketImpl;
1275 
1277  typedef typename LocalIndexT::Attribute Attribute;
1278  // Get our rank
1279  int rank = comm.communicator().rank();
1280  // load local matrix
1281  auto [pureFilename, extension] = MatrixMarketImpl::splitFilename(filename);
1282  std::string rfilename;
1283  std::ifstream file;
1284  if (extension != "") {
1285  rfilename = pureFilename + "_" + std::to_string(rank) + extension;
1286  file.open(rfilename.c_str(), std::ios::in);
1287  dverb<< rfilename <<std::endl;
1288  if(!file)
1289  DUNE_THROW(IOError, "Could not open file: " << rfilename.c_str());
1290  }
1291  else {
1292  // try both .mm and .mtx
1293  rfilename = pureFilename + "_" + std::to_string(rank) + ".mm";
1294  file.open(rfilename.c_str(), std::ios::in);
1295  if(!file) {
1296  rfilename = pureFilename + "_" + std::to_string(rank) + ".mtx";
1297  file.open(rfilename.c_str(), std::ios::in);
1298  dverb<< rfilename <<std::endl;
1299  if(!file)
1300  DUNE_THROW(IOError, "Could not open file: " << rfilename.c_str());
1301  }
1302  }
1303  readMatrixMarket(matrix,file);
1304  file.close();
1305 
1306  if(!readIndices)
1307  return;
1308 
1309  // read indices
1310  typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet IndexSet;
1311  IndexSet& pis=comm.pis;
1312  rfilename = pureFilename + "_" + std::to_string(rank) + ".idx";
1313  file.open(rfilename.c_str());
1314  if(!file)
1315  DUNE_THROW(IOError, "Could not open file: " << rfilename.c_str());
1316  if(pis.size()!=0)
1317  DUNE_THROW(InvalidIndexSetState, "Index set is not empty!");
1318 
1319  pis.beginResize();
1320  while(!file.eof() && file.peek()!='n') {
1321  G g;
1322  file >>g;
1323  std::size_t l;
1324  file >>l;
1325  int c;
1326  file >>c;
1327  bool b;
1328  file >> b;
1329  pis.add(g,LocalIndexT(l,Attribute(c),b));
1330  lineFeed(file);
1331  }
1332  pis.endResize();
1333  if(!file.eof()) {
1334  // read neighbours
1335  std::string s;
1336  file>>s;
1337  if(s!="neighbours:")
1338  DUNE_THROW(MatrixMarketFormatError, "was expecting the string: \"neighbours:\"");
1339  std::set<int> nb;
1340  while(!file.eof()) {
1341  int i;
1342  file >> i;
1343  nb.insert(i);
1344  }
1345  file.close();
1346  comm.ri.setNeighbours(nb);
1347  }
1348  comm.ri.template rebuild<false>();
1349  }
1350 
1351  #endif
1352 
1363  template<typename M>
1364  void loadMatrixMarket(M& matrix,
1365  const std::string& filename)
1366  {
1367  auto [pureFilename, extension] = MatrixMarketImpl::splitFilename(filename);
1368  std::string rfilename;
1369  std::ifstream file;
1370  if (extension != "") {
1371  rfilename = pureFilename + extension;
1372  file.open(rfilename.c_str());
1373  if(!file)
1374  DUNE_THROW(IOError, "Could not open file: " << rfilename.c_str());
1375  }
1376  else {
1377  // try both .mm and .mtx
1378  rfilename = pureFilename + ".mm";
1379  file.open(rfilename.c_str(), std::ios::in);
1380  if(!file) {
1381  rfilename = pureFilename + ".mtx";
1382  file.open(rfilename.c_str(), std::ios::in);
1383  if(!file)
1384  DUNE_THROW(IOError, "Could not open file: " << rfilename.c_str());
1385  }
1386  }
1387  readMatrixMarket(matrix,file);
1388  file.close();
1389  }
1390 
1392 }
1393 #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:119
Col col
Definition: matrixmatrix.hh:351
void readMatrixMarket(Dune::BlockVector< T, A > &vector, std::istream &istr)
Reads a BlockVector from a matrix market file.
Definition: matrixmarket.hh:930
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:1151
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:1269
std::size_t countEntries(const BlockVector< T, A > &vector)
Definition: matrixmarket.hh:1076
void writeMatrixMarket(const V &vector, std::ostream &ostr, const std::integral_constant< int, 0 > &)
Definition: matrixmarket.hh:1089
void mm_print_vector_entry(const V &entry, std::ostream &ostr, const std::integral_constant< int, 1 > &, size_t lane)
Definition: matrixmarket.hh:1049
static const int default_precision
Definition: matrixmarket.hh:1138
void mm_read_vector_entries(Dune::BlockVector< T, A > &vector, std::size_t size, std::istream &istr, size_t lane)
Definition: matrixmarket.hh:900
void mm_read_header(std::size_t &rows, std::size_t &cols, MatrixMarketImpl::MMHeader &header, std::istream &istr, bool isVector)
Definition: matrixmarket.hh:871
void mm_print_entry(const B &entry, std::size_t rowidx, std::size_t colidx, std::ostream &ostr)
Definition: matrixmarket.hh:1030
Definition: allocator.hh:11
std::istream & operator>>(std::istream &is, NumericWrapper< T > &num)
Definition: matrixmarket.hh:614
bool operator<(const IndexData< T > &i1, const IndexData< T > &i2)
LessThan operator.
Definition: matrixmarket.hh:630
LineType
Definition: matrixmarket.hh:296
@ DATA
Definition: matrixmarket.hh:296
@ MM_HEADER
Definition: matrixmarket.hh:296
@ MM_ISTLSTRUCT
Definition: matrixmarket.hh:296
bool readMatrixMarketBanner(std::istream &file, MMHeader &mmHeader)
Definition: matrixmarket.hh:353
std::enable_if_t< is_complex< T >::value, T > conj(const T &r)
Definition: matrixmarket.hh:738
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:547
void readSparseEntries(Dune::BCRSMatrix< T, A > &matrix, std::istream &file, std::size_t entries, const MMHeader &mmHeader, const D &)
Definition: matrixmarket.hh:765
MM_TYPE
Definition: matrixmarket.hh:299
@ array_type
Definition: matrixmarket.hh:299
@ coordinate_type
Definition: matrixmarket.hh:299
@ unknown_type
Definition: matrixmarket.hh:299
std::tuple< std::string, std::string > splitFilename(const std::string &filename)
Definition: matrixmarket.hh:852
std::enable_if_t<!is_complex< T >::value, T > conj(const T &r)
Definition: matrixmarket.hh:733
void skipComments(std::istream &file)
Definition: matrixmarket.hh:339
bool lineFeed(std::istream &file)
Definition: matrixmarket.hh:315
@ MM_MAX_LINE_LENGTH
Definition: matrixmarket.hh:297
MM_STRUCTURE
Definition: matrixmarket.hh:303
@ skew_symmetric
Definition: matrixmarket.hh:303
@ general
Definition: matrixmarket.hh:303
@ hermitian
Definition: matrixmarket.hh:303
@ unknown_structure
Definition: matrixmarket.hh:303
@ symmetric
Definition: matrixmarket.hh:303
MM_CTYPE
Definition: matrixmarket.hh:301
@ unknown_ctype
Definition: matrixmarket.hh:301
@ pattern
Definition: matrixmarket.hh:301
@ complex_type
Definition: matrixmarket.hh:301
@ double_type
Definition: matrixmarket.hh:301
@ integer_type
Definition: matrixmarket.hh:301
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:466
Iterator begin()
Get iterator to first row.
Definition: bcrsmatrix.hh:675
Iterator end()
Get iterator to one beyond last row.
Definition: bcrsmatrix.hh:681
CreateIterator createend()
get create iterator pointing to one after the last block
Definition: bcrsmatrix.hh:1103
size_type M() const
number of columns (counted in blocks)
Definition: bcrsmatrix.hh:1978
CreateIterator createbegin()
get initial create iterator
Definition: bcrsmatrix.hh:1097
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1972
void setBuildMode(BuildMode bm)
Sets the build mode of the matrix.
Definition: bcrsmatrix.hh:833
void setSize(size_type rows, size_type columns, size_type nnz=0)
Set the size of the matrix.
Definition: bcrsmatrix.hh:861
A vector of blocks with memory management.
Definition: bvector.hh:395
void resize(size_type size)
Resize the vector.
Definition: bvector.hh:503
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: bvector.hh:401
A generic dynamic dense matrix.
Definition: matrix.hh:561
Helper metaprogram to get the matrix market string representation of the numeric type.
Definition: matrixmarket.hh:76
@ is_numeric
Whether T is a supported numeric type.
Definition: matrixmarket.hh:81
static std::string str()
Definition: matrixmarket.hh:95
static std::string str()
Definition: matrixmarket.hh:111
static std::string str()
Definition: matrixmarket.hh:127
static std::string str()
Definition: matrixmarket.hh:143
static std::string str()
Definition: matrixmarket.hh:159
Meta program to write the correct Matrix Market header.
Definition: matrixmarket.hh:174
static void print(std::ostream &os)
Definition: matrixmarket.hh:179
static void print(std::ostream &os)
Definition: matrixmarket.hh:189
static void print(std::ostream &os)
Definition: matrixmarket.hh:199
static void print(std::ostream &os)
Definition: matrixmarket.hh:209
Metaprogram for writing the ISTL block structure header.
Definition: matrixmarket.hh:225
static void print(std::ostream &os, const M &)
Definition: matrixmarket.hh:233
BlockVector< T, A > M
Definition: matrixmarket.hh:230
static void print(std::ostream &os, const M &)
Definition: matrixmarket.hh:245
BCRSMatrix< T, A > M
Definition: matrixmarket.hh:255
static void print(std::ostream &os, const M &)
Definition: matrixmarket.hh:258
static void print(std::ostream &os, const M &)
Definition: matrixmarket.hh:270
static void print(std::ostream &os, const M &m)
Definition: matrixmarket.hh:283
FieldMatrix< T, i, j > M
Definition: matrixmarket.hh:281
static void print(std::ostream &os, const M &m)
Definition: matrixmarket.hh:292
FieldVector< T, i > M
Definition: matrixmarket.hh:290
Definition: matrixmarket.hh:306
MM_STRUCTURE structure
Definition: matrixmarket.hh:312
MM_TYPE type
Definition: matrixmarket.hh:310
MMHeader()
Definition: matrixmarket.hh:307
MM_CTYPE ctype
Definition: matrixmarket.hh:311
Definition: matrixmarket.hh:578
std::size_t index
Definition: matrixmarket.hh:579
a wrapper class of numeric values.
Definition: matrixmarket.hh:595
T number
Definition: matrixmarket.hh:596
Utility class for marking the pattern type of the MatrixMarket matrices.
Definition: matrixmarket.hh:607
Functor to the data values of the matrix.
Definition: matrixmarket.hh:677
void operator()(const std::vector< std::set< IndexData< D > > > &rows, BCRSMatrix< T > &matrix)
Sets the matrix values.
Definition: matrixmarket.hh:684
void operator()(const std::vector< std::set< IndexData< D > > > &rows, BCRSMatrix< FieldMatrix< T, brows, bcols > > &matrix)
Sets the matrix values.
Definition: matrixmarket.hh:702
void operator()(const std::vector< std::set< IndexData< PatternDummy > > > &rows, M &matrix)
Definition: matrixmarket.hh:723
Definition: matrixmarket.hh:728
Definition: matrixmarket.hh:744
Definition: matrixmarket.hh:868
Definition: matrixutils.hh:27
Test whether a type is an ISTL Matrix.
Definition: matrixutils.hh:504
A class setting up standard communication for a two-valued attribute set with owner/overlap/copy sema...
Definition: owneroverlapcopy.hh:174
const ParallelIndexSet & indexSet() const
Get the underlying parallel index set.
Definition: owneroverlapcopy.hh:462
const Communication< MPI_Comm > & communicator() const
Definition: owneroverlapcopy.hh:299
const RemoteIndices & remoteIndices() const
Get the underlying remote indices.
Definition: owneroverlapcopy.hh:471
Dune::ParallelIndexSet< GlobalIdType, LI, 512 > ParallelIndexSet
The type of the parallel index set.
Definition: owneroverlapcopy.hh:449