5#ifndef DUNE_ISTL_OVERLAPPINGSCHWARZ_HH
6#define DUNE_ISTL_OVERLAPPINGSCHWARZ_HH
16#include <dune/istl/bccsmatrixinitializer.hh>
39 template<
class M,
class X,
class TM,
class TD,
class TA>
40 class SeqOverlappingSchwarz;
45 template<
class I,
class S,
class D>
52 typedef I InitializerList;
53 typedef typename InitializerList::value_type AtomInitializer;
54 typedef typename AtomInitializer::Matrix Matrix;
55 typedef typename Matrix::const_iterator Iter;
56 typedef typename Matrix::row_type::const_iterator CIter;
59 typedef typename IndexSet::size_type size_type;
62 const IndexSet& indices,
66 void addRowNnz(
const Iter& row);
70 void countEntries(
const Iter& row,
const CIter& col)
const;
72 void calcColstart()
const;
74 void copyValue(
const Iter& row,
const CIter& col)
const;
76 void createMatrix()
const;
81 typedef typename S::size_type size_type;
82 typedef std::map<size_type,size_type> Map;
83 typedef typename Map::iterator iterator;
84 typedef typename Map::const_iterator const_iterator;
88 void insert(size_type grow);
90 const_iterator find(size_type grow)
const;
92 iterator find(size_type grow);
96 const_iterator begin()
const;
100 const_iterator end()
const;
103 std::map<size_type,size_type> map_;
108 typedef typename InitializerList::iterator InitIterator;
109 typedef typename IndexSet::const_iterator IndexIteratur;
110 InitializerList* initializers;
112 mutable std::vector<IndexMap> indexMaps;
139 template<
class M,
class X,
class Y>
143 template<
class K,
class Al,
class X,
class Y>
149 typedef typename std::remove_const<M>::type matrix_type;
150 typedef typename X::field_type field_type;
151 typedef typename std::remove_const<M>::type rilu_type;
153 typedef X domain_type;
155 typedef Y range_type;
156 static constexpr size_t n = std::decay_t<decltype(Impl::asMatrix(std::declval<K>()))>::rows;
164 assert(v.size() > 0);
165 assert(v.size() == d.size());
166 assert(A.rows() <= v.size());
167 assert(A.cols() <= v.size());
168 size_t sz = A.rows();
184 void setSubMatrix(
const M& BCRS, S& rowset)
186 size_t sz = rowset.size();
188 typedef typename S::const_iterator SIter;
190 for(SIter rowIdx = rowset.begin(), rowEnd=rowset.end();
191 rowIdx!= rowEnd; ++rowIdx, r++)
194 for(SIter colIdx = rowset.begin(), colEnd=rowset.end();
195 colIdx!= colEnd; ++colIdx, c++)
197 if (BCRS[*rowIdx].find(*colIdx) == BCRS[*rowIdx].end())
199 for (
size_t i=0; i<n; i++)
201 for (
size_t j=0; j<n; j++)
203 A[r*n+i][c*n+j] = Impl::asMatrix(BCRS[*rowIdx][*colIdx])[i][j];
213 template<
typename T,
bool tag>
214 class OverlappingAssignerHelper
218 using OverlappingAssigner = OverlappingAssignerHelper<T, Dune::StoresColumnCompressed<T>::value>;
221 template<
class K,
class Al,
class X,
class Y>
222 class OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix<K, Al>, X, Y >,false>
225 typedef BCRSMatrix< K, Al> matrix_type;
226 typedef typename X::field_type field_type;
227 typedef Y range_type;
228 typedef typename range_type::block_type block_type;
229 typedef typename matrix_type::size_type size_type;
230 static constexpr size_t n = std::decay_t<decltype(Impl::asMatrix(std::declval<K>()))>::rows;
238 OverlappingAssignerHelper(std::size_t maxlength,
const BCRSMatrix<K, Al>& mat_,
const X& b_, Y& x_);
250 void resetIndexForNextDomain();
257 DynamicVector<field_type> & lhs();
264 DynamicVector<field_type> & rhs();
271 void relaxResult(field_type relax);
277 void operator()(
const size_type& domainIndex);
287 void assignResult(block_type& res);
293 const matrix_type* mat;
296 DynamicVector<field_type> * rhs_;
299 DynamicVector<field_type> * lhs_;
307 std::size_t maxlength_;
310#if HAVE_SUPERLU || HAVE_SUITESPARSE_UMFPACK
311 template<
template<
class>
class S,
typename T,
typename A>
312 struct OverlappingAssignerHelper<S<BCRSMatrix<T, A>>, true>
314 typedef BCRSMatrix<T, A> matrix_type;
315 typedef typename S<BCRSMatrix<T, A>>::range_type range_type;
316 typedef typename range_type::field_type field_type;
317 typedef typename range_type::block_type block_type;
319 typedef typename matrix_type::size_type size_type;
321 static constexpr size_t n = std::decay_t<decltype(Impl::asMatrix(std::declval<T>()))>::rows;
322 static constexpr size_t m = std::decay_t<decltype(Impl::asMatrix(std::declval<T>()))>::cols;
330 OverlappingAssignerHelper(std::size_t maxlength,
const matrix_type& mat,
331 const range_type& b, range_type& x);
342 void resetIndexForNextDomain();
360 void relaxResult(field_type relax);
366 void operator()(
const size_type& domain);
375 void assignResult(block_type& res);
381 const matrix_type* mat;
393 std::size_t maxlength_;
398 template<
class M,
class X,
class Y>
399 class OverlappingAssignerILUBase
402 typedef M matrix_type;
404 typedef typename Y::field_type field_type;
406 typedef typename Y::block_type block_type;
408 typedef typename matrix_type::size_type size_type;
416 OverlappingAssignerILUBase(std::size_t maxlength,
const M& mat,
428 void resetIndexForNextDomain();
446 void relaxResult(field_type relax);
452 void operator()(
const size_type& domain);
461 void assignResult(block_type& res);
481 template<
class M,
class X,
class Y>
482 class OverlappingAssignerHelper<ILU0SubdomainSolver<M,X,Y>, false>
483 :
public OverlappingAssignerILUBase<M,X,Y>
493 OverlappingAssignerHelper(std::size_t maxlength,
const M& mat,
495 : OverlappingAssignerILUBase<M,X,Y>(maxlength, mat,b,x)
500 template<
class M,
class X,
class Y>
501 class OverlappingAssignerHelper<ILUNSubdomainSolver<M,X,Y>,false>
502 :
public OverlappingAssignerILUBase<M,X,Y>
512 OverlappingAssignerHelper(std::size_t maxlength,
const M& mat,
514 : OverlappingAssignerILUBase<M,X,Y>(maxlength, mat,b,x)
518 template<
typename S,
typename T>
522 template<
typename S,
typename T,
typename A>
523 struct AdditiveAdder<S, BlockVector<T,A> >
525 typedef typename A::size_type size_type;
526 typedef typename std::decay_t<decltype(Impl::asVector(std::declval<T>()))>::field_type field_type;
527 AdditiveAdder(BlockVector<T,A>& v, BlockVector<T,A>& x,
528 OverlappingAssigner<S>& assigner,
const field_type& relax_);
529 void operator()(
const size_type& domain);
531 static constexpr size_t n = std::decay_t<decltype(Impl::asVector(std::declval<T>()))>::dimension;
536 OverlappingAssigner<S>* assigner;
540 template<
typename S,
typename T>
541 struct MultiplicativeAdder
544 template<
typename S,
typename T,
typename A>
545 struct MultiplicativeAdder<S, BlockVector<T,A> >
547 typedef typename A::size_type size_type;
548 typedef typename std::decay_t<decltype(Impl::asVector(std::declval<T>()))>::field_type field_type;
549 MultiplicativeAdder(BlockVector<T,A>& v, BlockVector<T,A>& x,
550 OverlappingAssigner<S>& assigner_,
const field_type& relax_);
551 void operator()(
const size_type& domain);
553 static constexpr size_t n = std::decay_t<decltype(Impl::asVector(std::declval<T>()))>::dimension;
557 OverlappingAssigner<S>* assigner;
570 template<
typename T,
class X,
class S>
574 template<
class X,
class S>
577 typedef AdditiveAdder<S,X> Adder;
580 template<
class X,
class S>
581 struct AdderSelector<MultiplicativeSchwarzMode,X,S>
583 typedef MultiplicativeAdder<S,X> Adder;
586 template<
class X,
class S>
587 struct AdderSelector<SymmetricMultiplicativeSchwarzMode,X,S>
589 typedef MultiplicativeAdder<S,X> Adder;
603 template<
typename T1,
typename T2,
bool forward>
606 typedef T1 solver_vector;
607 typedef typename solver_vector::iterator solver_iterator;
608 typedef T2 subdomain_vector;
609 typedef typename subdomain_vector::const_iterator domain_iterator;
611 static solver_iterator begin(solver_vector& sv)
616 static solver_iterator end(solver_vector& sv)
620 static domain_iterator begin(
const subdomain_vector& sv)
625 static domain_iterator end(
const subdomain_vector& sv)
631 template<
typename T1,
typename T2>
634 typedef T1 solver_vector;
635 typedef typename solver_vector::reverse_iterator solver_iterator;
636 typedef T2 subdomain_vector;
637 typedef typename subdomain_vector::const_reverse_iterator domain_iterator;
639 static solver_iterator begin(solver_vector& sv)
644 static solver_iterator end(solver_vector& sv)
648 static domain_iterator begin(
const subdomain_vector& sv)
653 static domain_iterator end(
const subdomain_vector& sv)
671 typedef typename smoother::range_type range_type;
673 static void apply(smoother& sm, range_type& v,
const range_type& b)
675 sm.template apply<true>(v, b);
679 template<
class M,
class X,
class TD,
class TA>
683 typedef typename smoother::range_type range_type;
685 static void apply(smoother& sm, range_type& v,
const range_type& b)
687 sm.template apply<true>(v, b);
688 sm.template apply<false>(v, b);
692 template<
class T,
bool tag>
693 struct SeqOverlappingSchwarzAssemblerHelper
697 using SeqOverlappingSchwarzAssembler = SeqOverlappingSchwarzAssemblerHelper<T,Dune::StoresColumnCompressed<T>::value>;
699 template<
class K,
class Al,
class X,
class Y>
700 struct SeqOverlappingSchwarzAssemblerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,false>
702 typedef BCRSMatrix< K, Al> matrix_type;
703 static constexpr size_t n = std::decay_t<decltype(Impl::asMatrix(std::declval<K>()))>::rows;
704 template<
class RowToDomain,
class Solvers,
class SubDomains>
705 static std::size_t assembleLocalProblems(
const RowToDomain& rowToDomain,
const matrix_type& mat,
706 Solvers& solvers,
const SubDomains& domains,
710 template<
template<
class>
class S,
typename T,
typename A>
711 struct SeqOverlappingSchwarzAssemblerHelper<S<BCRSMatrix<T,A>>,true>
713 typedef BCRSMatrix<T,A> matrix_type;
714 static constexpr size_t n = std::decay_t<decltype(Impl::asMatrix(std::declval<T>()))>::rows;
715 template<
class RowToDomain,
class Solvers,
class SubDomains>
716 static std::size_t assembleLocalProblems(
const RowToDomain& rowToDomain,
const matrix_type& mat,
717 Solvers& solvers,
const SubDomains& domains,
721 template<
class M,
class X,
class Y>
722 struct SeqOverlappingSchwarzAssemblerILUBase
724 typedef M matrix_type;
725 template<
class RowToDomain,
class Solvers,
class SubDomains>
726 static std::size_t assembleLocalProblems(
const RowToDomain& rowToDomain,
const matrix_type& mat,
727 Solvers& solvers,
const SubDomains& domains,
731 template<
class M,
class X,
class Y>
732 struct SeqOverlappingSchwarzAssemblerHelper<ILU0SubdomainSolver<M,X,Y>,false>
733 :
public SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>
736 template<
class M,
class X,
class Y>
737 struct SeqOverlappingSchwarzAssemblerHelper<ILUNSubdomainSolver<M,X,Y>,false>
738 :
public SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>
751 template<
class M,
class X,
class TM=AdditiveSchwarzMode,
752 class TD=ILU0SubdomainSolver<M,X,X>,
class TA=std::allocator<X> >
792 typedef std::set<size_type, std::less<size_type>,
793 typename std::allocator_traits<TA>::template rebind_alloc<size_type> >
797 typedef std::vector<subdomain_type, typename std::allocator_traits<TA>::template rebind_alloc<subdomain_type> >
subdomain_vector;
803 typedef std::vector<subdomain_list, typename std::allocator_traits<TA>::template rebind_alloc<subdomain_list> >
rowtodomain_vector;
809 typedef std::vector<slu, typename std::allocator_traits<TA>::template rebind_alloc<slu> >
slu_vector;
825 field_type relaxationFactor=1,
bool onTheFly_=
true);
839 field_type relaxationFactor=1,
bool onTheFly_=
true);
846 void pre ([[maybe_unused]] X& x, [[maybe_unused]] X& b)
override
854 void apply (X& v,
const X& d)
override;
861 void post ([[maybe_unused]] X& x)
override
864 template<
bool forward>
879 typename M::size_type maxlength;
886 template<
class I,
class S,
class D>
887 OverlappingSchwarzInitializer<I,S,D>::OverlappingSchwarzInitializer(InitializerList& il,
889 const subdomain_vector& domains_)
890 : initializers(&il), indices(&idx), indexMaps(il.
size()), domains(domains_)
894 template<
class I,
class S,
class D>
895 void OverlappingSchwarzInitializer<I,S,D>::addRowNnz(
const Iter& row)
897 typedef typename IndexSet::value_type::const_iterator iterator;
898 for(iterator domain=(*indices)[row.index()].begin(); domain != (*indices)[row.index()].end(); ++domain) {
899 (*initializers)[*domain].addRowNnz(row, domains[*domain]);
900 indexMaps[*domain].insert(row.index());
904 template<
class I,
class S,
class D>
905 void OverlappingSchwarzInitializer<I,S,D>::allocate()
907 for(
auto&& i: *initializers)
908 i.allocateMatrixStorage();
909 for(
auto&& i: *initializers)
913 template<
class I,
class S,
class D>
914 void OverlappingSchwarzInitializer<I,S,D>::countEntries(
const Iter& row,
const CIter& col)
const
916 typedef typename IndexSet::value_type::const_iterator iterator;
917 for(iterator domain=(*indices)[row.index()].begin(); domain != (*indices)[row.index()].end(); ++domain) {
918 typename std::map<size_type,size_type>::const_iterator v = indexMaps[*domain].
find(col.index());
919 if(v!= indexMaps[*domain].end()) {
920 (*initializers)[*domain].countEntries(indexMaps[*domain].find(col.index())->second);
925 template<
class I,
class S,
class D>
926 void OverlappingSchwarzInitializer<I,S,D>::calcColstart()
const
928 for(
auto&& i : *initializers)
932 template<
class I,
class S,
class D>
933 void OverlappingSchwarzInitializer<I,S,D>::copyValue(
const Iter& row,
const CIter& col)
const
935 typedef typename IndexSet::value_type::const_iterator iterator;
936 for(iterator domain=(*indices)[row.index()].begin(); domain!= (*indices)[row.index()].end(); ++domain) {
937 typename std::map<size_type,size_type>::const_iterator v = indexMaps[*domain].find(col.index());
938 if(v!= indexMaps[*domain].end()) {
939 assert(indexMaps[*domain].end()!=indexMaps[*domain].find(row.index()));
940 (*initializers)[*domain].copyValue(col, indexMaps[*domain].find(row.index())->second,
946 template<
class I,
class S,
class D>
947 void OverlappingSchwarzInitializer<I,S,D>::createMatrix()
const
949 std::vector<IndexMap>().swap(indexMaps);
950 for(
auto&& i: *initializers)
954 template<
class I,
class S,
class D>
955 OverlappingSchwarzInitializer<I,S,D>::IndexMap::IndexMap()
959 template<
class I,
class S,
class D>
960 void OverlappingSchwarzInitializer<I,S,D>::IndexMap::insert(size_type grow)
962 assert(map_.find(grow)==map_.end());
963 map_.insert(std::make_pair(grow, row++));
966 template<
class I,
class S,
class D>
967 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator
968 OverlappingSchwarzInitializer<I,S,D>::IndexMap::find(size_type grow)
const
970 return map_.find(grow);
973 template<
class I,
class S,
class D>
974 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator
975 OverlappingSchwarzInitializer<I,S,D>::IndexMap::find(size_type grow)
977 return map_.find(grow);
980 template<
class I,
class S,
class D>
981 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator
982 OverlappingSchwarzInitializer<I,S,D>::IndexMap::end()
const
987 template<
class I,
class S,
class D>
988 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator
989 OverlappingSchwarzInitializer<I,S,D>::IndexMap::end()
994 template<
class I,
class S,
class D>
995 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator
996 OverlappingSchwarzInitializer<I,S,D>::IndexMap::begin()
const
1001 template<
class I,
class S,
class D>
1002 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator
1003 OverlappingSchwarzInitializer<I,S,D>::IndexMap::begin()
1005 return map_.begin();
1008 template<
class M,
class X,
class TM,
class TD,
class TA>
1011 : mat(mat_), relax(relaxationFactor), onTheFly(fly)
1013 typedef typename rowtodomain_vector::const_iterator RowDomainIterator;
1015#ifdef DUNE_ISTL_WITH_CHECKING
1016 assert(rowToDomain.size()==mat.N());
1017 assert(rowToDomain.size()==mat.M());
1019 for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter)
1020 assert(iter->size()>0);
1025 for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter)
1026 for(DomainIterator d=iter->
begin(); d != iter->
end(); ++d)
1027 domains=std::max(domains, *d);
1030 solvers.resize(domains);
1031 subDomains.resize(domains);
1035 for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter, ++row)
1036 for(DomainIterator d=iter->
begin(); d != iter->
end(); ++d)
1037 subDomains[*d].insert(row);
1039#ifdef DUNE_ISTL_WITH_CHECKING
1041 typedef typename subdomain_vector::const_iterator iterator;
1042 for(iterator iter=subDomains.begin(); iter != subDomains.end(); ++iter) {
1043 typedef typename subdomain_type::const_iterator entry_iterator;
1045 for(entry_iterator entry = iter->begin(); entry != iter->end(); ++entry) {
1051 maxlength = SeqOverlappingSchwarzAssembler<slu>
1052 ::assembleLocalProblems(rowToDomain, mat, solvers, subDomains, onTheFly);
1055 template<
class M,
class X,
class TM,
class TD,
class TA>
1060 : mat(mat_), solvers(sd.
size()), subDomains(sd), relax(relaxationFactor),
1063 typedef typename subdomain_vector::const_iterator DomainIterator;
1065#ifdef DUNE_ISTL_WITH_CHECKING
1068 for(DomainIterator d=sd.
begin(); d != sd.
end(); ++d,++i) {
1070 assert(d->size()>0);
1071 typedef typename DomainIterator::value_type::const_iterator entry_iterator;
1073 for(entry_iterator entry = d->
begin(); entry != d->
end(); ++entry) {
1086 for(DomainIterator domain=sd.begin(); domain != sd.end(); ++domain, ++domainId) {
1087 typedef typename subdomain_type::const_iterator iterator;
1088 for(iterator row=domain->begin(); row != domain->end(); ++row)
1089 rowToDomain[*row].push_back(domainId);
1092 maxlength = SeqOverlappingSchwarzAssembler<slu>
1093 ::assembleLocalProblems(rowToDomain, mat, solvers, subDomains, onTheFly);
1105 template<
typename T,
typename A>
1108 static constexpr size_t n = std::decay_t<decltype(Impl::asMatrix(std::declval<T>()))>::rows;
1109 static constexpr size_t m = std::decay_t<decltype(Impl::asMatrix(std::declval<T>()))>::cols;
1110 template<
class Domain>
1111 static int size(
const Domain & d)
1118 template<
class K,
class Al,
class X,
class Y>
1119 template<
class RowToDomain,
class Solvers,
class SubDomains>
1121 SeqOverlappingSchwarzAssemblerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>::
1122 assembleLocalProblems([[maybe_unused]]
const RowToDomain& rowToDomain,
1123 [[maybe_unused]]
const matrix_type& mat,
1124 [[maybe_unused]] Solvers& solvers,
1125 const SubDomains& subDomains,
1126 [[maybe_unused]]
bool onTheFly)
1128 typedef typename SubDomains::const_iterator DomainIterator;
1129 std::size_t maxlength = 0;
1133 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end(); ++domain)
1134 maxlength=std::max(maxlength, domain->size());
1140#if HAVE_SUPERLU || HAVE_SUITESPARSE_UMFPACK
1141 template<
template<
class>
class S,
typename T,
typename A>
1142 template<
class RowToDomain,
class Solvers,
class SubDomains>
1143 std::size_t SeqOverlappingSchwarzAssemblerHelper<S<BCRSMatrix<T,A>>,
true>::assembleLocalProblems(
const RowToDomain& rowToDomain,
1144 const matrix_type& mat,
1146 const SubDomains& subDomains,
1149 typedef typename S<BCRSMatrix<T,A>>::MatrixInitializer MatrixInitializer;
1150 typedef typename std::vector<MatrixInitializer>::iterator InitializerIterator;
1151 typedef typename SubDomains::const_iterator DomainIterator;
1152 typedef typename Solvers::iterator SolverIterator;
1153 std::size_t maxlength = 0;
1156 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end(); ++domain)
1157 maxlength=std::max(maxlength, domain->size());
1158 maxlength*=Impl::asMatrix(*mat[0].begin()).N();
1161 DomainIterator domain=subDomains.begin();
1164 std::vector<MatrixInitializer> initializers(subDomains.size());
1166 SolverIterator solver=solvers.begin();
1167 for(InitializerIterator initializer=initializers.begin(); initializer!=initializers.end();
1168 ++initializer, ++solver, ++domain) {
1169 solver->getInternalMatrix().N_=SeqOverlappingSchwarzDomainSize<matrix_type>::size(*domain);
1170 solver->getInternalMatrix().M_=SeqOverlappingSchwarzDomainSize<matrix_type>::size(*domain);
1172 *initializer=MatrixInitializer(solver->getInternalMatrix());
1176 typedef OverlappingSchwarzInitializer<std::vector<MatrixInitializer>,
1177 RowToDomain, SubDomains> Initializer;
1179 Initializer initializer(initializers, rowToDomain, subDomains);
1180 copyToBCCSMatrix(initializer, mat);
1183 for(
auto&& s: solvers)
1185 for (SolverIterator solverIt = solvers.begin(); solverIt != solvers.end(); ++solverIt)
1187 assert(solverIt->getInternalMatrix().N() == solverIt->getInternalMatrix().M());
1188 maxlength = std::max(maxlength, solverIt->getInternalMatrix().N());
1196 template<
class M,
class X,
class Y>
1197 template<
class RowToDomain,
class Solvers,
class SubDomains>
1198 std::size_t SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>::assembleLocalProblems([[maybe_unused]]
const RowToDomain& rowToDomain,
1199 const matrix_type& mat,
1201 const SubDomains& subDomains,
1204 typedef typename SubDomains::const_iterator DomainIterator;
1205 typedef typename Solvers::iterator SolverIterator;
1206 std::size_t maxlength = 0;
1209 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end(); ++domain)
1210 maxlength=std::max(maxlength, domain->size());
1213 SolverIterator solver=solvers.begin();
1214 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end();
1215 ++domain, ++solver) {
1216 solver->setSubMatrix(mat, *domain);
1217 maxlength=std::max(maxlength, domain->size());
1226 template<
class M,
class X,
class TM,
class TD,
class TA>
1232 template<
class M,
class X,
class TM,
class TD,
class TA>
1233 template<
bool forward>
1236 typedef slu_vector solver_vector;
1237 typedef typename IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::solver_iterator iterator;
1238 typedef typename IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::domain_iterator
1241 OverlappingAssigner<TD> assigner(maxlength, mat, b, x);
1249 Adder adder(v, x, assigner, relax);
1253 std::for_each(domain->begin(), domain->end(), assigner);
1254 assigner.resetIndexForNextDomain();
1258 sdsolver.setSubMatrix(mat, *domain);
1260 sdsolver.apply(assigner.lhs(), assigner.rhs());
1262 solver->apply(assigner.lhs(), assigner.rhs());
1267 std::for_each(domain->begin(), domain->end(), adder);
1268 assigner.resetIndexForNextDomain();
1273 assigner.deallocate();
1276 template<
class K,
class Al,
class X,
class Y>
1277 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1279 const X& b_, Y& x_) :
1286 maxlength_(maxlength)
1289 template<
class K,
class Al,
class X,
class Y>
1291 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1298 template<
class K,
class Al,
class X,
class Y>
1300 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1301 ::resetIndexForNextDomain()
1306 template<
class K,
class Al,
class X,
class Y>
1308 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1314 template<
class K,
class Al,
class X,
class Y>
1316 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1322 template<
class K,
class Al,
class X,
class Y>
1324 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1325 ::relaxResult(field_type relax)
1330 template<
class K,
class Al,
class X,
class Y>
1332 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1333 ::operator()(
const size_type& domainIndex)
1338 for(size_type j=0; j<n; ++j, ++i) {
1339 assert(i<maxlength_);
1340 rhs()[i]=(*b)[domainIndex][j];
1347 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1348 block_type tmp(0.0);
1349 (*col).mv((*x)[col.index()], tmp);
1351 for(size_type j=0; j<n; ++j, ++i) {
1352 assert(i<maxlength_);
1358 for(size_type j=0; j<n; ++j, ++i) {
1359 assert(i<maxlength_);
1360 rhs()[i]=Impl::asVector((*b)[domainIndex])[j];
1366 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1367 for(size_type k=0; k<n; ++k) {
1368 rhs()[i]-=Impl::asMatrix(*col)[j][k] * Impl::asVector((*x)[col.index()])[k];
1375 template<
class K,
class Al,
class X,
class Y>
1377 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1378 ::assignResult(block_type& res)
1381 for(size_type j=0; j<n; ++j, ++i) {
1382 assert(i<maxlength_);
1383 Impl::asVector(res)[j]+=lhs()[i];
1387#if HAVE_SUPERLU || HAVE_SUITESPARSE_UMFPACK
1389 template<
template<
class>
class S,
typename T,
typename A>
1390 OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>
1391 ::OverlappingAssignerHelper(std::size_t maxlength,
1393 const range_type& b_,
1397 x(&x_), i(0), maxlength_(maxlength)
1399 rhs_ =
new field_type[maxlength];
1400 lhs_ =
new field_type[maxlength];
1404 template<
template<
class>
class S,
typename T,
typename A>
1405 void OverlappingAssignerHelper<S<BCRSMatrix<T,A> >,
true>::deallocate()
1411 template<
template<
class>
class S,
typename T,
typename A>
1412 void OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::operator()(
const size_type& domainIndex)
1417 for(size_type j=0; j<n; ++j, ++i) {
1418 assert(i<maxlength_);
1419 rhs_[i]=Impl::asVector((*b)[domainIndex])[j];
1427 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1429 Impl::asMatrix(*col).mv((*x)[col.index()], tmp);
1431 for(size_type j=0; j<n; ++j, ++i) {
1432 assert(i<maxlength_);
1433 rhs_[i]-=Impl::asVector(tmp)[j];
1440 template<
template<
class>
class S,
typename T,
typename A>
1441 void OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::relaxResult(field_type relax)
1443 for(size_type j=i+n; i<j; ++i) {
1444 assert(i<maxlength_);
1450 template<
template<
class>
class S,
typename T,
typename A>
1451 void OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::assignResult(block_type& res)
1454 for(size_type j=0; j<n; ++j, ++i) {
1455 assert(i<maxlength_);
1456 Impl::asVector(res)[j]+=lhs_[i];
1460 template<
template<
class>
class S,
typename T,
typename A>
1461 void OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::resetIndexForNextDomain()
1466 template<
template<
class>
class S,
typename T,
typename A>
1467 typename OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::field_type*
1468 OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::lhs()
1473 template<
template<
class>
class S,
typename T,
typename A>
1474 typename OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::field_type*
1475 OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::rhs()
1482 template<
class M,
class X,
class Y>
1483 OverlappingAssignerILUBase<M,X,Y>::OverlappingAssignerILUBase(std::size_t maxlength,
1491 rhs_=
new Y(maxlength);
1492 lhs_ =
new X(maxlength);
1495 template<
class M,
class X,
class Y>
1496 void OverlappingAssignerILUBase<M,X,Y>::deallocate()
1502 template<
class M,
class X,
class Y>
1503 void OverlappingAssignerILUBase<M,X,Y>::operator()(
const size_type& domainIndex)
1505 (*rhs_)[i]=(*b)[domainIndex];
1508 typedef typename matrix_type::ConstColIterator col_iterator;
1511 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1512 Impl::asMatrix(*col).mmv((*x)[col.index()], (*rhs_)[i]);
1518 template<
class M,
class X,
class Y>
1519 void OverlappingAssignerILUBase<M,X,Y>::relaxResult(field_type relax)
1524 template<
class M,
class X,
class Y>
1525 void OverlappingAssignerILUBase<M,X,Y>::assignResult(block_type& res)
1530 template<
class M,
class X,
class Y>
1531 X& OverlappingAssignerILUBase<M,X,Y>::lhs()
1536 template<
class M,
class X,
class Y>
1537 Y& OverlappingAssignerILUBase<M,X,Y>::rhs()
1542 template<
class M,
class X,
class Y>
1543 void OverlappingAssignerILUBase<M,X,Y>::resetIndexForNextDomain()
1548 template<
typename S,
typename T,
typename A>
1551 OverlappingAssigner<S>& assigner_,
1552 const field_type& relax_)
1553 : v(&v_), x(&x_), assigner(&assigner_), relax(relax_)
1556 template<
typename S,
typename T,
typename A>
1557 void AdditiveAdder<S,BlockVector<T,A> >::operator()(
const size_type& domainIndex)
1560 assigner->assignResult((*v)[domainIndex]);
1564 template<
typename S,
typename T,
typename A>
1565 void AdditiveAdder<S,BlockVector<T,A> >::axpy()
1572 template<
typename S,
typename T,
typename A>
1573 MultiplicativeAdder<S,BlockVector<T,A> >
1574 ::MultiplicativeAdder([[maybe_unused]] BlockVector<T,A>& v_,
1575 BlockVector<T,A>& x_,
1576 OverlappingAssigner<S>& assigner_,
const field_type& relax_)
1577 : x(&x_), assigner(&assigner_), relax(relax_)
1581 template<
typename S,
typename T,
typename A>
1582 void MultiplicativeAdder<S,BlockVector<T,A> >::operator()(
const size_type& domainIndex)
1585 assigner->relaxResult(relax);
1586 assigner->assignResult((*x)[domainIndex]);
1590 template<
typename S,
typename T,
typename A>
1591 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:466
row_type::ConstIterator ConstColIterator
Const iterator to the entries of a row.
Definition: bcrsmatrix.hh:737
Iterator begin()
begin iterator
Definition: densevector.hh:347
Iterator end()
end iterator
Definition: densevector.hh:353
Iterator find(size_type i)
return iterator to given element or end()
Definition: densevector.hh:373
Exact subdomain solver using Dune::DynamicMatrix<T>::solve.
Definition: overlappingschwarz.hh:140
size_type capacity() const
Number of elements for which memory has been allocated.
Definition: dynvector.hh:137
Index Set Interface base class.
Definition: indexidset.hh:78
Initializer for SuperLU Matrices representing the subdomains.
Definition: overlappingschwarz.hh:47
D subdomain_vector
The vector type containing the subdomain to row index mapping.
Definition: overlappingschwarz.hh:50
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:33
A constant iterator for the SLList.
Definition: sllist.hh:371
A single linked list.
Definition: sllist.hh:44
Sequential overlapping Schwarz preconditioner.
Definition: overlappingschwarz.hh:755
X::field_type field_type
The field type of the preconditioner.
Definition: overlappingschwarz.hh:783
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:800
std::vector< slu, typename std::allocator_traits< TA >::template rebind_alloc< slu > > slu_vector
The vector type containing subdomain solvers.
Definition: overlappingschwarz.hh:809
SolverCategory::Category category() const override
Category of the preconditioner (see SolverCategory::Category)
Definition: overlappingschwarz.hh:868
M matrix_type
The type of the matrix to precondition.
Definition: overlappingschwarz.hh:760
TM Mode
The mode (additive or multiplicative) of the Schwarz method.
Definition: overlappingschwarz.hh:778
void post(X &x) override
Postprocess the preconditioner.
Definition: overlappingschwarz.hh:861
X range_type
The range type of the preconditioner.
Definition: overlappingschwarz.hh:770
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:794
X domain_type
The domain type of the preconditioner.
Definition: overlappingschwarz.hh:765
TD slu
The type for the subdomain solver in use.
Definition: overlappingschwarz.hh:806
TA allocator
The allocator to use.
Definition: overlappingschwarz.hh:789
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:797
void pre(X &x, X &b) override
Prepare the preconditioner.
Definition: overlappingschwarz.hh:846
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:803
matrix_type::size_type size_type
The return type of the size method.
Definition: overlappingschwarz.hh:786
Classes for using UMFPack with ISTL matrices.
This file implements a dense matrix with dynamic numbers of rows and columns.
SeqOverlappingSchwarz(const matrix_type &mat, const subdomain_vector &subDomains, field_type relaxationFactor=1, bool onTheFly_=true)
Construct the overlapping Schwarz method.
Definition: overlappingschwarz.hh:1056
void apply(X &v, const X &d) override
Apply the preconditioner.
Definition: overlappingschwarz.hh:1227
SeqOverlappingSchwarz(const matrix_type &mat, const rowtodomain_vector &rowToDomain, field_type relaxationFactor=1, bool onTheFly_=true)
Definition: overlappingschwarz.hh:1009
DVVerbType dvverb(std::cout)
stream for very verbose output.
Definition: stdstreams.hh:96
Various local subdomain solvers based on ILU for SeqOverlappingSchwarz.
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
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:572
Tag that the tells the Schwarz method to be additive.
Definition: overlappingschwarz.hh:120
Helper template meta program for application of overlapping Schwarz.
Definition: overlappingschwarz.hh:605
Tag that tells the Schwarz method to be multiplicative.
Definition: overlappingschwarz.hh:126
Helper template meta program for application of overlapping Schwarz.
Definition: overlappingschwarz.hh:669
Definition: overlappingschwarz.hh:1103
Category
Definition: solvercategory.hh:23
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:25
Tag that tells the Schwarz method to be multiplicative and symmetric.
Definition: overlappingschwarz.hh:133
Classes for using SuperLU with ISTL matrices.