dune-istl  2.7.1
About: dune-istl - DUNE (Distributed and Unified Numerics Environment) is a modular toolbox for solving partial differential equations (PDEs) with grid-based methods: DUNE Iterative Solver Template Library.
  Fossies Dox: dune-istl-2.7.1.tar.gz  ("unofficial" and yet experimental doxygen-generated source code documentation)  

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 
37  /**
38  * @defgroup ISTL_IO IO for matrices and vectors.
39  * @ingroup ISTL_SPMV
40  * @brief Provides methods for reading and writing matrices and vectors
41  * in various formats.
42  *
43  *
44  * Routine printmatrix prints a (sparse matrix with all entries (even zeroes).
45  * Function printvector prints a vector to a stream.
46  * PrintSparseMatrix prints a sparse matrix omitting all nonzeroes.
47  * With writeMatrixToMatlab one can write a matrix in a Matlab readable format.
48  * Using storeMartrixMarket and loadMatrixMarket one can store and load a parallel ISTL
49  * matrix in MatrixMarket format. The latter can even read a matrix written with
50  * writeMatrixToMatlab.
51  *
52  *
53  * @addtogroup ISTL_IO
54  * @{
55  */
56 
57  /** @file
58  * @author Markus Blatt
59  * @brief Provides classes for reading and writing MatrixMarket Files with
60  * an extension for parallel matrices.
61  */
62  namespace MatrixMarketImpl
63  {
64  /**
65  * @brief Helper metaprogram to get
66  * the matrix market string representation
67  * of the numeric type.
68  *
69  * Member function mm_numeric_type::c_str()
70  * returns the string representation of
71  * the type.
72  */
73  template<class T>
74  struct mm_numeric_type {
75  enum {
76  /**
77  * @brief Whether T is a supported numeric type.
78  */
79  is_numeric=false
80  };
81  };
82 
83  template<>
84  struct mm_numeric_type<int>
85  {
86  enum {
87  /**
88  * @brief Whether T is a supported numeric type.
89  */
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 {
103  /**
104  * @brief Whether T is a supported numeric type.
105  */
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 {
119  /**
120  * @brief Whether T is a supported numeric type.
121  */
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 {
135  /**
136  * @brief Whether T is a supported numeric type.
137  */
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 {
151  /**
152  * @brief Whether T is a supported numeric type.
153  */
154  is_numeric=true
155  };
156 
157  static std::string str()
158  {
159  return "complex";
160  }
161  };
162 
163  /**
164  * @brief Meta program to write the correct Matrix Market
165  * header.
166  *
167  * Member function mm_header::operator() writes the header.
168  *
169  * @tparam M The matrix type.
170  */
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 
214  /**
215  * @brief Metaprogram for writing the ISTL block
216  * structure header.
217  *
218  * Member function mm_block_structure_header::print(os, mat) writes
219  * the corresponding header to the specified ostream.
220  * @tparam The type of the matrix to generate the header for.
221  */
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 
581  /**
582  * @brief a wrapper class of numeric values.
583  *
584  * Use for template specialization to catch the pattern
585  * type MatrixMarket matrices. Uses Empty base class optimization
586  * in the pattern case.
587  *
588  * @tparam T Either a NumericWrapper of the numeric type or PatternDummy
589  * for MatrixMarket pattern case.
590  */
591  template<typename T>
593  {
595  operator T&()
596  {
597  return number;
598  }
599  };
600 
601  /**
602  * @brief Utility class for marking the pattern type of the MatrixMarket matrices.
603  */
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 
623  /**
624  * @brief LessThan operator.
625  *
626  * It simply compares the index.
627  */
628  template<typename T>
629  bool operator<(const IndexData<T>& i1, const IndexData<T>& i2)
630  {
631  return i1.index<i2.index;
632  }
633 
634  /**
635  * @brief Read IndexData from a stream.
636  * @param is The input stream we read.
637  * @param data Where to store the read data.
638  */
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 
648  /**
649  * @brief Read IndexData from a stream. Specialization for std::complex
650  * @param is The input stream we read.
651  * @param data Where to store the read data.
652  */
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 
668  /**
669  * @brief Functor to the data values of the matrix.
670  *
671  * This is specialized for PatternDummy. The specialization does not
672  * set anything.
673  */
674  template<typename D, int brows, int bcols>
676  {
677  /**
678  * @brief Sets the matrix values.
679  * @param row The row data as read from file.
680  * @param matrix The matrix whose data we set.
681  */
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 
695  /**
696  * @brief Sets the matrix values.
697  * @param row The row data as read from file.
698  * @param matrix The matrix whose data we set.
699  */
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 
906  /**
907  * @brief Reads a BlockVector from a matrix market file.
908  * @param vector The vector to store the data in.
909  * @param istr The input stream to read the data from.
910  * @warning Not all formats are supported!
911  */
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 
946  /**
947  * @brief Reads a sparse matrix from a matrix market file.
948  * @param matrix The matrix to store the data in.
949  * @param istr The input stream to read the data from.
950  * @warning Not all formats are supported!
951  */
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 
1094  /**
1095  * @brief writes a ISTL matrix or vector to a stream in matrix market format.
1096  */
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 
1111  /**
1112  * @brief Stores a parallel matrix/vector in matrix market format in a file.
1113  *
1114  * More about the matrix market exchange format can be found
1115  * <a href="http://math.nist.gov/MatrixMarket/formats.html">here</a>.
1116  *
1117  * @param matrix The matrix/vector to store.
1118  * @param filename the name of the filename (without suffix and rank!)
1119  * rank i will write the file filename_i.mm
1120  */
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
1132  /**
1133  * @brief Stores a parallel matrix/vector in matrix market format in a file.
1134  *
1135  * More about the matrix market exchange format can be found
1136  * <a href="http://math.nist.gov/MatrixMarket/formats.html">here</a>.
1137  *
1138  * @param matrix The matrix/vector to store.
1139  * @param filename the name of the filename (without suffix and rank!)
1140  * rank i will write the file filename_i.mm
1141  * @param comm The information about the data distribution.
1142  * @param storeIndices Whether to store the parallel index information.
1143  * If true rank i writes the index information to file filename_i.idx.
1144  */
1145  template<typename M, typename G, typename L>
1146  void storeMatrixMarket(const M& matrix,
1147  std::string filename,
1148  const OwnerOverlapCopyCommunication<G,L>& comm,
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 
1187  /**
1188  * @brief Load a parallel matrix/vector stored in matrix market format.
1189  *
1190  * More about the matrix market exchange format can be found
1191  * <a href="http://math.nist.gov/MatrixMarket/formats.html">here</a>.
1192  *
1193  * @param matrix Where to store the matrix/vector.
1194  * @param filename the name of the filename (without suffix and rank!)
1195  * rank i will read the file filename_i.mm
1196  * @param comm The information about the data distribution.
1197  * @param readIndices Whether to read the parallel index information.
1198  * If true rank i reads the index information form file filename_i.idx
1199  * And builds the remote indices information.
1200  */
1201  template<typename M, typename G, typename L>
1202  void loadMatrixMarket(M& matrix,
1203  const std::string& filename,
1204  OwnerOverlapCopyCommunication<G,L>& comm,
1205  bool readIndices=true)
1206  {
1207  using namespace MatrixMarketImpl;
1208 
1209  typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet::LocalIndex LocalIndex;
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 
1271  /**
1272  * @brief Load a matrix/vector stored in matrix market format.
1273  *
1274  * More about the matrix market exchange format can be found
1275  * <a href="http://math.nist.gov/MatrixMarket/formats.html">here</a>.
1276  *
1277  * @param matrix Where to store the matrix/vector.
1278  * @param filename the name of the filename (without suffix and rank!)
1279  * rank i will read the file filename_i.mm
1280  */
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 
1293  /** @} */
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...
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
size_type size() const
number of blocks in the array (are of size 1 here)
Definition: basearray.hh:250
A generic dynamic dense matrix.
Definition: matrix.hh:558
void readMatrixMarket(Dune::BlockVector< T, A > &vector, std::istream &istr)
Reads a BlockVector from a matrix market file.
std::size_t countEntries(const BlockVector< T, A > &vector)
void writeMatrixMarket(const V &vector, std::ostream &ostr, const std::integral_constant< int, 0 > &)
void mm_print_vector_entry(const V &entry, std::ostream &ostr, const std::integral_constant< int, 1 > &)
void storeMatrixMarket(const M &matrix, std::string filename)
Stores a parallel matrix/vector in matrix market format in a file.
void loadMatrixMarket(M &matrix, const std::string &filename)
Load a matrix/vector stored in matrix market format.
void mm_read_vector_entries(Dune::BlockVector< T, A > &vector, std::size_t size, std::istream &istr)
void mm_read_header(std::size_t &rows, std::size_t &cols, MatrixMarketImpl::MMHeader &header, std::istream &istr, bool isVector)
void mm_print_entry(const B &entry, std::size_t rowidx, std::size_t colidx, std::ostream &ostr)
Col col
Row row
Matrix & A
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
Some handy generic functions for ISTL matrices.
std::istream & operator>>(std::istream &is, NumericWrapper< T > &num)
bool operator<(const IndexData< T > &i1, const IndexData< T > &i2)
LessThan operator.
bool readMatrixMarketBanner(std::istream &file, MMHeader &mmHeader)
std::enable_if_t< is_complex< T >::value, T > conj(const T &r)
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)
void readSparseEntries(Dune::BCRSMatrix< T, A > &matrix, std::istream &file, std::size_t entries, const MMHeader &mmHeader, const D &)
std::enable_if_t<!is_complex< T >::value, T > conj(const T &r)
void skipComments(std::istream &file)
bool lineFeed(std::istream &file)
Definition: allocator.hh:7
Classes providing communication interfaces for overlapping Schwarz methods.
Test whether a type is an ISTL Matrix.
Definition: matrixutils.hh:502
void operator()(const std::vector< std::set< IndexData< PatternDummy > > > &rows, M &matrix)
Functor to the data values of the matrix.
void operator()(const std::vector< std::set< IndexData< D > > > &rows, BCRSMatrix< T > &matrix)
Sets the matrix values.
void operator()(const std::vector< std::set< IndexData< D > > > &rows, BCRSMatrix< FieldMatrix< T, brows, bcols > > &matrix)
Sets the matrix values.
a wrapper class of numeric values.
Utility class for marking the pattern type of the MatrixMarket matrices.
Metaprogram for writing the ISTL block structure header.
Meta program to write the correct Matrix Market header.
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