3 #ifndef DUNE_PETSCMATRIXBACKEND_HH
4 #define DUNE_PETSCMATRIXBACKEND_HH
22 class PetscMatrixBackend;
23 class PetscMatrixContainer;
25 class PetscNestedMatrixBackend;
26 class PetscNestedMatrixContainer;
28 class PetscMatrixAccessorBase;
30 template<
typename LFSV,
typename LFSU>
31 class PetscMatrixAccessor;
33 template<
typename LFSV,
typename LFSU>
34 class PetscNestedMatrixAccessor;
37 class Pattern :
public std::vector< std::set<T> >
40 typedef std::vector< std::set<T> > BaseT;
48 void add_link (T i, T j)
69 struct MapOrdering<std::complex<T> >
70 :
public std::binary_function<
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()));
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;
94 template<
typename GFSV,
typename GFSVB,
typename GFSU,
typename GFSUB>
95 struct petsc_matrix_builder
100 template<
typename GFSV,
typename GFSU>
101 petsc_types::MatrixPtr build_matrix(
const GFSV& gfsv,
const GFSU& gfsu,
const petsc_types::Pattern& pattern)
103 return petsc_matrix_builder<GFSV,typename GFSV::Traits::BackendType,GFSU,typename GFSU::Traits::BackendType>::build(gfsv,gfsu,pattern);
106 template<
typename GFSV,
typename GFSU>
107 struct petsc_matrix_builder<GFSV,PetscVectorBackend,GFSU,PetscVectorBackend>
110 static MatrixPtr build(
const GFSV& gfsv,
const GFSU& gfsu,
const Pattern& pattern)
112 const size_type M = gfsv.size();
113 const size_type N = gfsu.size();
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)
120 return make_shared<PetscMatrixContainer>(M,N,nnz);
124 template<
typename GFSV,
typename GFSU, petsc_types::
size_type row = 0, petsc_types::
size_type col = 0>
125 struct recursive_matrix_builder;
127 template<
typename GFSV,
typename GFSU>
128 struct petsc_matrix_builder<GFSV,PetscNestedVectorBackend,GFSU,PetscNestedVectorBackend>
131 static MatrixPtr build(
const GFSV& gfsv,
const GFSU& gfsu,
const Pattern& pattern)
133 const size_type M = GFSV::CHILDREN;
134 const size_type N = GFSU::CHILDREN;
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);
142 struct recursive_matrix_builder_base
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)
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)
154 for (ColIterator cit = rit->begin(); cit != rit->end(); ++cit)
156 if (*cit >= c_l && *cit < c_h)
157 outit->insert(*cit-c_l);
160 return std::move(out);
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
170 static const size_type M = GFSV::CHILDREN;
171 static const size_type N = GFSU::CHILDREN;
173 static MatrixPtr build_child(
const GFSV& gfsv,
const GFSU& gfsu,
const Pattern& pattern)
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);
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)
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);
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)
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);
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)
210 class PetscMatrixContainer
214 friend class PetscMatrixBackend;
216 friend class PetscMatrixAccessorBase;
218 template<
typename LFSV,
typename LFSU>
219 friend class PetscMatrixAccessor;
221 friend class PetscNestedMatrixContainer;
231 typedef PetscScalar ElementType;
232 typedef ElementType E;
233 typedef Mat ContainerType;
235 typedef E block_type;
236 typedef block_type field_type;
238 typedef PetscMatrixBackend Backend;
243 template<
typename GO>
244 PetscMatrixContainer (
const GO& go)
246 , _accessorState(clean)
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;
259 PetscMatrixContainer (
const size_type M,
const size_type N,
const std::vector<int>& nnz)
260 : _accessorState(clean)
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));
270 PetscMatrixContainer (
const PetscMatrixContainer& rhs)
271 : _accessorState(clean)
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));
280 PetscMatrixContainer (Mat m,
bool managed =
true)
282 , _accessorState(clean)
289 PetscMatrixContainer (
bool managed =
true)
291 , _accessorState(clean)
298 virtual ~PetscMatrixContainer()
301 PETSC_CALL(MatDestroy(&_m));
304 PetscMatrixContainer& operator=(
const PetscMatrixContainer& rhs)
306 assert(_accessorState == clean);
307 assert(_rowsToClear.empty());
308 PETSC_CALL(MatCopy(rhs._m,_m, DIFFERENT_NONZERO_PATTERN));
317 const Mat& base ()
const
327 operator const Mat&()
const
335 PETSC_CALL(MatGetLocalSize(_m,&size,PETSC_NULL));
342 PETSC_CALL(MatGetLocalSize(_m,PETSC_NULL,&size));
346 PetscMatrixContainer& operator*= (
const E&
e)
348 PETSC_CALL(MatScale(_m,e));
354 AccessorState _accessorState;
356 std::map<PetscScalar,std::vector<int>,MapOrdering<PetscScalar> > _rowsToClear;
362 void enqueue_row_clear(size_type i, PetscScalar diagonal_value)
364 _rowsToClear[diagonal_value].push_back(i);
367 void access(AccessorState state)
369 if (_accessorState == clean)
370 _accessorState = state;
371 else if (_accessorState != state)
373 DUNE_THROW(InvalidStateException,
"PETSc matrices do not allow mixing different access types between calls to flush()");
377 virtual void flush(MatAssemblyType assemblyType)
379 if (_accessorState == addValues || _accessorState == setValues || assemblyType == MAT_FINAL_ASSEMBLY)
381 PETSC_CALL(MatAssemblyBegin(_m,assemblyType));
382 PETSC_CALL(MatAssemblyEnd(_m,assemblyType));
384 _accessorState = clean;
385 if (_rowsToClear.size() > 0)
387 for (
auto it = _rowsToClear.begin(); it != _rowsToClear.end(); ++it)
389 PETSC_CALL(MatZeroRows(_m,it->second.size(),&(it->second[0]),it->first,PETSC_NULL,PETSC_NULL));
391 _rowsToClear.clear();
407 void checkout()
const
423 class PetscMatrixAccessorBase
427 template<
typename,
typename>
428 friend class PetscNestedMatrixAccessor;
430 void setup(PetscMatrixContainer::AccessorState state)
432 if (_state == PetscMatrixContainer::clean)
435 if (state == PetscMatrixContainer::readValues)
437 PETSC_CALL(MatGetValues(_m.base(),_M,&(_rows[0]),_N,&(_cols[0]),&(_vals[0])));
440 }
else if (_state != state)
442 DUNE_THROW(InvalidStateException,
"Attempt to mix different access patterns within accessor");
448 PetscMatrixAccessorBase(PetscMatrixContainer& m, size_type M, size_type N, size_type row_offset = 0, size_type col_offset = 0)
455 , _state(PetscMatrixContainer::clean)
456 , _row_offset(row_offset)
457 , _col_offset(col_offset)
460 ~PetscMatrixAccessorBase()
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));
473 PetscScalar
get(size_type i, size_type j)
475 setup(PetscMatrixContainer::readValues);
476 return _vals[i * _cols.size() + j];
479 void set(size_type i, size_type j, PetscScalar v)
481 setup(PetscMatrixContainer::setValues);
482 _vals[i * _cols.size() + j] = v;
485 void add(size_type i, size_type j, PetscScalar v)
487 setup(PetscMatrixContainer::addValues);
488 _vals[i * _cols.size() + j] += v;
491 void setGlobal(size_type gi, size_type gj,
const typename Matrix::field_type& v)
493 DUNE_THROW(NotImplemented,
"Non-Dirichlet constraint is not yet implemented for PETSc backend.");
496 void addGlobal(size_type gi, size_type gj,
const typename Matrix::field_type& v)
498 DUNE_THROW(NotImplemented,
"Non-Dirichlet constraint is not yet implemented for PETSc backend.");
501 typename Matrix::field_type getGlobal(size_type gi, size_type gj)
const
503 DUNE_THROW(NotImplemented,
"Non-Dirichlet constraint is not yet implemented for PETSc backend.");
508 PetscMatrixContainer& _m;
513 std::vector<int> _rows;
514 std::vector<int> _cols;
517 std::vector<PetscScalar> _vals;
518 PetscMatrixContainer::AccessorState _state;
519 const size_type _row_offset;
520 const size_type _col_offset;
526 template<
typename LFSV,
typename LFSU>
527 class PetscMatrixAccessor
528 :
public PetscMatrixAccessorBase
531 friend class PetscMatrixBackend;
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)
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);
553 class PetscNestedMatrixContainer
554 :
public PetscMatrixContainer
557 friend class PetscNestedMatrixBackend;
559 template<
typename LFSV,
typename LFSU>
560 friend class PetscMatrixAccessor;
563 typedef PetscNestedMatrixBackend Backend;
567 PetscNestedMatrixContainer(
const size_type M,
const size_type N,
const MatrixArray& children)
570 for (size_type i = 0; i < N*M; ++i)
571 matrices[i] = children[i]->base();
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));
579 template<
typename GO>
580 PetscNestedMatrixContainer (
const GO& go)
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;
589 PETSC_CALL(MatSetOption(_m,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE));
590 PETSC_CALL(MatSetOption(_m,MAT_NO_OFF_PROC_ZERO_ROWS,PETSC_TRUE));
593 virtual ~PetscNestedMatrixContainer()
597 PetscMatrixContainer subMatrix(size_type i, size_type j)
600 PETSC_CALL(MatNestGetSubMat(_m,i,j,&sub_mat));
601 return PetscMatrixContainer(sub_mat,
false);
604 virtual void flush(MatAssemblyType assemblyType)
606 if (_accessorState == addValues || _accessorState == setValues || assemblyType == MAT_FINAL_ASSEMBLY)
608 PETSC_CALL(MatAssemblyBegin(_m,assemblyType));
609 PETSC_CALL(MatAssemblyEnd(_m,assemblyType));
611 _accessorState = clean;
612 if (_rowsToClear.size() > 0)
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;
622 for (
int i = 0; i < M; ++i)
625 PETSC_CALL(MatGetLocalSize(sm[i][0],&MM,PETSC_NULL));
626 row_offsets[i+1] = row_offsets[i] + MM;
629 for (
int j = 0; j < N; ++j)
632 PETSC_CALL(MatGetLocalSize(sm[0][j],PETSC_NULL,&NN));
633 col_offsets[j+1] = col_offsets[j] + NN;
636 for (
auto it = _rowsToClear.begin(); it != _rowsToClear.end(); ++it)
638 auto rows = it->second;
639 std::sort(rows.begin(),rows.end());
640 auto rowit = rows.begin();
644 rowend = std::find_if(rowit,rows.end(),std::bind1st(std::less_equal<int>(),row_offsets[row_block+1]));
645 for (; rowit != rowend; ++rowit)
647 int local_row = (*rowit) - row_offsets[row_block];
648 for (
int j = 0; j < N; ++j)
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])
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));
660 }
while (rowit != rows.end());
662 PETSC_CALL(MatAssemblyBegin(_m,MAT_FINAL_ASSEMBLY));
663 PETSC_CALL(MatAssemblyEnd(_m,MAT_FINAL_ASSEMBLY));
664 _rowsToClear.clear();
672 class PetscMatrixBackend
676 template<
typename LFSV,
typename LFSU,
typename E>
678 :
public PetscMatrixAccessor<LFSV,LFSU>
683 Accessor(PetscMatrixContainer& m,
const LFSV& lfsv,
const LFSU& lfsu)
684 : PetscMatrixAccessor<LFSV,LFSU>(m,lfsv,lfsu)
691 :
public PetscMatrixContainer
697 template<
typename Gr
idOperator>
698 explicit Matrix(
const GridOperator& go)
699 : PetscMatrixContainer(go)
703 typedef petsc_types::Pattern Pattern;
708 typedef typename C::ElementType Type;
712 typedef typename std::size_t size_type;
714 static void clear_row (size_type i, PetscMatrixContainer& c, PetscScalar diagonal_value)
716 c.enqueue_row_clear(i,diagonal_value);
719 static void flush(PetscMatrixContainer& c)
721 c.flush(MAT_FLUSH_ASSEMBLY);
724 static void finalize(PetscMatrixContainer& c)
726 c.flush(MAT_FINAL_ASSEMBLY);
733 template<
typename LFSV,
typename LFSU>
734 class PetscNestedMatrixAccessor
738 struct get_child_offsets
739 :
public TypeTree::DirectChildrenVisitor
740 ,
public TypeTree::DynamicTraversal
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)
748 template<
typename T,
typename Child,
typename TreePath,
typename ChildIndex>
749 void beforeChild(
const T& t,
const Child& child, TreePath treePath, ChildIndex childIndex)
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];
758 std::vector<std::size_t>& _local_offsets;
759 std::vector<std::size_t>& _global_offsets;
760 std::vector<std::vector<int> >& _indices;
763 friend class PetscNestedMatrixBackend;
765 static const std::size_t M = LFSV::CHILDREN;
766 static const std::size_t N = LFSU::CHILDREN;
768 PetscNestedMatrixAccessor(PetscNestedMatrixContainer& m,
const LFSV& lfsv,
const LFSU& lfsu)
774 , _local_row_offsets(M+1)
775 , _local_col_offsets(N+1)
776 , _global_row_offsets(M+1)
777 , _global_col_offsets(N+1)
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);
787 PetscMatrixAccessorBase& accessor(size_type& i, size_type& j)
789 auto row = std::find_if(_local_row_offsets.begin(),_local_row_offsets.end(),std::bind1st(std::less<size_type>(),i));
791 auto col = std::find_if(_local_col_offsets.begin(),_local_col_offsets.end(),std::bind1st(std::less<size_type>(),j));
793 const size_type mi = row - _local_row_offsets.begin();
794 const size_type mj = col - _local_col_offsets.begin();
796 shared_ptr<PetscMatrixAccessorBase> a;
797 if (!(a = _accessors[mi*N + mj]))
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;
815 typedef PetscNestedMatrixContainer::size_type size_type;
817 PetscScalar
get(size_type i, size_type j)
819 PetscMatrixAccessorBase& a = accessor(i,j);
823 void set(size_type i, size_type j, PetscScalar v)
825 PetscMatrixAccessorBase& a = accessor(i,j);
829 void add(size_type i, size_type j, PetscScalar v)
831 PetscMatrixAccessorBase& a = accessor(i,j);
837 PetscNestedMatrixContainer& _m;
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;
854 class PetscNestedMatrixBackend
855 :
public PetscMatrixBackend
861 :
public PetscNestedMatrixContainer
867 template<
typename Gr
idOperator>
868 Matrix(
const GridOperator& go)
869 : PetscNestedMatrixContainer(go)
874 template<
typename LFSV,
typename LFSU,
typename E>
876 :
public PetscNestedMatrixAccessor<LFSV,LFSU>
882 Accessor(PetscNestedMatrixContainer& m,
const LFSV& lfsv,
const LFSU& lfsu)
883 : PetscNestedMatrixAccessor<LFSV,LFSU>(m,lfsv,lfsu)
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