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  *
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
166  *
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>
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
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>
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>
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
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
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
352  {
353  std::string buffer;
354  char c;
355  file >> buffer;
356  c=buffer[0];
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  }
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  }
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  }
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  }
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  }
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  }
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  }
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  }
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  }
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  }
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>
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
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,
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;
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
796  {
797  case general:
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
858  bool isVector)
859  {
860  using namespace MatrixMarketImpl;
861
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)
870  }
871
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>
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
919  std::size_t rows, cols;
921  if(cols!=1)
922  DUNE_THROW(MatrixMarketFormatError, "cols!=1, therefore this is no vector!");
923
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');
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
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  }
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
1000  DUNE_THROW(Dune::NotImplemented, "Array format currently not supported for matrices!");
1001
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)
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
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.
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>
1203  const std::string& filename,
1204  OwnerOverlapCopyCommunication<G,L>& comm,
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();
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  //
1223  file.close();
1224
1226  return;
1227
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;
1248  lineFeed(file);
1249  }
1250  pis.endResize();
1251  if(!file.eof()) {
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>
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);
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_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.
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)
std::enable_if_t<!is_complex< T >::value, T > conj(const T &r)
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