3#ifndef DUNE_ISTL_OVERLAPPINGSCHWARZ_HH
4#define DUNE_ISTL_OVERLAPPINGSCHWARZ_HH
14#include <dune/istl/bccsmatrixinitializer.hh>
37 template<
class M,
class X,
class TM,
class TD,
class TA>
38 class SeqOverlappingSchwarz;
43 template<
class I,
class S,
class D>
50 typedef I InitializerList;
51 typedef typename InitializerList::value_type AtomInitializer;
52 typedef typename AtomInitializer::Matrix Matrix;
53 typedef typename Matrix::const_iterator Iter;
54 typedef typename Matrix::row_type::const_iterator CIter;
57 typedef typename IndexSet::size_type size_type;
60 const IndexSet& indices,
64 void addRowNnz(
const Iter& row);
68 void countEntries(
const Iter& row,
const CIter& col)
const;
70 void calcColstart()
const;
72 void copyValue(
const Iter& row,
const CIter& col)
const;
74 void createMatrix()
const;
79 typedef typename S::size_type size_type;
80 typedef std::map<size_type,size_type> Map;
81 typedef typename Map::iterator iterator;
82 typedef typename Map::const_iterator const_iterator;
86 void insert(size_type grow);
88 const_iterator find(size_type grow)
const;
90 iterator find(size_type grow);
94 const_iterator begin()
const;
98 const_iterator end()
const;
101 std::map<size_type,size_type> map_;
106 typedef typename InitializerList::iterator InitIterator;
107 typedef typename IndexSet::const_iterator IndexIteratur;
108 InitializerList* initializers;
110 mutable std::vector<IndexMap> indexMaps;
137 template<
class M,
class X,
class Y>
141 template<
class K,
class Al,
class X,
class Y>
147 typedef typename std::remove_const<M>::type matrix_type;
148 typedef typename X::field_type field_type;
149 typedef typename std::remove_const<M>::type rilu_type;
151 typedef X domain_type;
153 typedef Y range_type;
154 static constexpr size_t n = std::decay_t<decltype(Impl::asMatrix(std::declval<K>()))>::rows;
162 assert(v.size() > 0);
163 assert(v.size() == d.size());
164 assert(A.rows() <= v.size());
165 assert(A.cols() <= v.size());
166 size_t sz = A.rows();
182 void setSubMatrix(
const M& BCRS, S& rowset)
184 size_t sz = rowset.size();
186 typedef typename S::const_iterator SIter;
188 for(SIter rowIdx = rowset.begin(), rowEnd=rowset.end();
189 rowIdx!= rowEnd; ++rowIdx, r++)
192 for(SIter colIdx = rowset.begin(), colEnd=rowset.end();
193 colIdx!= colEnd; ++colIdx, c++)
195 if (BCRS[*rowIdx].find(*colIdx) == BCRS[*rowIdx].end())
197 for (
size_t i=0; i<n; i++)
199 for (
size_t j=0; j<n; j++)
201 A[r*n+i][c*n+j] = Impl::asMatrix(BCRS[*rowIdx][*colIdx])[i][j];
211 template<
typename T,
bool tag>
212 class OverlappingAssignerHelper
216 using OverlappingAssigner = OverlappingAssignerHelper<T, Dune::StoresColumnCompressed<T>::value>;
219 template<
class K,
class Al,
class X,
class Y>
220 class OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix<K, Al>, X, Y >,false>
223 typedef BCRSMatrix< K, Al> matrix_type;
224 typedef typename X::field_type field_type;
225 typedef Y range_type;
226 typedef typename range_type::block_type block_type;
227 typedef typename matrix_type::size_type size_type;
228 static constexpr size_t n = std::decay_t<decltype(Impl::asMatrix(std::declval<K>()))>::rows;
236 OverlappingAssignerHelper(std::size_t maxlength,
const BCRSMatrix<K, Al>& mat_,
const X& b_, Y& x_);
248 void resetIndexForNextDomain();
255 DynamicVector<field_type> & lhs();
262 DynamicVector<field_type> & rhs();
269 void relaxResult(field_type relax);
275 void operator()(
const size_type& domainIndex);
285 void assignResult(block_type& res);
291 const matrix_type* mat;
294 DynamicVector<field_type> * rhs_;
297 DynamicVector<field_type> * lhs_;
305 std::size_t maxlength_;
308#if HAVE_SUPERLU || HAVE_SUITESPARSE_UMFPACK
309 template<
template<
class>
class S,
typename T,
typename A>
310 struct OverlappingAssignerHelper<S<BCRSMatrix<T, A>>, true>
312 typedef BCRSMatrix<T, A> matrix_type;
313 typedef typename S<BCRSMatrix<T, A>>::range_type range_type;
314 typedef typename range_type::field_type field_type;
315 typedef typename range_type::block_type block_type;
317 typedef typename matrix_type::size_type size_type;
319 static constexpr size_t n = std::decay_t<decltype(Impl::asMatrix(std::declval<T>()))>::rows;
320 static constexpr size_t m = std::decay_t<decltype(Impl::asMatrix(std::declval<T>()))>::cols;
328 OverlappingAssignerHelper(std::size_t maxlength,
const matrix_type& mat,
329 const range_type& b, range_type& x);
340 void resetIndexForNextDomain();
358 void relaxResult(field_type relax);
364 void operator()(
const size_type& domain);
373 void assignResult(block_type& res);
379 const matrix_type* mat;
391 std::size_t maxlength_;
396 template<
class M,
class X,
class Y>
397 class OverlappingAssignerILUBase
400 typedef M matrix_type;
402 typedef typename Y::field_type field_type;
404 typedef typename Y::block_type block_type;
406 typedef typename matrix_type::size_type size_type;
414 OverlappingAssignerILUBase(std::size_t maxlength,
const M& mat,
426 void resetIndexForNextDomain();
444 void relaxResult(field_type relax);
450 void operator()(
const size_type& domain);
459 void assignResult(block_type& res);
479 template<
class M,
class X,
class Y>
480 class OverlappingAssignerHelper<ILU0SubdomainSolver<M,X,Y>, false>
481 :
public OverlappingAssignerILUBase<M,X,Y>
491 OverlappingAssignerHelper(std::size_t maxlength,
const M& mat,
493 : OverlappingAssignerILUBase<M,X,Y>(maxlength, mat,b,x)
498 template<
class M,
class X,
class Y>
499 class OverlappingAssignerHelper<ILUNSubdomainSolver<M,X,Y>,false>
500 :
public OverlappingAssignerILUBase<M,X,Y>
510 OverlappingAssignerHelper(std::size_t maxlength,
const M& mat,
512 : OverlappingAssignerILUBase<M,X,Y>(maxlength, mat,b,x)
516 template<
typename S,
typename T>
520 template<
typename S,
typename T,
typename A>
521 struct AdditiveAdder<S, BlockVector<T,A> >
523 typedef typename A::size_type size_type;
524 typedef typename std::decay_t<decltype(Impl::asVector(std::declval<T>()))>::field_type field_type;
525 AdditiveAdder(BlockVector<T,A>& v, BlockVector<T,A>& x,
526 OverlappingAssigner<S>& assigner,
const field_type& relax_);
527 void operator()(
const size_type& domain);
529 static constexpr size_t n = std::decay_t<decltype(Impl::asVector(std::declval<T>()))>::dimension;
534 OverlappingAssigner<S>* assigner;
538 template<
typename S,
typename T>
539 struct MultiplicativeAdder
542 template<
typename S,
typename T,
typename A>
543 struct MultiplicativeAdder<S, BlockVector<T,A> >
545 typedef typename A::size_type size_type;
546 typedef typename std::decay_t<decltype(Impl::asVector(std::declval<T>()))>::field_type field_type;
547 MultiplicativeAdder(BlockVector<T,A>& v, BlockVector<T,A>& x,
548 OverlappingAssigner<S>& assigner_,
const field_type& relax_);
549 void operator()(
const size_type& domain);
551 static constexpr size_t n = std::decay_t<decltype(Impl::asVector(std::declval<T>()))>::dimension;
555 OverlappingAssigner<S>* assigner;
568 template<
typename T,
class X,
class S>
572 template<
class X,
class S>
575 typedef AdditiveAdder<S,X> Adder;
578 template<
class X,
class S>
579 struct AdderSelector<MultiplicativeSchwarzMode,X,S>
581 typedef MultiplicativeAdder<S,X> Adder;
584 template<
class X,
class S>
585 struct AdderSelector<SymmetricMultiplicativeSchwarzMode,X,S>
587 typedef MultiplicativeAdder<S,X> Adder;
601 template<
typename T1,
typename T2,
bool forward>
604 typedef T1 solver_vector;
605 typedef typename solver_vector::iterator solver_iterator;
606 typedef T2 subdomain_vector;
607 typedef typename subdomain_vector::const_iterator domain_iterator;
609 static solver_iterator begin(solver_vector& sv)
614 static solver_iterator end(solver_vector& sv)
618 static domain_iterator begin(
const subdomain_vector& sv)
623 static domain_iterator end(
const subdomain_vector& sv)
629 template<
typename T1,
typename T2>
632 typedef T1 solver_vector;
633 typedef typename solver_vector::reverse_iterator solver_iterator;
634 typedef T2 subdomain_vector;
635 typedef typename subdomain_vector::const_reverse_iterator domain_iterator;
637 static solver_iterator begin(solver_vector& sv)
642 static solver_iterator end(solver_vector& sv)
646 static domain_iterator begin(
const subdomain_vector& sv)
651 static domain_iterator end(
const subdomain_vector& sv)
669 typedef typename smoother::range_type range_type;
671 static void apply(smoother& sm, range_type& v,
const range_type& b)
673 sm.template apply<true>(v, b);
677 template<
class M,
class X,
class TD,
class TA>
681 typedef typename smoother::range_type range_type;
683 static void apply(smoother& sm, range_type& v,
const range_type& b)
685 sm.template apply<true>(v, b);
686 sm.template apply<false>(v, b);
690 template<
class T,
bool tag>
691 struct SeqOverlappingSchwarzAssemblerHelper
695 using SeqOverlappingSchwarzAssembler = SeqOverlappingSchwarzAssemblerHelper<T,Dune::StoresColumnCompressed<T>::value>;
697 template<
class K,
class Al,
class X,
class Y>
698 struct SeqOverlappingSchwarzAssemblerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,false>
700 typedef BCRSMatrix< K, Al> matrix_type;
701 static constexpr size_t n = std::decay_t<decltype(Impl::asMatrix(std::declval<K>()))>::rows;
702 template<
class RowToDomain,
class Solvers,
class SubDomains>
703 static std::size_t assembleLocalProblems(
const RowToDomain& rowToDomain,
const matrix_type& mat,
704 Solvers& solvers,
const SubDomains& domains,
708 template<
template<
class>
class S,
typename T,
typename A>
709 struct SeqOverlappingSchwarzAssemblerHelper<S<BCRSMatrix<T,A>>,true>
711 typedef BCRSMatrix<T,A> matrix_type;
712 static constexpr size_t n = std::decay_t<decltype(Impl::asMatrix(std::declval<T>()))>::rows;
713 template<
class RowToDomain,
class Solvers,
class SubDomains>
714 static std::size_t assembleLocalProblems(
const RowToDomain& rowToDomain,
const matrix_type& mat,
715 Solvers& solvers,
const SubDomains& domains,
719 template<
class M,
class X,
class Y>
720 struct SeqOverlappingSchwarzAssemblerILUBase
722 typedef M matrix_type;
723 template<
class RowToDomain,
class Solvers,
class SubDomains>
724 static std::size_t assembleLocalProblems(
const RowToDomain& rowToDomain,
const matrix_type& mat,
725 Solvers& solvers,
const SubDomains& domains,
729 template<
class M,
class X,
class Y>
730 struct SeqOverlappingSchwarzAssemblerHelper<ILU0SubdomainSolver<M,X,Y>,false>
731 :
public SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>
734 template<
class M,
class X,
class Y>
735 struct SeqOverlappingSchwarzAssemblerHelper<ILUNSubdomainSolver<M,X,Y>,false>
736 :
public SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>
749 template<
class M,
class X,
class TM=AdditiveSchwarzMode,
750 class TD=ILU0SubdomainSolver<M,X,X>,
class TA=std::allocator<X> >
790 typedef std::set<size_type, std::less<size_type>,
791 typename std::allocator_traits<TA>::template rebind_alloc<size_type> >
795 typedef std::vector<subdomain_type, typename std::allocator_traits<TA>::template rebind_alloc<subdomain_type> >
subdomain_vector;
801 typedef std::vector<subdomain_list, typename std::allocator_traits<TA>::template rebind_alloc<subdomain_list> >
rowtodomain_vector;
807 typedef std::vector<slu, typename std::allocator_traits<TA>::template rebind_alloc<slu> >
slu_vector;
823 field_type relaxationFactor=1,
bool onTheFly_=
true);
837 field_type relaxationFactor=1,
bool onTheFly_=
true);
844 virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] X& b)
852 virtual void apply (X& v,
const X& d);
859 virtual void post ([[maybe_unused]] X& x)
862 template<
bool forward>
877 typename M::size_type maxlength;
884 template<
class I,
class S,
class D>
885 OverlappingSchwarzInitializer<I,S,D>::OverlappingSchwarzInitializer(InitializerList& il,
887 const subdomain_vector& domains_)
888 : initializers(&il), indices(&idx), indexMaps(il.size()), domains(domains_)
892 template<
class I,
class S,
class D>
893 void OverlappingSchwarzInitializer<I,S,D>::addRowNnz(
const Iter& row)
895 typedef typename IndexSet::value_type::const_iterator iterator;
896 for(iterator domain=(*indices)[row.index()].begin(); domain != (*indices)[row.index()].end(); ++domain) {
897 (*initializers)[*domain].addRowNnz(row, domains[*domain]);
898 indexMaps[*domain].insert(row.index());
902 template<
class I,
class S,
class D>
903 void OverlappingSchwarzInitializer<I,S,D>::allocate()
905 for(
auto&& i: *initializers)
906 i.allocateMatrixStorage();
907 for(
auto&& i: *initializers)
911 template<
class I,
class S,
class D>
912 void OverlappingSchwarzInitializer<I,S,D>::countEntries(
const Iter& row,
const CIter& col)
const
914 typedef typename IndexSet::value_type::const_iterator iterator;
915 for(iterator domain=(*indices)[row.index()].begin(); domain != (*indices)[row.index()].end(); ++domain) {
916 typename std::map<size_type,size_type>::const_iterator v = indexMaps[*domain].
find(col.index());
917 if(v!= indexMaps[*domain].end()) {
918 (*initializers)[*domain].countEntries(indexMaps[*domain].find(col.index())->second);
923 template<
class I,
class S,
class D>
924 void OverlappingSchwarzInitializer<I,S,D>::calcColstart()
const
926 for(
auto&& i : *initializers)
930 template<
class I,
class S,
class D>
931 void OverlappingSchwarzInitializer<I,S,D>::copyValue(
const Iter& row,
const CIter& col)
const
933 typedef typename IndexSet::value_type::const_iterator iterator;
934 for(iterator domain=(*indices)[row.index()].begin(); domain!= (*indices)[row.index()].end(); ++domain) {
935 typename std::map<size_type,size_type>::const_iterator v = indexMaps[*domain].find(col.index());
936 if(v!= indexMaps[*domain].end()) {
937 assert(indexMaps[*domain].end()!=indexMaps[*domain].find(row.index()));
938 (*initializers)[*domain].copyValue(col, indexMaps[*domain].find(row.index())->second,
944 template<
class I,
class S,
class D>
945 void OverlappingSchwarzInitializer<I,S,D>::createMatrix()
const
947 std::vector<IndexMap>().swap(indexMaps);
948 for(
auto&& i: *initializers)
952 template<
class I,
class S,
class D>
953 OverlappingSchwarzInitializer<I,S,D>::IndexMap::IndexMap()
957 template<
class I,
class S,
class D>
958 void OverlappingSchwarzInitializer<I,S,D>::IndexMap::insert(size_type grow)
960 assert(map_.find(grow)==map_.end());
961 map_.insert(std::make_pair(grow, row++));
964 template<
class I,
class S,
class D>
965 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator
966 OverlappingSchwarzInitializer<I,S,D>::IndexMap::find(size_type grow)
const
968 return map_.find(grow);
971 template<
class I,
class S,
class D>
972 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator
973 OverlappingSchwarzInitializer<I,S,D>::IndexMap::find(size_type grow)
975 return map_.find(grow);
978 template<
class I,
class S,
class D>
979 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator
980 OverlappingSchwarzInitializer<I,S,D>::IndexMap::end()
const
985 template<
class I,
class S,
class D>
986 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator
987 OverlappingSchwarzInitializer<I,S,D>::IndexMap::end()
992 template<
class I,
class S,
class D>
993 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator
994 OverlappingSchwarzInitializer<I,S,D>::IndexMap::begin()
const
999 template<
class I,
class S,
class D>
1000 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator
1001 OverlappingSchwarzInitializer<I,S,D>::IndexMap::begin()
1003 return map_.begin();
1006 template<
class M,
class X,
class TM,
class TD,
class TA>
1009 : mat(mat_), relax(relaxationFactor), onTheFly(fly)
1011 typedef typename rowtodomain_vector::const_iterator RowDomainIterator;
1013#ifdef DUNE_ISTL_WITH_CHECKING
1014 assert(rowToDomain.size()==mat.N());
1015 assert(rowToDomain.size()==mat.M());
1017 for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter)
1018 assert(iter->size()>0);
1023 for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter)
1024 for(DomainIterator d=iter->
begin(); d != iter->
end(); ++d)
1028 solvers.resize(domains);
1029 subDomains.resize(domains);
1033 for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter, ++row)
1034 for(DomainIterator d=iter->
begin(); d != iter->
end(); ++d)
1035 subDomains[*d].insert(row);
1037#ifdef DUNE_ISTL_WITH_CHECKING
1039 typedef typename subdomain_vector::const_iterator iterator;
1040 for(iterator iter=subDomains.begin(); iter != subDomains.end(); ++iter) {
1041 typedef typename subdomain_type::const_iterator entry_iterator;
1043 for(entry_iterator entry = iter->begin(); entry != iter->end(); ++entry) {
1049 maxlength = SeqOverlappingSchwarzAssembler<slu>
1050 ::assembleLocalProblems(rowToDomain, mat, solvers, subDomains, onTheFly);
1053 template<
class M,
class X,
class TM,
class TD,
class TA>
1058 : mat(mat_), solvers(sd.size()), subDomains(sd), relax(relaxationFactor),
1061 typedef typename subdomain_vector::const_iterator DomainIterator;
1063#ifdef DUNE_ISTL_WITH_CHECKING
1066 for(DomainIterator d=sd.
begin(); d != sd.
end(); ++d,++i) {
1068 assert(d->size()>0);
1069 typedef typename DomainIterator::value_type::const_iterator entry_iterator;
1071 for(entry_iterator entry = d->
begin(); entry != d->
end(); ++entry) {
1084 for(DomainIterator domain=sd.begin(); domain != sd.end(); ++domain, ++domainId) {
1085 typedef typename subdomain_type::const_iterator iterator;
1086 for(iterator row=domain->begin(); row != domain->end(); ++row)
1087 rowToDomain[*row].push_back(domainId);
1090 maxlength = SeqOverlappingSchwarzAssembler<slu>
1091 ::assembleLocalProblems(rowToDomain, mat, solvers, subDomains, onTheFly);
1103 template<
typename T,
typename A>
1106 static constexpr size_t n = std::decay_t<decltype(Impl::asMatrix(std::declval<T>()))>::rows;
1107 static constexpr size_t m = std::decay_t<decltype(Impl::asMatrix(std::declval<T>()))>::cols;
1108 template<
class Domain>
1109 static int size(
const Domain & d)
1116 template<
class K,
class Al,
class X,
class Y>
1117 template<
class RowToDomain,
class Solvers,
class SubDomains>
1119 SeqOverlappingSchwarzAssemblerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>::
1120 assembleLocalProblems([[maybe_unused]]
const RowToDomain& rowToDomain,
1121 [[maybe_unused]]
const matrix_type& mat,
1122 [[maybe_unused]] Solvers& solvers,
1123 const SubDomains& subDomains,
1124 [[maybe_unused]]
bool onTheFly)
1126 typedef typename SubDomains::const_iterator DomainIterator;
1127 std::size_t maxlength = 0;
1131 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end(); ++domain)
1132 maxlength=
std::max(maxlength, domain->size());
1138#if HAVE_SUPERLU || HAVE_SUITESPARSE_UMFPACK
1139 template<
template<
class>
class S,
typename T,
typename A>
1140 template<
class RowToDomain,
class Solvers,
class SubDomains>
1141 std::size_t SeqOverlappingSchwarzAssemblerHelper<S<BCRSMatrix<T,A>>,
true>::assembleLocalProblems(
const RowToDomain& rowToDomain,
1142 const matrix_type& mat,
1144 const SubDomains& subDomains,
1147 typedef typename S<BCRSMatrix<T,A>>::MatrixInitializer MatrixInitializer;
1148 typedef typename std::vector<MatrixInitializer>::iterator InitializerIterator;
1149 typedef typename SubDomains::const_iterator DomainIterator;
1150 typedef typename Solvers::iterator SolverIterator;
1151 std::size_t maxlength = 0;
1154 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end(); ++domain)
1155 maxlength=
std::max(maxlength, domain->size());
1156 maxlength*=Impl::asMatrix(*mat[0].begin()).N();
1159 DomainIterator domain=subDomains.begin();
1162 std::vector<MatrixInitializer> initializers(subDomains.size());
1164 SolverIterator solver=solvers.begin();
1165 for(InitializerIterator initializer=initializers.begin(); initializer!=initializers.end();
1166 ++initializer, ++solver, ++domain) {
1167 solver->getInternalMatrix().N_=SeqOverlappingSchwarzDomainSize<matrix_type>::size(*domain);
1168 solver->getInternalMatrix().M_=SeqOverlappingSchwarzDomainSize<matrix_type>::size(*domain);
1170 *initializer=MatrixInitializer(solver->getInternalMatrix());
1174 typedef OverlappingSchwarzInitializer<std::vector<MatrixInitializer>,
1175 RowToDomain, SubDomains> Initializer;
1177 Initializer initializer(initializers, rowToDomain, subDomains);
1178 copyToBCCSMatrix(initializer, mat);
1181 for(
auto&& s: solvers)
1183 for (SolverIterator solverIt = solvers.begin(); solverIt != solvers.end(); ++solverIt)
1185 assert(solverIt->getInternalMatrix().N() == solverIt->getInternalMatrix().M());
1186 maxlength =
std::max(maxlength, solverIt->getInternalMatrix().N());
1194 template<
class M,
class X,
class Y>
1195 template<
class RowToDomain,
class Solvers,
class SubDomains>
1196 std::size_t SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>::assembleLocalProblems([[maybe_unused]]
const RowToDomain& rowToDomain,
1197 const matrix_type& mat,
1199 const SubDomains& subDomains,
1202 typedef typename SubDomains::const_iterator DomainIterator;
1203 typedef typename Solvers::iterator SolverIterator;
1204 std::size_t maxlength = 0;
1207 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end(); ++domain)
1208 maxlength=
std::max(maxlength, domain->size());
1211 SolverIterator solver=solvers.begin();
1212 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end();
1213 ++domain, ++solver) {
1214 solver->setSubMatrix(mat, *domain);
1215 maxlength=
std::max(maxlength, domain->size());
1224 template<
class M,
class X,
class TM,
class TD,
class TA>
1230 template<
class M,
class X,
class TM,
class TD,
class TA>
1231 template<
bool forward>
1234 typedef slu_vector solver_vector;
1235 typedef typename IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::solver_iterator iterator;
1236 typedef typename IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::domain_iterator
1239 OverlappingAssigner<TD> assigner(maxlength, mat, b, x);
1247 Adder adder(v, x, assigner, relax);
1251 std::for_each(domain->begin(), domain->end(), assigner);
1252 assigner.resetIndexForNextDomain();
1256 sdsolver.setSubMatrix(mat, *domain);
1258 sdsolver.apply(assigner.lhs(), assigner.rhs());
1260 solver->apply(assigner.lhs(), assigner.rhs());
1265 std::for_each(domain->begin(), domain->end(), adder);
1266 assigner.resetIndexForNextDomain();
1271 assigner.deallocate();
1274 template<
class K,
class Al,
class X,
class Y>
1275 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1277 const X& b_, Y& x_) :
1284 maxlength_(maxlength)
1287 template<
class K,
class Al,
class X,
class Y>
1289 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1296 template<
class K,
class Al,
class X,
class Y>
1298 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1299 ::resetIndexForNextDomain()
1304 template<
class K,
class Al,
class X,
class Y>
1306 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1312 template<
class K,
class Al,
class X,
class Y>
1314 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1320 template<
class K,
class Al,
class X,
class Y>
1322 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1323 ::relaxResult(field_type relax)
1328 template<
class K,
class Al,
class X,
class Y>
1330 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1331 ::operator()(
const size_type& domainIndex)
1336 for(size_type j=0; j<n; ++j, ++i) {
1337 assert(i<maxlength_);
1338 rhs()[i]=(*b)[domainIndex][j];
1345 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1346 block_type tmp(0.0);
1347 (*col).mv((*x)[col.index()], tmp);
1349 for(size_type j=0; j<n; ++j, ++i) {
1350 assert(i<maxlength_);
1356 for(size_type j=0; j<n; ++j, ++i) {
1357 assert(i<maxlength_);
1358 rhs()[i]=Impl::asVector((*b)[domainIndex])[j];
1364 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1365 for(size_type k=0; k<n; ++k) {
1366 rhs()[i]-=Impl::asMatrix(*col)[j][k] * Impl::asVector((*x)[col.index()])[k];
1373 template<
class K,
class Al,
class X,
class Y>
1375 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1376 ::assignResult(block_type& res)
1379 for(size_type j=0; j<n; ++j, ++i) {
1380 assert(i<maxlength_);
1381 Impl::asVector(res)[j]+=lhs()[i];
1385#if HAVE_SUPERLU || HAVE_SUITESPARSE_UMFPACK
1387 template<
template<
class>
class S,
typename T,
typename A>
1388 OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>
1389 ::OverlappingAssignerHelper(std::size_t maxlength,
1391 const range_type& b_,
1395 x(&x_), i(0), maxlength_(maxlength)
1397 rhs_ =
new field_type[maxlength];
1398 lhs_ =
new field_type[maxlength];
1402 template<
template<
class>
class S,
typename T,
typename A>
1403 void OverlappingAssignerHelper<S<BCRSMatrix<T,A> >,
true>::deallocate()
1409 template<
template<
class>
class S,
typename T,
typename A>
1410 void OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::operator()(
const size_type& domainIndex)
1415 for(size_type j=0; j<n; ++j, ++i) {
1416 assert(i<maxlength_);
1417 rhs_[i]=Impl::asVector((*b)[domainIndex])[j];
1425 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1427 Impl::asMatrix(*col).mv((*x)[col.index()], tmp);
1429 for(size_type j=0; j<n; ++j, ++i) {
1430 assert(i<maxlength_);
1431 rhs_[i]-=Impl::asVector(tmp)[j];
1438 template<
template<
class>
class S,
typename T,
typename A>
1439 void OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::relaxResult(field_type relax)
1441 for(size_type j=i+n; i<j; ++i) {
1442 assert(i<maxlength_);
1448 template<
template<
class>
class S,
typename T,
typename A>
1449 void OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::assignResult(block_type& res)
1452 for(size_type j=0; j<n; ++j, ++i) {
1453 assert(i<maxlength_);
1454 Impl::asVector(res)[j]+=lhs_[i];
1458 template<
template<
class>
class S,
typename T,
typename A>
1459 void OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::resetIndexForNextDomain()
1464 template<
template<
class>
class S,
typename T,
typename A>
1465 typename OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::field_type*
1466 OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::lhs()
1471 template<
template<
class>
class S,
typename T,
typename A>
1472 typename OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::field_type*
1473 OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::rhs()
1480 template<
class M,
class X,
class Y>
1481 OverlappingAssignerILUBase<M,X,Y>::OverlappingAssignerILUBase(std::size_t maxlength,
1489 rhs_=
new Y(maxlength);
1490 lhs_ =
new X(maxlength);
1493 template<
class M,
class X,
class Y>
1494 void OverlappingAssignerILUBase<M,X,Y>::deallocate()
1500 template<
class M,
class X,
class Y>
1501 void OverlappingAssignerILUBase<M,X,Y>::operator()(
const size_type& domainIndex)
1503 (*rhs_)[i]=(*b)[domainIndex];
1506 typedef typename matrix_type::ConstColIterator col_iterator;
1509 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1510 Impl::asMatrix(*col).mmv((*x)[col.index()], (*rhs_)[i]);
1516 template<
class M,
class X,
class Y>
1517 void OverlappingAssignerILUBase<M,X,Y>::relaxResult(field_type relax)
1522 template<
class M,
class X,
class Y>
1523 void OverlappingAssignerILUBase<M,X,Y>::assignResult(block_type& res)
1528 template<
class M,
class X,
class Y>
1529 X& OverlappingAssignerILUBase<M,X,Y>::lhs()
1534 template<
class M,
class X,
class Y>
1535 Y& OverlappingAssignerILUBase<M,X,Y>::rhs()
1540 template<
class M,
class X,
class Y>
1541 void OverlappingAssignerILUBase<M,X,Y>::resetIndexForNextDomain()
1546 template<
typename S,
typename T,
typename A>
1549 OverlappingAssigner<S>& assigner_,
1550 const field_type& relax_)
1551 : v(&v_), x(&x_), assigner(&assigner_), relax(relax_)
1554 template<
typename S,
typename T,
typename A>
1555 void AdditiveAdder<S,BlockVector<T,A> >::operator()(
const size_type& domainIndex)
1558 assigner->assignResult((*v)[domainIndex]);
1562 template<
typename S,
typename T,
typename A>
1563 void AdditiveAdder<S,BlockVector<T,A> >::axpy()
1570 template<
typename S,
typename T,
typename A>
1571 MultiplicativeAdder<S,BlockVector<T,A> >
1572 ::MultiplicativeAdder([[maybe_unused]] BlockVector<T,A>& v_,
1573 BlockVector<T,A>& x_,
1574 OverlappingAssigner<S>& assigner_,
const field_type& relax_)
1575 : x(&x_), assigner(&assigner_), relax(relax_)
1579 template<
typename S,
typename T,
typename A>
1580 void MultiplicativeAdder<S,BlockVector<T,A> >::operator()(
const size_type& domainIndex)
1583 assigner->relaxResult(relax);
1584 assigner->assignResult((*x)[domainIndex]);
1588 template<
typename S,
typename T,
typename A>
1589 void MultiplicativeAdder<S,BlockVector<T,A> >::axpy()
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:464
row_type::ConstIterator ConstColIterator
Const iterator to the entries of a row.
Definition: bcrsmatrix.hh:739
Iterator begin()
begin iterator
Definition: densevector.hh:348
Iterator end()
end iterator
Definition: densevector.hh:354
Iterator find(size_type i)
return iterator to given element or end()
Definition: densevector.hh:374
Exact subdomain solver using Dune::DynamicMatrix<T>::solve.
Definition: overlappingschwarz.hh:138
size_type capacity() const
Number of elements for which memory has been allocated.
Definition: dynvector.hh:135
Index Set Interface base class.
Definition: indexidset.hh:76
Initializer for SuperLU Matrices representing the subdomains.
Definition: overlappingschwarz.hh:45
D subdomain_vector
The vector type containing the subdomain to row index mapping.
Definition: overlappingschwarz.hh:48
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:30
A constant iterator for the SLList.
Definition: sllist.hh:369
A single linked list.
Definition: sllist.hh:42
Sequential overlapping Schwarz preconditioner.
Definition: overlappingschwarz.hh:753
X::field_type field_type
The field type of the preconditioner.
Definition: overlappingschwarz.hh:781
void apply(X &v, const X &d)
Apply one step of the preconditioner to the system A(v)=d.
SLList< size_type, typename std::allocator_traits< TA >::template rebind_alloc< size_type > > subdomain_list
The type for the row to subdomain mapping.
Definition: overlappingschwarz.hh:798
std::vector< slu, typename std::allocator_traits< TA >::template rebind_alloc< slu > > slu_vector
The vector type containing subdomain solvers.
Definition: overlappingschwarz.hh:807
M matrix_type
The type of the matrix to precondition.
Definition: overlappingschwarz.hh:758
TM Mode
The mode (additive or multiplicative) of the Schwarz method.
Definition: overlappingschwarz.hh:776
X range_type
The range type of the preconditioner.
Definition: overlappingschwarz.hh:768
std::set< size_type, std::less< size_type >, typename std::allocator_traits< TA >::template rebind_alloc< size_type > > subdomain_type
The type for the subdomain to row index mapping.
Definition: overlappingschwarz.hh:792
X domain_type
The domain type of the preconditioner.
Definition: overlappingschwarz.hh:763
TD slu
The type for the subdomain solver in use.
Definition: overlappingschwarz.hh:804
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: overlappingschwarz.hh:866
virtual void pre(X &x, X &b)
Prepare the preconditioner.
Definition: overlappingschwarz.hh:844
virtual void post(X &x)
Postprocess the preconditioner.
Definition: overlappingschwarz.hh:859
TA allocator
The allocator to use.
Definition: overlappingschwarz.hh:787
std::vector< subdomain_type, typename std::allocator_traits< TA >::template rebind_alloc< subdomain_type > > subdomain_vector
The vector type containing the subdomain to row index mapping.
Definition: overlappingschwarz.hh:795
std::vector< subdomain_list, typename std::allocator_traits< TA >::template rebind_alloc< subdomain_list > > rowtodomain_vector
The vector type containing the row index to subdomain mapping.
Definition: overlappingschwarz.hh:801
matrix_type::size_type size_type
The return type of the size method.
Definition: overlappingschwarz.hh:784
This file implements a dense matrix with dynamic numbers of rows and columns.
virtual void apply(X &v, const X &d)
Apply the precondtioner.
Definition: overlappingschwarz.hh:1225
SeqOverlappingSchwarz(const matrix_type &mat, const subdomain_vector &subDomains, field_type relaxationFactor=1, bool onTheFly_=true)
Construct the overlapping Schwarz method.
Definition: overlappingschwarz.hh:1054
SeqOverlappingSchwarz(const matrix_type &mat, const rowtodomain_vector &rowToDomain, field_type relaxationFactor=1, bool onTheFly_=true)
Definition: overlappingschwarz.hh:1007
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
DVVerbType dvverb(std::cout)
stream for very verbose output.
Definition: stdstreams.hh:93
Various local subdomain solvers based on ILU for SeqOverlappingSchwarz.
Dune namespace.
Definition: alignedallocator.hh:11
Define general preconditioner interface.
Implements a singly linked list together with the necessary iterators.
Templates characterizing the type of a solver.
template meta program for choosing how to add the correction.
Definition: overlappingschwarz.hh:570
Tag that the tells the Schwarz method to be additive.
Definition: overlappingschwarz.hh:118
Helper template meta program for application of overlapping Schwarz.
Definition: overlappingschwarz.hh:603
Tag that tells the Schwarz method to be multiplicative.
Definition: overlappingschwarz.hh:124
Helper template meta program for application of overlapping Schwarz.
Definition: overlappingschwarz.hh:667
Definition: overlappingschwarz.hh:1101
Category
Definition: solvercategory.hh:21
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:23
Tag that tells the Schwarz method to be multiplicative and symmetric.
Definition: overlappingschwarz.hh:131
Classes for using SuperLU with ISTL matrices.
Classes for using UMFPack with ISTL matrices.