dune-pdelab  2.0.0
petscmatrixbackend.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_PETSCMATRIXBACKEND_HH
4 #define DUNE_PETSCMATRIXBACKEND_HH
5 
6 #if HAVE_PETSC
7 
8 #include<vector>
9 #include <functional>
10 #include <algorithm>
11 #include <complex>
12 
16 
17 #include <petscmat.h>
18 
19 namespace Dune {
20  namespace PDELab {
21 
22  class PetscMatrixBackend;
23  class PetscMatrixContainer;
24 
25  class PetscNestedMatrixBackend;
26  class PetscNestedMatrixContainer;
27 
28  class PetscMatrixAccessorBase;
29 
30  template<typename LFSV, typename LFSU>
31  class PetscMatrixAccessor;
32 
33  template<typename LFSV, typename LFSU>
34  class PetscNestedMatrixAccessor;
35 
36  template<typename T>
37  class Pattern : public std::vector< std::set<T> >
38  {
39 
40  typedef std::vector< std::set<T> > BaseT;
41 
42  public:
43 
44  Pattern (T m_, T n_)
45  : BaseT(m_)
46  {}
47 
48  void add_link (T i, T j)
49  {
50  (*this)[i].insert(j);
51  }
52  };
53 
54  namespace {
55 
56  // Custom ordering functor for row-zeroing map.
57  // Required for complex support, as complex is not ordered by default.
58 
59 
60  // By default, just use std::less
61  template <class T>
62  struct MapOrdering
63  : public std::less<T>
64  {};
65 
66 
67  // For std::complex, implement a lexicographic ordering
68  template<typename T>
69  struct MapOrdering<std::complex<T> >
70  : public std::binary_function<
71  std::complex<T>,
72  std::complex<T>,
73  bool>
74  {
75 
76  bool operator() (const std::complex<T>& lhs, const std::complex<T>& rhs) {
77  return (lhs.real() < rhs.real()) ||
78  (!(rhs.real() < lhs.real()) && (lhs.imag() < rhs.imag()));
79  }
80 
81  };
82 
83  } // anonymous namespace
84 
85 
86  struct petsc_types
87  {
88  typedef std::size_t size_type;
89  typedef Dune::PDELab::Pattern<size_type> Pattern;
90  typedef shared_ptr<PetscMatrixContainer> MatrixPtr;
91  typedef std::vector<MatrixPtr> MatrixArray;
92  };
93 
94  template<typename GFSV, typename GFSVB, typename GFSU, typename GFSUB>
95  struct petsc_matrix_builder
96  {
97  dune_static_assert((AlwaysFalse<GFSV>::value),"Unsupported combination of trial and test space blocking");
98  };
99 
100  template<typename GFSV, typename GFSU>
101  petsc_types::MatrixPtr build_matrix(const GFSV& gfsv, const GFSU& gfsu, const petsc_types::Pattern& pattern)
102  {
103  return petsc_matrix_builder<GFSV,typename GFSV::Traits::BackendType,GFSU,typename GFSU::Traits::BackendType>::build(gfsv,gfsu,pattern);
104  }
105 
106  template<typename GFSV,typename GFSU>
107  struct petsc_matrix_builder<GFSV,PetscVectorBackend,GFSU,PetscVectorBackend>
108  : public petsc_types
109  {
110  static MatrixPtr build(const GFSV& gfsv, const GFSU& gfsu, const Pattern& pattern)
111  {
112  const size_type M = gfsv.size();
113  const size_type N = gfsu.size();
114 
115  std::vector<int> nnz(M);
116  std::vector<int>::iterator oit = nnz.begin();
117  for (Pattern::const_iterator it = pattern.begin(); it != pattern.end(); ++it, ++oit)
118  (*oit) = it->size();
119 
120  return make_shared<PetscMatrixContainer>(M,N,nnz);
121  }
122  };
123 
124  template<typename GFSV, typename GFSU, petsc_types::size_type row = 0, petsc_types::size_type col = 0>
125  struct recursive_matrix_builder;
126 
127  template<typename GFSV,typename GFSU>
128  struct petsc_matrix_builder<GFSV,PetscNestedVectorBackend,GFSU,PetscNestedVectorBackend>
129  : public petsc_types
130  {
131  static MatrixPtr build(const GFSV& gfsv, const GFSU& gfsu, const Pattern& pattern)
132  {
133  const size_type M = GFSV::CHILDREN;
134  const size_type N = GFSU::CHILDREN;
135 
136  MatrixArray matrices(N*M);
137  recursive_matrix_builder<GFSV,GFSU>::build_children(gfsv,gfsu,pattern,matrices);
138  return make_shared<PetscNestedMatrixContainer>(M,N,matrices);
139  }
140  };
141 
142  struct recursive_matrix_builder_base
143  : public petsc_types
144  {
145 
146  static Pattern extract_pattern(const Pattern& p, const size_type r_l, const size_type r_h, const size_type c_l, const size_type c_h)
147  {
148  Pattern out(r_h-r_l,c_h-c_l);
149  Pattern::const_iterator endit = p.begin() + r_h;
150  Pattern::iterator outit = out.begin();
151  typedef Pattern::value_type::const_iterator ColIterator;
152  for (Pattern::const_iterator rit = p.begin() + r_l; rit != endit; ++rit, ++outit)
153  {
154  for (ColIterator cit = rit->begin(); cit != rit->end(); ++cit)
155  {
156  if (*cit >= c_l && *cit < c_h)
157  outit->insert(*cit-c_l);
158  }
159  }
160  return std::move(out);
161  }
162 
163  };
164 
165  template<typename GFSV, typename GFSU, petsc_types::size_type row, petsc_types::size_type col>
166  struct recursive_matrix_builder
167  : public recursive_matrix_builder_base
168  {
169 
170  static const size_type M = GFSV::CHILDREN;
171  static const size_type N = GFSU::CHILDREN;
172 
173  static MatrixPtr build_child(const GFSV& gfsv, const GFSU& gfsu, const Pattern& pattern)
174  {
175  auto cgfsv = gfsv.template child<row>();
176  auto cgfsu = gfsu.template child<col>();
177  const size_type row_l = gfsv.ordering().subMap(row,0);
178  const size_type row_h = gfsv.ordering().subMap(row,cgfsv.ordering().size()-1) + 1;
179  const size_type col_l = gfsu.ordering().subMap(col,0);
180  const size_type col_h = gfsu.ordering().subMap(col,cgfsu.ordering().size()-1) + 1;
181  Pattern child_pattern = extract_pattern(pattern,row_l,row_h,col_l,col_h);
182  return build_matrix(cgfsv,cgfsu,child_pattern);
183  }
184 
185  template<std::size_t r = row>
186  static typename enable_if<(r < M) && (col < N-1)>::type
187  build_children(const GFSV& gfsv, const GFSU& gfsu, const Pattern& pattern, MatrixArray& matrices)
188  {
189  matrices[row*N + col] = build_child(gfsv,gfsu,pattern);
190  recursive_matrix_builder<GFSV,GFSU,row,col+1>::build_children(gfsv,gfsu,pattern,matrices);
191  }
192 
193  template<std::size_t r = row>
194  static typename enable_if<(r < M) && (col == N-1)>::type
195  build_children(const GFSV& gfsv, const GFSU& gfsu, const Pattern& pattern, MatrixArray& matrices)
196  {
197  matrices[row*N + col] = build_child(gfsv,gfsu,pattern);
198  recursive_matrix_builder<GFSV,GFSU,row+1,0>::build_children(gfsv,gfsu,pattern,matrices);
199  }
200 
201  template<std::size_t r = row>
202  static typename enable_if<(r == M)>::type
203  build_children(const GFSV& gfsv, const GFSU& gfsu, const Pattern& pattern, MatrixArray& matrices)
204  {
205  // end recursion
206  }
207 
208  };
209 
210  class PetscMatrixContainer
211  : public petsc_types
212  {
213 
214  friend class PetscMatrixBackend;
215 
216  friend class PetscMatrixAccessorBase;
217 
218  template<typename LFSV, typename LFSU>
219  friend class PetscMatrixAccessor;
220 
221  friend class PetscNestedMatrixContainer;
222 
223  enum AccessorState {
224  clean,
225  readValues,
226  setValues,
227  addValues,
228  };
229 
230  public:
231  typedef PetscScalar ElementType;
232  typedef ElementType E;
233  typedef Mat ContainerType;
234 
235  typedef E block_type;
236  typedef block_type field_type;
237 
238  typedef PetscMatrixBackend Backend;
239 
240 
241 
243  template<typename GO>
244  PetscMatrixContainer (const GO& go)
245  : _m(NULL)
246  , _accessorState(clean)
247  , _rowsToClear()
248  , _managed(true)
249  {
250  dune_static_assert((is_same<typename GO::Traits::MatrixBackend,PetscMatrixBackend>::value),"Wrong matrix backend type");
251 
252  Pattern pattern(go.globalSizeV(),go.globalSizeU());
253  go.fill_pattern(pattern);
254  MatrixPtr matrix(build_matrix(go.testGridFunctionSpace(),go.trialGridFunctionSpace(),pattern));
255  matrix->_managed = false;
256  _m = matrix->_m;
257  }
258 
259  PetscMatrixContainer (const size_type M, const size_type N, const std::vector<int>& nnz)
260  : _accessorState(clean)
261  , _rowsToClear()
262  , _managed(true)
263  {
264  PetscInt nz = 0;
265  PETSC_CALL(MatCreateSeqAIJ(PETSC_COMM_SELF,M,N,nz,&(nnz[0]),&_m));
266  PETSC_CALL(MatSetOption(_m,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE));
267  PETSC_CALL(MatSetOption(_m,MAT_NO_OFF_PROC_ZERO_ROWS,PETSC_TRUE)); // we only ever zero our own rows
268  }
269 
270  PetscMatrixContainer (const PetscMatrixContainer& rhs)
271  : _accessorState(clean)
272  , _rowsToClear()
273  , _managed(true)
274  {
275  PETSC_CALL(MatDuplicate(rhs._m,MAT_COPY_VALUES,&_m));
276  PETSC_CALL(MatSetOption(_m,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE));
277  PETSC_CALL(MatSetOption(_m,MAT_NO_OFF_PROC_ZERO_ROWS,PETSC_TRUE)); // we only ever zero our own rows
278  }
279 
280  PetscMatrixContainer (Mat m, bool managed = true)
281  : _m(m)
282  , _accessorState(clean)
283  , _rowsToClear()
284  , _managed(managed)
285  {}
286 
287  protected:
288 
289  PetscMatrixContainer (bool managed = true)
290  : _m(NULL)
291  , _accessorState(clean)
292  , _rowsToClear()
293  , _managed(managed)
294  {}
295 
296  public:
297 
298  virtual ~PetscMatrixContainer()
299  {
300  if (_managed)
301  PETSC_CALL(MatDestroy(&_m));
302  }
303 
304  PetscMatrixContainer& operator=(const PetscMatrixContainer& rhs)
305  {
306  assert(_accessorState == clean);
307  assert(_rowsToClear.empty());
308  PETSC_CALL(MatCopy(rhs._m,_m, DIFFERENT_NONZERO_PATTERN));
309  _managed = true;
310  }
311 
312  Mat& base ()
313  {
314  return _m;
315  }
316 
317  const Mat& base () const
318  {
319  return _m;
320  }
321 
322  operator Mat&()
323  {
324  return _m;
325  }
326 
327  operator const Mat&() const
328  {
329  return _m;
330  }
331 
332  size_type M() const
333  {
334  int size = 0;
335  PETSC_CALL(MatGetLocalSize(_m,&size,PETSC_NULL));
336  return size;
337  }
338 
339  size_type N() const
340  {
341  int size = 0;
342  PETSC_CALL(MatGetLocalSize(_m,PETSC_NULL,&size));
343  return size;
344  }
345 
346  PetscMatrixContainer& operator*= (const E& e)
347  {
348  PETSC_CALL(MatScale(_m,e));
349  return *this;
350  }
351 
352  protected:
353  Mat _m;
354  AccessorState _accessorState;
355 
356  std::map<PetscScalar,std::vector<int>,MapOrdering<PetscScalar> > _rowsToClear;
357 
358  bool _managed;
359 
360  private:
361 
362  void enqueue_row_clear(size_type i, PetscScalar diagonal_value)
363  {
364  _rowsToClear[diagonal_value].push_back(i);
365  }
366 
367  void access(AccessorState state)
368  {
369  if (_accessorState == clean)
370  _accessorState = state;
371  else if (_accessorState != state)
372  {
373  DUNE_THROW(InvalidStateException,"PETSc matrices do not allow mixing different access types between calls to flush()");
374  }
375  }
376 
377  virtual void flush(MatAssemblyType assemblyType)
378  {
379  if (_accessorState == addValues || _accessorState == setValues || assemblyType == MAT_FINAL_ASSEMBLY)
380  {
381  PETSC_CALL(MatAssemblyBegin(_m,assemblyType));
382  PETSC_CALL(MatAssemblyEnd(_m,assemblyType));
383  }
384  _accessorState = clean;
385  if (_rowsToClear.size() > 0) // TODO: communicate in MPI case
386  {
387  for (auto it = _rowsToClear.begin(); it != _rowsToClear.end(); ++it)
388  {
389  PETSC_CALL(MatZeroRows(_m,it->second.size(),&(it->second[0]),it->first,PETSC_NULL,PETSC_NULL));
390  }
391  _rowsToClear.clear();
392  }
393  }
394 
395 
396  void checkin() const
397  {
398  /*
399  if (_data)
400  {
401  std::cout << "restoring " << _data << std::endl;
402  PETSC_CALL(VecRestoreArray(_v,&_data));
403  }
404  */
405  }
406 
407  void checkout() const
408  {
409  /*
410  if (!_data)
411  {
412  PETSC_CALL(VecGetArray(_v,&_data));
413  std::cout << "got " << _data << std::endl;
414  }
415  */
416  }
417 
418 
419  };
420 
421 
422 
423  class PetscMatrixAccessorBase
424  : petsc_types
425  {
426 
427  template<typename,typename>
428  friend class PetscNestedMatrixAccessor;
429 
430  void setup(PetscMatrixContainer::AccessorState state)
431  {
432  if (_state == PetscMatrixContainer::clean)
433  {
434  _m.access(state);
435  if (state == PetscMatrixContainer::readValues)
436  {
437  PETSC_CALL(MatGetValues(_m.base(),_M,&(_rows[0]),_N,&(_cols[0]),&(_vals[0])));
438  }
439  _state = state;
440  } else if (_state != state)
441  {
442  DUNE_THROW(InvalidStateException,"Attempt to mix different access patterns within accessor");
443  }
444  }
445 
446  public:
447 
448  PetscMatrixAccessorBase(PetscMatrixContainer& m, size_type M, size_type N, size_type row_offset = 0, size_type col_offset = 0)
449  : _m(m)
450  , _M(M)
451  , _N(N)
452  , _rows(_M)
453  , _cols(_N)
454  , _vals(_M * _N)
455  , _state(PetscMatrixContainer::clean)
456  , _row_offset(row_offset)
457  , _col_offset(col_offset)
458  {}
459 
460  ~PetscMatrixAccessorBase()
461  {
462  switch (_state)
463  {
464  case PetscMatrixContainer::setValues:
465  PETSC_CALL(MatSetValues(_m.base(),_M,&(_rows[0]),_N,&(_cols[0]),&(_vals[0]),INSERT_VALUES));
466  case PetscMatrixContainer::addValues:
467  PETSC_CALL(MatSetValues(_m.base(),_M,&(_rows[0]),_N,&(_cols[0]),&(_vals[0]),ADD_VALUES));
468  default:
469  break;
470  }
471  }
472 
473  PetscScalar get(size_type i, size_type j)
474  {
475  setup(PetscMatrixContainer::readValues);
476  return _vals[i * _cols.size() + j];
477  }
478 
479  void set(size_type i, size_type j, PetscScalar v)
480  {
481  setup(PetscMatrixContainer::setValues);
482  _vals[i * _cols.size() + j] = v;
483  }
484 
485  void add(size_type i, size_type j, PetscScalar v)
486  {
487  setup(PetscMatrixContainer::addValues);
488  _vals[i * _cols.size() + j] += v;
489  }
490 
491  void setGlobal(size_type gi, size_type gj, const typename Matrix::field_type& v)
492  {
493  DUNE_THROW(NotImplemented,"Non-Dirichlet constraint is not yet implemented for PETSc backend.");
494  }
495 
496  void addGlobal(size_type gi, size_type gj, const typename Matrix::field_type& v)
497  {
498  DUNE_THROW(NotImplemented,"Non-Dirichlet constraint is not yet implemented for PETSc backend.");
499  }
500 
501  typename Matrix::field_type getGlobal(size_type gi, size_type gj) const
502  {
503  DUNE_THROW(NotImplemented,"Non-Dirichlet constraint is not yet implemented for PETSc backend.");
504  }
505 
506  private:
507 
508  PetscMatrixContainer& _m;
509 
510  protected:
511  const size_type _M;
512  const size_type _N;
513  std::vector<int> _rows;
514  std::vector<int> _cols;
515 
516  private:
517  std::vector<PetscScalar> _vals;
518  PetscMatrixContainer::AccessorState _state;
519  const size_type _row_offset;
520  const size_type _col_offset;
521 
522  };
523 
524 
525 
526  template<typename LFSV,typename LFSU>
527  class PetscMatrixAccessor
528  : public PetscMatrixAccessorBase
529  {
530 
531  friend class PetscMatrixBackend;
532 
533  protected:
534 
535  PetscMatrixAccessor(PetscMatrixContainer& m, const LFSV& lfsv, const LFSU& lfsu, size_type row_offset = 0, size_type col_offset = 0)
536  : PetscMatrixAccessorBase(m,lfsv.localVectorSize(),lfsu.localVectorSize(),row_offset,col_offset)
537  {
538  int m_offset = 0;
539  int n_offset = 0;
540  PETSC_CALL(MatGetOwnershipRange(m.base(),&m_offset,PETSC_NULL));
541  PETSC_CALL(MatGetOwnershipRangeColumn(m.base(),&n_offset,PETSC_NULL));
542  m_offset -= row_offset;
543  n_offset -= col_offset;
544  for (size_type i = 0; i < _M; ++i)
545  _rows[i] = m_offset + lfsv.globalIndex(i);
546  for (size_type i = 0; i < _N; ++i)
547  _cols[i] = n_offset + lfsu.globalIndex(i);
548  }
549 
550  };
551 
552 
553  class PetscNestedMatrixContainer
554  : public PetscMatrixContainer
555  {
556 
557  friend class PetscNestedMatrixBackend;
558 
559  template<typename LFSV, typename LFSU>
560  friend class PetscMatrixAccessor;
561 
562  public:
563  typedef PetscNestedMatrixBackend Backend;
564 
565 
566 
567  PetscNestedMatrixContainer(const size_type M, const size_type N, const MatrixArray& children)
568  {
569  Mat matrices[N*M];
570  for (size_type i = 0; i < N*M; ++i)
571  matrices[i] = children[i]->base();
572 
573  PETSC_CALL(MatCreateNest(PETSC_COMM_SELF,M,PETSC_NULL,N,PETSC_NULL,matrices,&_m));
574  PETSC_CALL(MatSetOption(_m,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE));
575  PETSC_CALL(MatSetOption(_m,MAT_NO_OFF_PROC_ZERO_ROWS,PETSC_TRUE)); // we only ever zero our own rows
576  }
577 
579  template<typename GO>
580  PetscNestedMatrixContainer (const GO& go)
581  {
582  dune_static_assert((is_same<typename GO::Traits::MatrixBackend,PetscNestedMatrixBackend>::value),"Wrong matrix backend type");
583 
584  Pattern pattern(go.globalSizeV(),go.globalSizeU());
585  go.fill_pattern(pattern);
586  MatrixPtr matrix(build_matrix(go.testGridFunctionSpace(),go.trialGridFunctionSpace(),pattern));
587  matrix->_managed = false;
588  _m = matrix->_m;
589  PETSC_CALL(MatSetOption(_m,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE));
590  PETSC_CALL(MatSetOption(_m,MAT_NO_OFF_PROC_ZERO_ROWS,PETSC_TRUE)); // we only ever zero our own rows
591  }
592 
593  virtual ~PetscNestedMatrixContainer()
594  {
595  }
596 
597  PetscMatrixContainer subMatrix(size_type i, size_type j)
598  {
599  Mat sub_mat;
600  PETSC_CALL(MatNestGetSubMat(_m,i,j,&sub_mat));
601  return PetscMatrixContainer(sub_mat,false);
602  }
603 
604  virtual void flush(MatAssemblyType assemblyType)
605  {
606  if (_accessorState == addValues || _accessorState == setValues || assemblyType == MAT_FINAL_ASSEMBLY)
607  {
608  PETSC_CALL(MatAssemblyBegin(_m,assemblyType));
609  PETSC_CALL(MatAssemblyEnd(_m,assemblyType));
610  }
611  _accessorState = clean;
612  if (_rowsToClear.size() > 0) // TODO: communicate in MPI case
613  {
614  int M;
615  int N;
616  Mat** sm;
617  PETSC_CALL(MatNestGetSubMats(_m,&M,&N,&sm));
618  int row_offsets[M+1];
619  int col_offsets[N+1];
620  row_offsets[0] = col_offsets[0] = 0;
621 
622  for (int i = 0; i < M; ++i)
623  {
624  int MM;
625  PETSC_CALL(MatGetLocalSize(sm[i][0],&MM,PETSC_NULL));
626  row_offsets[i+1] = row_offsets[i] + MM;
627  }
628 
629  for (int j = 0; j < N; ++j)
630  {
631  int NN;
632  PETSC_CALL(MatGetLocalSize(sm[0][j],PETSC_NULL,&NN));
633  col_offsets[j+1] = col_offsets[j] + NN;
634  }
635 
636  for (auto it = _rowsToClear.begin(); it != _rowsToClear.end(); ++it)
637  {
638  auto rows = it->second;
639  std::sort(rows.begin(),rows.end());
640  auto rowit = rows.begin();
641  auto rowend = rowit;
642  int row_block = 0;
643  do {
644  rowend = std::find_if(rowit,rows.end(),std::bind1st(std::less_equal<int>(),row_offsets[row_block+1]));
645  for (; rowit != rowend; ++rowit)
646  {
647  int local_row = (*rowit) - row_offsets[row_block];
648  for (int j = 0; j < N; ++j)
649  {
650  PETSC_CALL(MatZeroRows(sm[row_block][j],1,&local_row,0.0,PETSC_NULL,PETSC_NULL));
651  if (col_offsets[j] <= *rowit && *rowit < col_offsets[j+1])
652  {
653  PETSC_CALL(MatSetValue(sm[row_block][j],local_row,*rowit-col_offsets[j],it->first,INSERT_VALUES));
654  PETSC_CALL(MatAssemblyBegin(sm[row_block][j],MAT_FINAL_ASSEMBLY));
655  PETSC_CALL(MatAssemblyEnd(sm[row_block][j],MAT_FINAL_ASSEMBLY));
656  }
657  }
658  }
659  ++row_block;
660  } while (rowit != rows.end());
661  }
662  PETSC_CALL(MatAssemblyBegin(_m,MAT_FINAL_ASSEMBLY));
663  PETSC_CALL(MatAssemblyEnd(_m,MAT_FINAL_ASSEMBLY));
664  _rowsToClear.clear();
665  }
666  }
667 
668  };
669 
670 
671 
672  class PetscMatrixBackend
673  {
674  public:
675 
676  template<typename LFSV, typename LFSU, typename E>
677  class Accessor
678  : public PetscMatrixAccessor<LFSV,LFSU>
679  {
680 
681  public:
682 
683  Accessor(PetscMatrixContainer& m, const LFSV& lfsv, const LFSU& lfsu)
684  : PetscMatrixAccessor<LFSV,LFSU>(m,lfsv,lfsu)
685  {}
686 
687  };
688 
689  template<typename E>
690  class Matrix
691  : public PetscMatrixContainer
692  {
693 
694 
695  public:
696 
697  template<typename GridOperator>
698  explicit Matrix(const GridOperator& go)
699  : PetscMatrixContainer(go)
700  {}
701  };
702 
703  typedef petsc_types::Pattern Pattern;
704 
705  template<class C>
706  struct Value
707  {
708  typedef typename C::ElementType Type;
709  };
710 
712  typedef typename std::size_t size_type;
713 
714  static void clear_row (size_type i, PetscMatrixContainer& c, PetscScalar diagonal_value)
715  {
716  c.enqueue_row_clear(i,diagonal_value);
717  }
718 
719  static void flush(PetscMatrixContainer& c)
720  {
721  c.flush(MAT_FLUSH_ASSEMBLY);
722  }
723 
724  static void finalize(PetscMatrixContainer& c)
725  {
726  c.flush(MAT_FINAL_ASSEMBLY);
727  }
728 
729  };
730 
731 
732 
733  template<typename LFSV,typename LFSU>
734  class PetscNestedMatrixAccessor
735  : public petsc_types
736  {
737 
738  struct get_child_offsets
739  : public TypeTree::DirectChildrenVisitor
740  , public TypeTree::DynamicTraversal
741  {
742  get_child_offsets(std::vector<std::size_t>& local_offsets, std::vector<std::size_t>& global_offsets, std::vector<std::vector<int> >& indices)
743  : _local_offsets(local_offsets)
744  , _global_offsets(global_offsets)
745  , _indices(indices)
746  {}
747 
748  template<typename T, typename Child, typename TreePath, typename ChildIndex>
749  void beforeChild(const T& t, const Child& child, TreePath treePath, ChildIndex childIndex)
750  {
751  _local_offsets[childIndex+1] = _local_offsets[childIndex] + child.size();
752  _global_offsets[childIndex+1] = _global_offsets[childIndex] + child.gridFunctionSpace().size();
753  _indices[childIndex].resize(child.size());
754  for (std::size_t i = 0; i < child.size(); ++i)
755  _indices[childIndex][i] = child.globalIndex(i) - _global_offsets[childIndex];
756  }
757 
758  std::vector<std::size_t>& _local_offsets;
759  std::vector<std::size_t>& _global_offsets;
760  std::vector<std::vector<int> >& _indices;
761  };
762 
763  friend class PetscNestedMatrixBackend;
764 
765  static const std::size_t M = LFSV::CHILDREN;
766  static const std::size_t N = LFSU::CHILDREN;
767 
768  PetscNestedMatrixAccessor(PetscNestedMatrixContainer& m, const LFSV& lfsv, const LFSU& lfsu)
769  : _m(m)
770  , _lfsv(lfsv)
771  , _lfsu(lfsu)
772  , _accessors(N*M)
773  , _matrices(N*M)
774  , _local_row_offsets(M+1)
775  , _local_col_offsets(N+1)
776  , _global_row_offsets(M+1)
777  , _global_col_offsets(N+1)
778  , _row_indices(M)
779  , _col_indices(N)
780  {
781  get_child_offsets rows(_local_row_offsets,_global_row_offsets,_row_indices);
782  TypeTree::applyToTree(lfsv,rows);
783  get_child_offsets cols(_local_col_offsets,_global_col_offsets,_col_indices);
784  TypeTree::applyToTree(lfsu,cols);
785  }
786 
787  PetscMatrixAccessorBase& accessor(size_type& i, size_type& j)
788  {
789  auto row = std::find_if(_local_row_offsets.begin(),_local_row_offsets.end(),std::bind1st(std::less<size_type>(),i));
790  --row;
791  auto col = std::find_if(_local_col_offsets.begin(),_local_col_offsets.end(),std::bind1st(std::less<size_type>(),j));
792  --col;
793  const size_type mi = row - _local_row_offsets.begin();
794  const size_type mj = col - _local_col_offsets.begin();
795 
796  shared_ptr<PetscMatrixAccessorBase> a;
797  if (!(a = _accessors[mi*N + mj]))
798  {
799  Mat submat;
800  PETSC_CALL(MatNestGetSubMat(_m.base(),mi,mj,&submat));
801  shared_ptr<PetscMatrixContainer> c(make_shared<PetscMatrixContainer>(submat,false));
802  a = make_shared<PetscMatrixAccessorBase>(*c,(*(row+1))-(*row),(*(col+1))-(*col),_global_row_offsets[mi],_global_col_offsets[mj]);
803  a->_rows = _row_indices[mi];
804  a->_cols = _col_indices[mj];
805  _accessors[mi*N + mj] = a;
806  _matrices[mi*N + mj] = c;
807  }
808  i -= *row;
809  j -= *col;
810  return *a;
811  }
812 
813  public:
814 
815  typedef PetscNestedMatrixContainer::size_type size_type;
816 
817  PetscScalar get(size_type i, size_type j)
818  {
819  PetscMatrixAccessorBase& a = accessor(i,j);
820  return a.get(i,j);
821  }
822 
823  void set(size_type i, size_type j, PetscScalar v)
824  {
825  PetscMatrixAccessorBase& a = accessor(i,j);
826  a.set(i,j,v);
827  }
828 
829  void add(size_type i, size_type j, PetscScalar v)
830  {
831  PetscMatrixAccessorBase& a = accessor(i,j);
832  a.add(i,j,v);
833  }
834 
835  private:
836 
837  PetscNestedMatrixContainer& _m;
838  const LFSV& _lfsv;
839  const LFSU& _lfsu;
840  std::vector<shared_ptr<PetscMatrixAccessorBase> > _accessors;
841  std::vector<shared_ptr<PetscMatrixContainer> > _matrices;
842  std::vector<std::size_t> _local_row_offsets;
843  std::vector<std::size_t> _local_col_offsets;
844  std::vector<std::size_t> _global_row_offsets;
845  std::vector<std::size_t> _global_col_offsets;
846  std::vector<std::vector<int> > _row_indices;
847  std::vector<std::vector<int> > _col_indices;
848 
849  };
850 
851 
852 
853 
854  class PetscNestedMatrixBackend
855  : public PetscMatrixBackend
856  {
857  public:
858 
859  template<typename E>
860  class Matrix
861  : public PetscNestedMatrixContainer
862  {
863 
864 
865  public:
866 
867  template<typename GridOperator>
868  Matrix(const GridOperator& go)
869  : PetscNestedMatrixContainer(go)
870  {}
871  };
872 
873 
874  template<typename LFSV, typename LFSU, typename E>
875  class Accessor
876  : public PetscNestedMatrixAccessor<LFSV,LFSU>
877  {
878 
879 
880  public:
881 
882  Accessor(PetscNestedMatrixContainer& m, const LFSV& lfsv, const LFSU& lfsu)
883  : PetscNestedMatrixAccessor<LFSV,LFSU>(m,lfsv,lfsu)
884  {}
885 
886  };
887 
888 
889  };
890 
891  } // namespace PDELab
892 } // namespace Dune
893 
894 #endif // HAVE_PETSC
895 
896 #endif // DUNE_PETSCMATRIXBACKEND_HH
static const unsigned int value
Definition: gridfunctionspace/tags.hh:175
R p(int i, D x)
Definition: qkdg.hh:62
const E & e
Definition: interpolate.hh:172