3#ifndef DUNE_ISTL_OVERLAPPINGSCHWARZ_HH
4#define DUNE_ISTL_OVERLAPPINGSCHWARZ_HH
35 template<
class M,
class X,
class TM,
class TD,
class TA>
36 class SeqOverlappingSchwarz;
41 template<
class I,
class S,
class D>
48 typedef I InitializerList;
49 typedef typename InitializerList::value_type AtomInitializer;
50 typedef typename AtomInitializer::Matrix Matrix;
51 typedef typename Matrix::const_iterator Iter;
52 typedef typename Matrix::row_type::const_iterator CIter;
55 typedef typename IndexSet::size_type size_type;
58 const IndexSet& indices,
62 void addRowNnz(
const Iter& row);
66 void countEntries(
const Iter& row,
const CIter& col)
const;
68 void calcColstart()
const;
70 void copyValue(
const Iter& row,
const CIter& col)
const;
72 void createMatrix()
const;
77 typedef typename S::size_type size_type;
78 typedef std::map<size_type,size_type> Map;
79 typedef typename Map::iterator iterator;
80 typedef typename Map::const_iterator const_iterator;
84 void insert(size_type grow);
86 const_iterator find(size_type grow)
const;
88 iterator find(size_type grow);
92 const_iterator begin()
const;
96 const_iterator end()
const;
99 std::map<size_type,size_type> map_;
104 typedef typename InitializerList::iterator InitIterator;
105 typedef typename IndexSet::const_iterator IndexIteratur;
106 InitializerList* initializers;
108 mutable std::vector<IndexMap> indexMaps;
135 template<
class M,
class X,
class Y>
139 template<
class K,
class Al,
class X,
class Y>
145 typedef typename std::remove_const<M>::type matrix_type;
146 typedef typename X::field_type field_type;
147 typedef typename std::remove_const<M>::type rilu_type;
149 typedef X domain_type;
151 typedef Y range_type;
152 static constexpr size_t n = std::decay_t<decltype(Impl::asMatrix(std::declval<K>()))>::rows;
160 assert(v.size() > 0);
161 assert(v.size() == d.size());
162 assert(A.rows() <= v.size());
163 assert(A.cols() <= v.size());
164 size_t sz = A.rows();
180 void setSubMatrix(
const M& BCRS, S& rowset)
182 size_t sz = rowset.size();
184 typedef typename S::const_iterator SIter;
186 for(SIter rowIdx = rowset.begin(), rowEnd=rowset.end();
187 rowIdx!= rowEnd; ++rowIdx, r++)
190 for(SIter colIdx = rowset.begin(), colEnd=rowset.end();
191 colIdx!= colEnd; ++colIdx, c++)
193 if (BCRS[*rowIdx].find(*colIdx) == BCRS[*rowIdx].end())
195 for (
size_t i=0; i<n; i++)
197 for (
size_t j=0; j<n; j++)
199 A[r*n+i][c*n+j] = Impl::asMatrix(BCRS[*rowIdx][*colIdx])[i][j];
209 template<
typename T,
bool tag>
210 class OverlappingAssignerHelper
214 using OverlappingAssigner = OverlappingAssignerHelper<T, Dune::StoresColumnCompressed<T>::value>;
217 template<
class K,
class Al,
class X,
class Y>
218 class OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix<K, Al>, X, Y >,false>
221 typedef BCRSMatrix< K, Al> matrix_type;
222 typedef typename X::field_type field_type;
223 typedef Y range_type;
224 typedef typename range_type::block_type block_type;
225 typedef typename matrix_type::size_type size_type;
226 static constexpr size_t n = std::decay_t<decltype(Impl::asMatrix(std::declval<K>()))>::rows;
234 OverlappingAssignerHelper(std::size_t maxlength,
const BCRSMatrix<K, Al>& mat_,
const X& b_, Y& x_);
246 void resetIndexForNextDomain();
253 DynamicVector<field_type> & lhs();
260 DynamicVector<field_type> & rhs();
267 void relaxResult(field_type relax);
273 void operator()(
const size_type& domainIndex);
283 void assignResult(block_type& res);
289 const matrix_type* mat;
292 DynamicVector<field_type> * rhs_;
295 DynamicVector<field_type> * lhs_;
303 std::size_t maxlength_;
306#if HAVE_SUPERLU || HAVE_SUITESPARSE_UMFPACK
307 template<
template<
class>
class S,
typename T,
typename A>
308 struct OverlappingAssignerHelper<S<BCRSMatrix<T, A>>, true>
310 typedef BCRSMatrix<T, A> matrix_type;
311 typedef typename S<BCRSMatrix<T, A>>::range_type range_type;
312 typedef typename range_type::field_type field_type;
313 typedef typename range_type::block_type block_type;
315 typedef typename matrix_type::size_type size_type;
317 static constexpr size_t n = std::decay_t<decltype(Impl::asMatrix(std::declval<T>()))>::rows;
318 static constexpr size_t m = std::decay_t<decltype(Impl::asMatrix(std::declval<T>()))>::cols;
326 OverlappingAssignerHelper(std::size_t maxlength,
const matrix_type& mat,
327 const range_type& b, range_type& x);
338 void resetIndexForNextDomain();
356 void relaxResult(field_type relax);
362 void operator()(
const size_type& domain);
371 void assignResult(block_type& res);
377 const matrix_type* mat;
389 std::size_t maxlength_;
394 template<
class M,
class X,
class Y>
395 class OverlappingAssignerILUBase
398 typedef M matrix_type;
400 typedef typename Y::field_type field_type;
402 typedef typename Y::block_type block_type;
404 typedef typename matrix_type::size_type size_type;
412 OverlappingAssignerILUBase(std::size_t maxlength,
const M& mat,
424 void resetIndexForNextDomain();
442 void relaxResult(field_type relax);
448 void operator()(
const size_type& domain);
457 void assignResult(block_type& res);
477 template<
class M,
class X,
class Y>
478 class OverlappingAssignerHelper<ILU0SubdomainSolver<M,X,Y>, false>
479 :
public OverlappingAssignerILUBase<M,X,Y>
489 OverlappingAssignerHelper(std::size_t maxlength,
const M& mat,
491 : OverlappingAssignerILUBase<M,X,Y>(maxlength, mat,b,x)
496 template<
class M,
class X,
class Y>
497 class OverlappingAssignerHelper<ILUNSubdomainSolver<M,X,Y>,false>
498 :
public OverlappingAssignerILUBase<M,X,Y>
508 OverlappingAssignerHelper(std::size_t maxlength,
const M& mat,
510 : OverlappingAssignerILUBase<M,X,Y>(maxlength, mat,b,x)
514 template<
typename S,
typename T>
518 template<
typename S,
typename T,
typename A>
519 struct AdditiveAdder<S, BlockVector<T,A> >
521 typedef typename A::size_type size_type;
522 typedef typename std::decay_t<decltype(Impl::asVector(std::declval<T>()))>::field_type field_type;
523 AdditiveAdder(BlockVector<T,A>& v, BlockVector<T,A>& x,
524 OverlappingAssigner<S>& assigner,
const field_type& relax_);
525 void operator()(
const size_type& domain);
527 static constexpr size_t n = std::decay_t<decltype(Impl::asVector(std::declval<T>()))>::dimension;
532 OverlappingAssigner<S>* assigner;
536 template<
typename S,
typename T>
537 struct MultiplicativeAdder
540 template<
typename S,
typename T,
typename A>
541 struct MultiplicativeAdder<S, BlockVector<T,A> >
543 typedef typename A::size_type size_type;
544 typedef typename std::decay_t<decltype(Impl::asVector(std::declval<T>()))>::field_type field_type;
545 MultiplicativeAdder(BlockVector<T,A>& v, BlockVector<T,A>& x,
546 OverlappingAssigner<S>& assigner_,
const field_type& relax_);
547 void operator()(
const size_type& domain);
549 static constexpr size_t n = std::decay_t<decltype(Impl::asVector(std::declval<T>()))>::dimension;
553 OverlappingAssigner<S>* assigner;
566 template<
typename T,
class X,
class S>
570 template<
class X,
class S>
573 typedef AdditiveAdder<S,X> Adder;
576 template<
class X,
class S>
577 struct AdderSelector<MultiplicativeSchwarzMode,X,S>
579 typedef MultiplicativeAdder<S,X> Adder;
582 template<
class X,
class S>
583 struct AdderSelector<SymmetricMultiplicativeSchwarzMode,X,S>
585 typedef MultiplicativeAdder<S,X> Adder;
599 template<
typename T1,
typename T2,
bool forward>
602 typedef T1 solver_vector;
603 typedef typename solver_vector::iterator solver_iterator;
604 typedef T2 subdomain_vector;
605 typedef typename subdomain_vector::const_iterator domain_iterator;
607 static solver_iterator begin(solver_vector& sv)
612 static solver_iterator end(solver_vector& sv)
616 static domain_iterator begin(
const subdomain_vector& sv)
621 static domain_iterator end(
const subdomain_vector& sv)
627 template<
typename T1,
typename T2>
630 typedef T1 solver_vector;
631 typedef typename solver_vector::reverse_iterator solver_iterator;
632 typedef T2 subdomain_vector;
633 typedef typename subdomain_vector::const_reverse_iterator domain_iterator;
635 static solver_iterator begin(solver_vector& sv)
640 static solver_iterator end(solver_vector& sv)
644 static domain_iterator begin(
const subdomain_vector& sv)
649 static domain_iterator end(
const subdomain_vector& sv)
667 typedef typename smoother::range_type range_type;
669 static void apply(smoother& sm, range_type& v,
const range_type& b)
671 sm.template apply<true>(v, b);
675 template<
class M,
class X,
class TD,
class TA>
679 typedef typename smoother::range_type range_type;
681 static void apply(smoother& sm, range_type& v,
const range_type& b)
683 sm.template apply<true>(v, b);
684 sm.template apply<false>(v, b);
688 template<
class T,
bool tag>
689 struct SeqOverlappingSchwarzAssemblerHelper
693 using SeqOverlappingSchwarzAssembler = SeqOverlappingSchwarzAssemblerHelper<T,Dune::StoresColumnCompressed<T>::value>;
695 template<
class K,
class Al,
class X,
class Y>
696 struct SeqOverlappingSchwarzAssemblerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,false>
698 typedef BCRSMatrix< K, Al> matrix_type;
699 static constexpr size_t n = std::decay_t<decltype(Impl::asMatrix(std::declval<K>()))>::rows;
700 template<
class RowToDomain,
class Solvers,
class SubDomains>
701 static std::size_t assembleLocalProblems(
const RowToDomain& rowToDomain,
const matrix_type& mat,
702 Solvers& solvers,
const SubDomains& domains,
706 template<
template<
class>
class S,
typename T,
typename A>
707 struct SeqOverlappingSchwarzAssemblerHelper<S<BCRSMatrix<T,A>>,true>
709 typedef BCRSMatrix<T,A> matrix_type;
710 static constexpr size_t n = std::decay_t<decltype(Impl::asMatrix(std::declval<T>()))>::rows;
711 template<
class RowToDomain,
class Solvers,
class SubDomains>
712 static std::size_t assembleLocalProblems(
const RowToDomain& rowToDomain,
const matrix_type& mat,
713 Solvers& solvers,
const SubDomains& domains,
717 template<
class M,
class X,
class Y>
718 struct SeqOverlappingSchwarzAssemblerILUBase
720 typedef M matrix_type;
721 template<
class RowToDomain,
class Solvers,
class SubDomains>
722 static std::size_t assembleLocalProblems(
const RowToDomain& rowToDomain,
const matrix_type& mat,
723 Solvers& solvers,
const SubDomains& domains,
727 template<
class M,
class X,
class Y>
728 struct SeqOverlappingSchwarzAssemblerHelper<ILU0SubdomainSolver<M,X,Y>,false>
729 :
public SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>
732 template<
class M,
class X,
class Y>
733 struct SeqOverlappingSchwarzAssemblerHelper<ILUNSubdomainSolver<M,X,Y>,false>
734 :
public SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>
747 template<
class M,
class X,
class TM=AdditiveSchwarzMode,
748 class TD=ILU0SubdomainSolver<M,X,X>,
class TA=std::allocator<X> >
788 typedef std::set<size_type, std::less<size_type>,
789 typename TA::template rebind<size_type>::other>
793 typedef std::vector<subdomain_type, typename TA::template rebind<subdomain_type>::other>
subdomain_vector;
799 typedef std::vector<subdomain_list, typename TA::template rebind<subdomain_list>::other >
rowtodomain_vector;
805 typedef std::vector<slu, typename TA::template rebind<slu>::other>
slu_vector;
821 field_type relaxationFactor=1,
bool onTheFly_=
true);
835 field_type relaxationFactor=1,
bool onTheFly_=
true);
842 virtual void pre (X& x, X& b)
853 virtual void apply (X& v,
const X& d);
865 template<
bool forward>
880 typename M::size_type maxlength;
887 template<
class I,
class S,
class D>
888 OverlappingSchwarzInitializer<I,S,D>::OverlappingSchwarzInitializer(InitializerList& il,
890 const subdomain_vector& domains_)
891 : initializers(&il), indices(&idx), indexMaps(il.size()), domains(domains_)
895 template<
class I,
class S,
class D>
896 void OverlappingSchwarzInitializer<I,S,D>::addRowNnz(
const Iter& row)
898 typedef typename IndexSet::value_type::const_iterator iterator;
899 for(iterator domain=(*indices)[row.index()].begin(); domain != (*indices)[row.index()].end(); ++domain) {
900 (*initializers)[*domain].addRowNnz(row, domains[*domain]);
901 indexMaps[*domain].insert(row.index());
905 template<
class I,
class S,
class D>
906 void OverlappingSchwarzInitializer<I,S,D>::allocate()
908 for(
auto&& i: *initializers)
909 i.allocateMatrixStorage();
910 for(
auto&& i: *initializers)
914 template<
class I,
class S,
class D>
915 void OverlappingSchwarzInitializer<I,S,D>::countEntries(
const Iter& row,
const CIter& col)
const
917 typedef typename IndexSet::value_type::const_iterator iterator;
918 for(iterator domain=(*indices)[row.index()].begin(); domain != (*indices)[row.index()].end(); ++domain) {
919 typename std::map<size_type,size_type>::const_iterator v = indexMaps[*domain].
find(col.index());
920 if(v!= indexMaps[*domain].end()) {
921 (*initializers)[*domain].countEntries(indexMaps[*domain].find(col.index())->second);
926 template<
class I,
class S,
class D>
927 void OverlappingSchwarzInitializer<I,S,D>::calcColstart()
const
929 for(
auto&& i : *initializers)
933 template<
class I,
class S,
class D>
934 void OverlappingSchwarzInitializer<I,S,D>::copyValue(
const Iter& row,
const CIter& col)
const
936 typedef typename IndexSet::value_type::const_iterator iterator;
937 for(iterator domain=(*indices)[row.index()].begin(); domain!= (*indices)[row.index()].end(); ++domain) {
938 typename std::map<size_type,size_type>::const_iterator v = indexMaps[*domain].find(col.index());
939 if(v!= indexMaps[*domain].end()) {
940 assert(indexMaps[*domain].end()!=indexMaps[*domain].find(row.index()));
941 (*initializers)[*domain].copyValue(col, indexMaps[*domain].find(row.index())->second,
947 template<
class I,
class S,
class D>
948 void OverlappingSchwarzInitializer<I,S,D>::createMatrix()
const
950 std::vector<IndexMap>().swap(indexMaps);
951 for(
auto&& i: *initializers)
955 template<
class I,
class S,
class D>
956 OverlappingSchwarzInitializer<I,S,D>::IndexMap::IndexMap()
960 template<
class I,
class S,
class D>
961 void OverlappingSchwarzInitializer<I,S,D>::IndexMap::insert(size_type grow)
963 assert(map_.find(grow)==map_.end());
964 map_.insert(std::make_pair(grow, row++));
967 template<
class I,
class S,
class D>
968 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator
969 OverlappingSchwarzInitializer<I,S,D>::IndexMap::find(size_type grow)
const
971 return map_.find(grow);
974 template<
class I,
class S,
class D>
975 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator
976 OverlappingSchwarzInitializer<I,S,D>::IndexMap::find(size_type grow)
978 return map_.find(grow);
981 template<
class I,
class S,
class D>
982 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator
983 OverlappingSchwarzInitializer<I,S,D>::IndexMap::end()
const
988 template<
class I,
class S,
class D>
989 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator
990 OverlappingSchwarzInitializer<I,S,D>::IndexMap::end()
995 template<
class I,
class S,
class D>
996 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator
997 OverlappingSchwarzInitializer<I,S,D>::IndexMap::begin()
const
1002 template<
class I,
class S,
class D>
1003 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator
1004 OverlappingSchwarzInitializer<I,S,D>::IndexMap::begin()
1006 return map_.begin();
1009 template<
class M,
class X,
class TM,
class TD,
class TA>
1012 : mat(mat_), relax(relaxationFactor), onTheFly(fly)
1014 typedef typename rowtodomain_vector::const_iterator RowDomainIterator;
1016#ifdef DUNE_ISTL_WITH_CHECKING
1017 assert(rowToDomain.size()==mat.N());
1018 assert(rowToDomain.size()==mat.M());
1020 for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter)
1021 assert(iter->size()>0);
1026 for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter)
1027 for(DomainIterator d=iter->
begin(); d != iter->
end(); ++d)
1031 solvers.resize(domains);
1032 subDomains.resize(domains);
1036 for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter, ++row)
1037 for(DomainIterator d=iter->
begin(); d != iter->
end(); ++d)
1038 subDomains[*d].insert(row);
1040#ifdef DUNE_ISTL_WITH_CHECKING
1042 typedef typename subdomain_vector::const_iterator iterator;
1043 for(iterator iter=subDomains.begin(); iter != subDomains.end(); ++iter) {
1044 typedef typename subdomain_type::const_iterator entry_iterator;
1046 for(entry_iterator entry = iter->begin(); entry != iter->end(); ++entry) {
1052 maxlength = SeqOverlappingSchwarzAssembler<slu>
1053 ::assembleLocalProblems(rowToDomain, mat, solvers, subDomains, onTheFly);
1056 template<
class M,
class X,
class TM,
class TD,
class TA>
1061 : mat(mat_), solvers(sd.size()), subDomains(sd), relax(relaxationFactor),
1064 typedef typename subdomain_vector::const_iterator DomainIterator;
1066#ifdef DUNE_ISTL_WITH_CHECKING
1069 for(DomainIterator d=sd.
begin(); d != sd.
end(); ++d,++i) {
1071 assert(d->size()>0);
1072 typedef typename DomainIterator::value_type::const_iterator entry_iterator;
1074 for(entry_iterator entry = d->
begin(); entry != d->
end(); ++entry) {
1087 for(DomainIterator domain=sd.begin(); domain != sd.end(); ++domain, ++domainId) {
1088 typedef typename subdomain_type::const_iterator iterator;
1089 for(iterator row=domain->begin(); row != domain->end(); ++row)
1090 rowToDomain[*row].push_back(domainId);
1093 maxlength = SeqOverlappingSchwarzAssembler<slu>
1094 ::assembleLocalProblems(rowToDomain, mat, solvers, subDomains, onTheFly);
1106 template<
typename T,
typename A>
1109 static constexpr size_t n = std::decay_t<decltype(Impl::asMatrix(std::declval<T>()))>::rows;
1110 static constexpr size_t m = std::decay_t<decltype(Impl::asMatrix(std::declval<T>()))>::cols;
1111 template<
class Domain>
1112 static int size(
const Domain & d)
1119 template<
class K,
class Al,
class X,
class Y>
1120 template<
class RowToDomain,
class Solvers,
class SubDomains>
1122 SeqOverlappingSchwarzAssemblerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>::
1123 assembleLocalProblems(
const RowToDomain& rowToDomain,
1124 const matrix_type& mat,
1126 const SubDomains& subDomains,
1133 typedef typename SubDomains::const_iterator DomainIterator;
1134 std::size_t maxlength = 0;
1138 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end(); ++domain)
1139 maxlength=
std::max(maxlength, domain->size());
1145#if HAVE_SUPERLU || HAVE_SUITESPARSE_UMFPACK
1146 template<
template<
class>
class S,
typename T,
typename A>
1147 template<
class RowToDomain,
class Solvers,
class SubDomains>
1148 std::size_t SeqOverlappingSchwarzAssemblerHelper<S<BCRSMatrix<T,A>>,
true>::assembleLocalProblems(
const RowToDomain& rowToDomain,
1149 const matrix_type& mat,
1151 const SubDomains& subDomains,
1154 typedef typename S<BCRSMatrix<T,A>>::MatrixInitializer MatrixInitializer;
1155 typedef typename std::vector<MatrixInitializer>::iterator InitializerIterator;
1156 typedef typename SubDomains::const_iterator DomainIterator;
1157 typedef typename Solvers::iterator SolverIterator;
1158 std::size_t maxlength = 0;
1161 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end(); ++domain)
1162 maxlength=
std::max(maxlength, domain->size());
1163 maxlength*=Impl::asMatrix(*mat[0].begin()).N();
1166 DomainIterator domain=subDomains.begin();
1169 std::vector<MatrixInitializer> initializers(subDomains.size());
1171 SolverIterator solver=solvers.begin();
1172 for(InitializerIterator initializer=initializers.begin(); initializer!=initializers.end();
1173 ++initializer, ++solver, ++domain) {
1174 solver->getInternalMatrix().N_=SeqOverlappingSchwarzDomainSize<matrix_type>::size(*domain);
1175 solver->getInternalMatrix().M_=SeqOverlappingSchwarzDomainSize<matrix_type>::size(*domain);
1177 *initializer=MatrixInitializer(solver->getInternalMatrix());
1181 typedef OverlappingSchwarzInitializer<std::vector<MatrixInitializer>,
1182 RowToDomain, SubDomains> Initializer;
1184 Initializer initializer(initializers, rowToDomain, subDomains);
1185 copyToColCompMatrix(initializer, mat);
1188 for(
auto&& s: solvers)
1190 for (SolverIterator solverIt = solvers.begin(); solverIt != solvers.end(); ++solverIt)
1192 assert(solverIt->getInternalMatrix().N() == solverIt->getInternalMatrix().M());
1193 maxlength =
std::max(maxlength, solverIt->getInternalMatrix().N());
1201 template<
class M,
class X,
class Y>
1202 template<
class RowToDomain,
class Solvers,
class SubDomains>
1203 std::size_t SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>::assembleLocalProblems(
const RowToDomain& rowToDomain,
1204 const matrix_type& mat,
1206 const SubDomains& subDomains,
1210 typedef typename SubDomains::const_iterator DomainIterator;
1211 typedef typename Solvers::iterator SolverIterator;
1212 std::size_t maxlength = 0;
1215 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end(); ++domain)
1216 maxlength=
std::max(maxlength, domain->size());
1219 SolverIterator solver=solvers.begin();
1220 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end();
1221 ++domain, ++solver) {
1222 solver->setSubMatrix(mat, *domain);
1223 maxlength=
std::max(maxlength, domain->size());
1232 template<
class M,
class X,
class TM,
class TD,
class TA>
1238 template<
class M,
class X,
class TM,
class TD,
class TA>
1239 template<
bool forward>
1242 typedef slu_vector solver_vector;
1243 typedef typename IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::solver_iterator iterator;
1244 typedef typename IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::domain_iterator
1247 OverlappingAssigner<TD> assigner(maxlength, mat, b, x);
1255 Adder adder(v, x, assigner, relax);
1259 std::for_each(domain->begin(), domain->end(), assigner);
1260 assigner.resetIndexForNextDomain();
1264 sdsolver.setSubMatrix(mat, *domain);
1266 sdsolver.apply(assigner.lhs(), assigner.rhs());
1268 solver->apply(assigner.lhs(), assigner.rhs());
1273 std::for_each(domain->begin(), domain->end(), adder);
1274 assigner.resetIndexForNextDomain();
1279 assigner.deallocate();
1282 template<
class K,
class Al,
class X,
class Y>
1283 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1285 const X& b_, Y& x_) :
1292 maxlength_(maxlength)
1295 template<
class K,
class Al,
class X,
class Y>
1297 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1304 template<
class K,
class Al,
class X,
class Y>
1306 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1307 ::resetIndexForNextDomain()
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>
1328 template<
class K,
class Al,
class X,
class Y>
1330 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1331 ::relaxResult(field_type relax)
1336 template<
class K,
class Al,
class X,
class Y>
1338 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1339 ::operator()(
const size_type& domainIndex)
1344 for(size_type j=0; j<n; ++j, ++i) {
1345 assert(i<maxlength_);
1346 rhs()[i]=(*b)[domainIndex][j];
1353 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1354 block_type tmp(0.0);
1355 (*col).mv((*x)[col.index()], tmp);
1357 for(size_type j=0; j<n; ++j, ++i) {
1358 assert(i<maxlength_);
1364 for(size_type j=0; j<n; ++j, ++i) {
1365 assert(i<maxlength_);
1366 rhs()[i]=Impl::asVector((*b)[domainIndex])[j];
1372 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1373 for(size_type k=0; k<n; ++k) {
1374 rhs()[i]-=Impl::asMatrix(*col)[j][k] * Impl::asVector((*x)[col.index()])[k];
1381 template<
class K,
class Al,
class X,
class Y>
1383 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< K, Al>, X, Y >,
false>
1384 ::assignResult(block_type& res)
1387 for(size_type j=0; j<n; ++j, ++i) {
1388 assert(i<maxlength_);
1389 Impl::asVector(res)[j]+=lhs()[i];
1393#if HAVE_SUPERLU || HAVE_SUITESPARSE_UMFPACK
1395 template<
template<
class>
class S,
typename T,
typename A>
1396 OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>
1397 ::OverlappingAssignerHelper(std::size_t maxlength,
1399 const range_type& b_,
1403 x(&x_), i(0), maxlength_(maxlength)
1405 rhs_ =
new field_type[maxlength];
1406 lhs_ =
new field_type[maxlength];
1410 template<
template<
class>
class S,
typename T,
typename A>
1411 void OverlappingAssignerHelper<S<BCRSMatrix<T,A> >,
true>::deallocate()
1417 template<
template<
class>
class S,
typename T,
typename A>
1418 void OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::operator()(
const size_type& domainIndex)
1423 for(size_type j=0; j<n; ++j, ++i) {
1424 assert(i<maxlength_);
1425 rhs_[i]=Impl::asVector((*b)[domainIndex])[j];
1433 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1435 Impl::asMatrix(*col).mv((*x)[col.index()], tmp);
1437 for(size_type j=0; j<n; ++j, ++i) {
1438 assert(i<maxlength_);
1439 rhs_[i]-=Impl::asVector(tmp)[j];
1446 template<
template<
class>
class S,
typename T,
typename A>
1447 void OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::relaxResult(field_type relax)
1449 for(size_type j=i+n; i<j; ++i) {
1450 assert(i<maxlength_);
1456 template<
template<
class>
class S,
typename T,
typename A>
1457 void OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::assignResult(block_type& res)
1460 for(size_type j=0; j<n; ++j, ++i) {
1461 assert(i<maxlength_);
1462 Impl::asVector(res)[j]+=lhs_[i];
1466 template<
template<
class>
class S,
typename T,
typename A>
1467 void OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::resetIndexForNextDomain()
1472 template<
template<
class>
class S,
typename T,
typename A>
1473 typename OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::field_type*
1474 OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::lhs()
1479 template<
template<
class>
class S,
typename T,
typename A>
1480 typename OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::field_type*
1481 OverlappingAssignerHelper<S<BCRSMatrix<T,A>>,
true>::rhs()
1488 template<
class M,
class X,
class Y>
1489 OverlappingAssignerILUBase<M,X,Y>::OverlappingAssignerILUBase(std::size_t maxlength,
1497 rhs_=
new Y(maxlength);
1498 lhs_ =
new X(maxlength);
1501 template<
class M,
class X,
class Y>
1502 void OverlappingAssignerILUBase<M,X,Y>::deallocate()
1508 template<
class M,
class X,
class Y>
1509 void OverlappingAssignerILUBase<M,X,Y>::operator()(
const size_type& domainIndex)
1511 (*rhs_)[i]=(*b)[domainIndex];
1514 typedef typename matrix_type::ConstColIterator col_iterator;
1517 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1518 Impl::asMatrix(*col).mmv((*x)[col.index()], (*rhs_)[i]);
1524 template<
class M,
class X,
class Y>
1525 void OverlappingAssignerILUBase<M,X,Y>::relaxResult(field_type relax)
1530 template<
class M,
class X,
class Y>
1531 void OverlappingAssignerILUBase<M,X,Y>::assignResult(block_type& res)
1536 template<
class M,
class X,
class Y>
1537 X& OverlappingAssignerILUBase<M,X,Y>::lhs()
1542 template<
class M,
class X,
class Y>
1543 Y& OverlappingAssignerILUBase<M,X,Y>::rhs()
1548 template<
class M,
class X,
class Y>
1549 void OverlappingAssignerILUBase<M,X,Y>::resetIndexForNextDomain()
1554 template<
typename S,
typename T,
typename A>
1557 OverlappingAssigner<S>& assigner_,
1558 const field_type& relax_)
1559 : v(&v_), x(&x_), assigner(&assigner_), relax(relax_)
1562 template<
typename S,
typename T,
typename A>
1563 void AdditiveAdder<S,BlockVector<T,A> >::operator()(
const size_type& domainIndex)
1566 assigner->assignResult((*v)[domainIndex]);
1570 template<
typename S,
typename T,
typename A>
1571 void AdditiveAdder<S,BlockVector<T,A> >::axpy()
1578 template<
typename S,
typename T,
typename A>
1579 MultiplicativeAdder<S,BlockVector<T,A> >
1580 ::MultiplicativeAdder(BlockVector<T,A>& v_,
1581 BlockVector<T,A>& x_,
1582 OverlappingAssigner<S>& assigner_,
const field_type& relax_)
1583 : x(&x_), assigner(&assigner_), relax(relax_)
1589 template<
typename S,
typename T,
typename A>
1590 void MultiplicativeAdder<S,BlockVector<T,A> >::operator()(
const size_type& domainIndex)
1593 assigner->relaxResult(relax);
1594 assigner->assignResult((*x)[domainIndex]);
1598 template<
typename S,
typename T,
typename A>
1599 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:425
row_type::ConstIterator ConstColIterator
Const iterator to the entries of a row.
Definition: bcrsmatrix.hh:699
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:136
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:43
D subdomain_vector
The vector type containing the subdomain to row index mapping.
Definition: overlappingschwarz.hh:46
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:751
X::field_type field_type
The field type of the preconditioner.
Definition: overlappingschwarz.hh:779
void apply(X &v, const X &d)
Apply one step of the preconditioner to the system A(v)=d.
M matrix_type
The type of the matrix to precondition.
Definition: overlappingschwarz.hh:756
TM Mode
The mode (additive or multiplicative) of the Schwarz method.
Definition: overlappingschwarz.hh:774
X range_type
The range type of the preconditioner.
Definition: overlappingschwarz.hh:766
X domain_type
The domain type of the preconditioner.
Definition: overlappingschwarz.hh:761
TD slu
The type for the subdomain solver in use.
Definition: overlappingschwarz.hh:802
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: overlappingschwarz.hh:869
virtual void pre(X &x, X &b)
Prepare the preconditioner.
Definition: overlappingschwarz.hh:842
virtual void post(X &x)
Postprocess the preconditioner.
Definition: overlappingschwarz.hh:860
std::vector< subdomain_type, typename TA::template rebind< subdomain_type >::other > subdomain_vector
The vector type containing the subdomain to row index mapping.
Definition: overlappingschwarz.hh:793
TA allocator
The allocator to use.
Definition: overlappingschwarz.hh:785
SLList< size_type, typename TA::template rebind< size_type >::other > subdomain_list
The type for the row to subdomain mapping.
Definition: overlappingschwarz.hh:796
std::vector< subdomain_list, typename TA::template rebind< subdomain_list >::other > rowtodomain_vector
The vector type containing the row index to subdomain mapping.
Definition: overlappingschwarz.hh:799
std::vector< slu, typename TA::template rebind< slu >::other > slu_vector
The vector type containing subdomain solvers.
Definition: overlappingschwarz.hh:805
matrix_type::size_type size_type
The return type of the size method.
Definition: overlappingschwarz.hh:782
std::set< size_type, std::less< size_type >, typename TA::template rebind< size_type >::other > subdomain_type
The type for the subdomain to row index mapping.
Definition: overlappingschwarz.hh:790
This file implements a dense matrix with dynamic numbers of rows and columns.
decltype(auto) apply(F &&f, ArgTuple &&args)
Apply function with arguments given as tuple.
Definition: apply.hh:46
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:25
virtual void apply(X &v, const X &d)
Apply the precondtioner.
Definition: overlappingschwarz.hh:1233
SeqOverlappingSchwarz(const matrix_type &mat, const subdomain_vector &subDomains, field_type relaxationFactor=1, bool onTheFly_=true)
Construct the overlapping Schwarz method.
Definition: overlappingschwarz.hh:1057
SeqOverlappingSchwarz(const matrix_type &mat, const rowtodomain_vector &rowToDomain, field_type relaxationFactor=1, bool onTheFly_=true)
Definition: overlappingschwarz.hh:1010
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:14
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:568
Tag that the tells the Schwarz method to be additive.
Definition: overlappingschwarz.hh:116
Helper template meta program for application of overlapping Schwarz.
Definition: overlappingschwarz.hh:601
Tag that tells the Schwarz method to be multiplicative.
Definition: overlappingschwarz.hh:122
Helper template meta program for application of overlapping Schwarz.
Definition: overlappingschwarz.hh:665
Definition: overlappingschwarz.hh:1104
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:129
Classes for using SuperLU with ISTL matrices.
Classes for using UMFPack with ISTL matrices.
Definition of the DUNE_UNUSED macro for the case that config.h is not available.