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;
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,
int n,
class Al,
class X,
class Y>
145 typedef typename Dune::remove_const<M>::type matrix_type;
146 typedef K field_type;
147 typedef typename Dune::remove_const<M>::type rilu_type;
149 typedef X domain_type;
151 typedef Y range_type;
159 assert(v.
size() > 0);
161 assert(A.rows() <= v.
size());
162 assert(A.cols() <= v.
size());
163 size_t sz = A.rows();
179 void setSubMatrix(
const M& BCRS, S& rowset)
181 size_t sz = rowset.size();
183 typedef typename S::const_iterator SIter;
185 for(SIter rowIdx = rowset.begin(), rowEnd=rowset.end();
186 rowIdx!= rowEnd; ++rowIdx, r++)
189 for(SIter colIdx = rowset.begin(), colEnd=rowset.end();
190 colIdx!= colEnd; ++colIdx, c++)
192 if (BCRS[*rowIdx].find(*colIdx) == BCRS[*rowIdx].end())
194 for (
size_t i=0; i<n; i++)
196 for (
size_t j=0; j<n; j++)
198 A[r*n+i][c*n+j] = BCRS[*rowIdx][*colIdx][i][j];
208 template<
typename T,
bool tag>
209 class OverlappingAssignerHelper
213 struct OverlappingAssigner :
public OverlappingAssignerHelper<T,Dune::StoresColumnCompressed<T>::value >
216 OverlappingAssigner(std::size_t maxlength,
const typename T::matrix_type& mat,
217 const typename T::range_type& b,
typename T::range_type& x)
218 : OverlappingAssignerHelper<T,
Dune::StoresColumnCompressed<T>::value >
224 template<
class K,
int n,
class Al,
class X,
class Y>
225 class OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
228 typedef BCRSMatrix< FieldMatrix<K,n,n>, Al> matrix_type;
229 typedef K field_type;
230 typedef Y range_type;
231 typedef typename range_type::block_type block_type;
232 typedef typename matrix_type::size_type size_type;
241 OverlappingAssignerHelper(std::size_t maxlength,
const BCRSMatrix<FieldMatrix<K,n,n>, Al>& mat_,
const X& b_, Y& x_);
253 void resetIndexForNextDomain();
260 DynamicVector<K> & lhs();
267 DynamicVector<K> & rhs();
274 void relaxResult(field_type relax);
280 void operator()(
const size_type& domainIndex);
290 void assignResult(block_type& res);
296 const matrix_type* mat;
299 DynamicVector<field_type> * rhs_;
302 DynamicVector<field_type> * lhs_;
310 std::size_t maxlength_;
313#if HAVE_SUPERLU || HAVE_UMFPACK
314 template<
template<
class>
class S,
int n,
int m,
typename T,
typename A>
315 struct OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>, A> >, true>
317 typedef BCRSMatrix<FieldMatrix<T,n,m>, A> matrix_type;
318 typedef typename S<BCRSMatrix<FieldMatrix<T,n,m>, A> >::range_type range_type;
319 typedef typename range_type::field_type field_type;
320 typedef typename range_type::block_type block_type;
322 typedef typename matrix_type::size_type size_type;
331 OverlappingAssignerHelper(std::size_t maxlength,
const matrix_type& mat,
332 const range_type& b, range_type& x);
343 void resetIndexForNextDomain();
361 void relaxResult(field_type relax);
367 void operator()(
const size_type& domain);
376 void assignResult(block_type& res);
382 const matrix_type* mat;
394 std::size_t maxlength_;
399 template<
class M,
class X,
class Y>
400 class OverlappingAssignerILUBase
403 typedef M matrix_type;
405 typedef typename M::field_type field_type;
407 typedef typename Y::block_type block_type;
409 typedef typename matrix_type::size_type size_type;
417 OverlappingAssignerILUBase(std::size_t maxlength,
const M& mat,
429 void resetIndexForNextDomain();
447 void relaxResult(field_type relax);
453 void operator()(
const size_type& domain);
462 void assignResult(block_type& res);
482 template<
class M,
class X,
class Y>
483 class OverlappingAssignerHelper<ILU0SubdomainSolver<M,X,Y>, false>
484 :
public OverlappingAssignerILUBase<M,X,Y>
494 OverlappingAssignerHelper(std::size_t maxlength,
const M& mat,
496 : OverlappingAssignerILUBase<M,X,Y>(maxlength, mat,b,x)
501 template<
class M,
class X,
class Y>
502 class OverlappingAssignerHelper<ILUNSubdomainSolver<M,X,Y>,false>
503 :
public OverlappingAssignerILUBase<M,X,Y>
513 OverlappingAssignerHelper(std::size_t maxlength,
const M& mat,
515 : OverlappingAssignerILUBase<M,X,Y>(maxlength, mat,b,x)
519 template<
typename S,
typename T>
523 template<
typename S,
typename T,
typename A,
int n>
524 struct AdditiveAdder<S, BlockVector<FieldVector<T,n>,A> >
526 typedef typename A::size_type size_type;
527 AdditiveAdder(BlockVector<FieldVector<T,n>,A>& v, BlockVector<FieldVector<T,n>,A>& x,
528 OverlappingAssigner<S>& assigner,
const T& relax_);
529 void operator()(
const size_type& domain);
533 BlockVector<FieldVector<T,n>,A>* v;
534 BlockVector<FieldVector<T,n>,A>* x;
535 OverlappingAssigner<S>* assigner;
539 template<
typename S,
typename T>
540 struct MultiplicativeAdder
543 template<
typename S,
typename T,
typename A,
int n>
544 struct MultiplicativeAdder<S, BlockVector<FieldVector<T,n>,A> >
546 typedef typename A::size_type size_type;
547 MultiplicativeAdder(BlockVector<FieldVector<T,n>,A>& v, BlockVector<FieldVector<T,n>,A>& x,
548 OverlappingAssigner<S>& assigner_,
const T& relax_);
549 void operator()(
const size_type& domain);
553 BlockVector<FieldVector<T,n>,A>* x;
554 OverlappingAssigner<S>* assigner;
567 template<
typename T,
class X,
class S>
571 template<
class X,
class S>
574 typedef AdditiveAdder<S,X> Adder;
577 template<
class X,
class S>
578 struct AdderSelector<MultiplicativeSchwarzMode,X,S>
580 typedef MultiplicativeAdder<S,X> Adder;
583 template<
class X,
class S>
584 struct AdderSelector<SymmetricMultiplicativeSchwarzMode,X,S>
586 typedef MultiplicativeAdder<S,X> Adder;
600 template<
typename T1,
typename T2,
bool forward>
603 typedef T1 solver_vector;
604 typedef typename solver_vector::iterator solver_iterator;
605 typedef T2 subdomain_vector;
606 typedef typename subdomain_vector::const_iterator domain_iterator;
608 static solver_iterator begin(solver_vector& sv)
613 static solver_iterator end(solver_vector& sv)
617 static domain_iterator begin(
const subdomain_vector& sv)
622 static domain_iterator end(
const subdomain_vector& sv)
628 template<
typename T1,
typename T2>
631 typedef T1 solver_vector;
632 typedef typename solver_vector::reverse_iterator solver_iterator;
633 typedef T2 subdomain_vector;
634 typedef typename subdomain_vector::const_reverse_iterator domain_iterator;
636 static solver_iterator begin(solver_vector& sv)
641 static solver_iterator end(solver_vector& sv)
645 static domain_iterator begin(
const subdomain_vector& sv)
650 static domain_iterator end(
const subdomain_vector& sv)
668 typedef typename smoother::range_type range_type;
670 static void apply(smoother& sm, range_type& v,
const range_type& b)
672 sm.template apply<true>(v, b);
676 template<
class M,
class X,
class TD,
class TA>
680 typedef typename smoother::range_type range_type;
682 static void apply(smoother& sm, range_type& v,
const range_type& b)
684 sm.template apply<true>(v, b);
685 sm.template apply<false>(v, b);
689 template<
class T,
bool tag>
690 struct SeqOverlappingSchwarzAssemblerHelper
694 struct SeqOverlappingSchwarzAssembler
695 :
public SeqOverlappingSchwarzAssemblerHelper<T,Dune::StoresColumnCompressed<T>::value>
699 template<
class K,
int n,
class Al,
class X,
class Y>
700 struct SeqOverlappingSchwarzAssemblerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
702 typedef BCRSMatrix< FieldMatrix<K,n,n>, Al> matrix_type;
703 template<
class RowToDomain,
class Solvers,
class SubDomains>
704 static std::size_t assembleLocalProblems(
const RowToDomain& rowToDomain,
const matrix_type& mat,
705 Solvers& solvers,
const SubDomains& domains,
709 template<
template<
class>
class S,
typename T,
typename A,
int m,
int n>
710 struct SeqOverlappingSchwarzAssemblerHelper<S<BCRSMatrix<FieldMatrix<T,m,n>,A> >,true>
712 typedef BCRSMatrix<FieldMatrix<T,m,n>,A> matrix_type;
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 TA::template rebind<size_type>::other>
795 typedef std::vector<subdomain_type, typename TA::template rebind<subdomain_type>::other>
subdomain_vector;
801 typedef std::vector<subdomain_list, typename TA::template rebind<subdomain_list>::other >
rowtodomain_vector;
807 typedef std::vector<slu, typename TA::template rebind<slu>::other>
slu_vector;
828 field_type relaxationFactor=1,
bool onTheFly_=
true);
842 field_type relaxationFactor=1,
bool onTheFly_=
true);
849 virtual void pre (X& x, X& b)
860 virtual void apply (X& v,
const X& d);
872 template<
bool forward>
881 typename M::size_type maxlength;
888 template<
class I,
class S,
class D>
891 const subdomain_vector& domains_)
892 : initializers(&il), indices(&idx), indexMaps(il.size()), domains(domains_)
896 template<
class I,
class S,
class D>
897 void OverlappingSchwarzInitializer<I,S,D>::addRowNnz(
const Iter& row)
899 typedef typename IndexSet::value_type::const_iterator iterator;
900 for(iterator domain=(*indices)[row.index()].begin(); domain != (*indices)[row.index()].end(); ++domain) {
901 (*initializers)[*domain].addRowNnz(row, domains[*domain]);
902 indexMaps[*domain].insert(row.index());
906 template<
class I,
class S,
class D>
907 void OverlappingSchwarzInitializer<I,S,D>::allocate()
909 std::for_each(initializers->begin(), initializers->end(),
910 std::mem_fun_ref(&AtomInitializer::allocateMatrixStorage));
911 std::for_each(initializers->begin(), initializers->end(),
912 std::mem_fun_ref(&AtomInitializer::allocateMarker));
915 template<
class I,
class S,
class D>
916 void OverlappingSchwarzInitializer<I,S,D>::countEntries(
const Iter& row,
const CIter& col)
const
918 typedef typename IndexSet::value_type::const_iterator iterator;
919 for(iterator domain=(*indices)[row.index()].begin(); domain != (*indices)[row.index()].end(); ++domain) {
920 typename std::map<size_type,size_type>::const_iterator v = indexMaps[*domain].
find(col.index());
921 if(v!= indexMaps[*domain].end()) {
922 (*initializers)[*domain].countEntries(indexMaps[*domain].find(col.index())->second);
927 template<
class I,
class S,
class D>
928 void OverlappingSchwarzInitializer<I,S,D>::calcColstart()
const
930 std::for_each(initializers->begin(), initializers->end(),
931 std::mem_fun_ref(&AtomInitializer::calcColstart));
934 template<
class I,
class S,
class D>
935 void OverlappingSchwarzInitializer<I,S,D>::copyValue(
const Iter& row,
const CIter& col)
const
937 typedef typename IndexSet::value_type::const_iterator iterator;
938 for(iterator domain=(*indices)[row.index()].begin(); domain!= (*indices)[row.index()].end(); ++domain) {
939 typename std::map<size_type,size_type>::const_iterator v = indexMaps[*domain].find(col.index());
940 if(v!= indexMaps[*domain].end()) {
941 assert(indexMaps[*domain].end()!=indexMaps[*domain].find(row.index()));
942 (*initializers)[*domain].copyValue(col, indexMaps[*domain].find(row.index())->second,
948 template<
class I,
class S,
class D>
949 void OverlappingSchwarzInitializer<I,S,D>::createMatrix()
const
951 std::vector<IndexMap>().swap(indexMaps);
952 std::for_each(initializers->begin(), initializers->end(),
953 std::mem_fun_ref(&AtomInitializer::createMatrix));
956 template<
class I,
class S,
class D>
957 OverlappingSchwarzInitializer<I,S,D>::IndexMap::IndexMap()
961 template<
class I,
class S,
class D>
962 void OverlappingSchwarzInitializer<I,S,D>::IndexMap::insert(size_type grow)
964 assert(map_.find(grow)==map_.end());
965 map_.insert(std::make_pair(grow, row++));
968 template<
class I,
class S,
class D>
969 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator
970 OverlappingSchwarzInitializer<I,S,D>::IndexMap::find(size_type grow)
const
972 return map_.find(grow);
975 template<
class I,
class S,
class D>
976 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator
977 OverlappingSchwarzInitializer<I,S,D>::IndexMap::find(size_type grow)
979 return map_.find(grow);
982 template<
class I,
class S,
class D>
983 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator
984 OverlappingSchwarzInitializer<I,S,D>::IndexMap::end()
const
989 template<
class I,
class S,
class D>
990 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator
991 OverlappingSchwarzInitializer<I,S,D>::IndexMap::end()
996 template<
class I,
class S,
class D>
997 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator
998 OverlappingSchwarzInitializer<I,S,D>::IndexMap::begin()
const
1000 return map_.begin();
1003 template<
class I,
class S,
class D>
1004 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator
1005 OverlappingSchwarzInitializer<I,S,D>::IndexMap::begin()
1007 return map_.begin();
1010 template<
class M,
class X,
class TM,
class TD,
class TA>
1013 : mat(mat_), relax(relaxationFactor), onTheFly(fly)
1015 typedef typename rowtodomain_vector::const_iterator RowDomainIterator;
1017#ifdef DUNE_ISTL_WITH_CHECKING
1018 assert(rowToDomain.size()==mat.N());
1019 assert(rowToDomain.size()==mat.M());
1021 for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter)
1022 assert(iter->size()>0);
1027 for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter)
1028 for(DomainIterator d=iter->
begin(); d != iter->
end(); ++d)
1029 domains=std::max(domains, *d);
1032 solvers.resize(domains);
1033 subDomains.resize(domains);
1037 for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter, ++row)
1038 for(DomainIterator d=iter->
begin(); d != iter->
end(); ++d)
1039 subDomains[*d].insert(row);
1041#ifdef DUNE_ISTL_WITH_CHECKING
1043 typedef typename subdomain_vector::const_iterator iterator;
1044 for(iterator iter=subDomains.begin(); iter != subDomains.end(); ++iter) {
1045 typedef typename subdomain_type::const_iterator entry_iterator;
1047 for(entry_iterator entry = iter->begin(); entry != iter->end(); ++entry) {
1053 maxlength = SeqOverlappingSchwarzAssembler<slu>
1054 ::assembleLocalProblems(rowToDomain, mat, solvers, subDomains, onTheFly);
1057 template<
class M,
class X,
class TM,
class TD,
class TA>
1062 : mat(mat_), solvers(sd.size()), subDomains(sd), relax(relaxationFactor),
1065 typedef typename subdomain_vector::const_iterator DomainIterator;
1067#ifdef DUNE_ISTL_WITH_CHECKING
1070 for(DomainIterator d=sd.
begin(); d != sd.
end(); ++d,++i) {
1072 assert(d->
size()>0);
1073 typedef typename DomainIterator::value_type::const_iterator entry_iterator;
1075 for(entry_iterator entry = d->
begin(); entry != d->
end(); ++entry) {
1088 for(DomainIterator domain=sd.begin(); domain != sd.end(); ++domain, ++domainId) {
1089 typedef typename subdomain_type::const_iterator iterator;
1090 for(iterator row=domain->begin(); row != domain->end(); ++row)
1091 rowToDomain[*row].push_back(domainId);
1094 maxlength = SeqOverlappingSchwarzAssembler<slu>
1095 ::assembleLocalProblems(rowToDomain, mat, solvers, subDomains, onTheFly);
1107 template<
typename T,
typename A,
int n,
int m>
1110 template<
class Domain>
1111 static int size(
const Domain & d)
1118 template<
class K,
int n,
class Al,
class X,
class Y>
1119 template<
class RowToDomain,
class Solvers,
class SubDomains>
1121 SeqOverlappingSchwarzAssemblerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,
false>::
1122 assembleLocalProblems(
const RowToDomain& rowToDomain,
1123 const matrix_type& mat,
1125 const SubDomains& subDomains,
1132 typedef typename SubDomains::const_iterator DomainIterator;
1133 std::size_t maxlength = 0;
1137 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end(); ++domain)
1138 maxlength=std::max(maxlength, domain->size());
1144#if HAVE_SUPERLU || HAVE_UMFPACK
1145 template<
template<
class>
class S,
typename T,
typename A,
int m,
int n>
1146 template<
class RowToDomain,
class Solvers,
class SubDomains>
1147 std::size_t SeqOverlappingSchwarzAssemblerHelper<S<BCRSMatrix<FieldMatrix<T,m,n>,A> >,
true>::assembleLocalProblems(
const RowToDomain& rowToDomain,
1148 const matrix_type& mat,
1150 const SubDomains& subDomains,
1153 typedef typename S<BCRSMatrix<FieldMatrix<T,m,n>,A> >::MatrixInitializer MatrixInitializer;
1154 typedef typename std::vector<MatrixInitializer>::iterator InitializerIterator;
1155 typedef typename SubDomains::const_iterator DomainIterator;
1156 typedef typename Solvers::iterator SolverIterator;
1157 std::size_t maxlength = 0;
1160 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end(); ++domain)
1161 maxlength=std::max(maxlength, domain->size());
1162 maxlength*=mat[0].
begin()->
N();
1165 DomainIterator domain=subDomains.begin();
1168 std::vector<MatrixInitializer> initializers(subDomains.size());
1170 SolverIterator solver=solvers.begin();
1171 for(InitializerIterator initializer=initializers.begin(); initializer!=initializers.end();
1172 ++initializer, ++solver, ++domain) {
1173 solver->getInternalMatrix().N_=SeqOverlappingSchwarzDomainSize<matrix_type>::size(*domain);
1174 solver->getInternalMatrix().M_=SeqOverlappingSchwarzDomainSize<matrix_type>::size(*domain);
1176 *initializer=MatrixInitializer(solver->getInternalMatrix());
1180 typedef OverlappingSchwarzInitializer<std::vector<MatrixInitializer>,
1181 RowToDomain, SubDomains> Initializer;
1183 Initializer initializer(initializers, rowToDomain, subDomains);
1184 copyToColCompMatrix(initializer, mat);
1187 std::for_each(solvers.begin(), solvers.end(), std::mem_fun_ref(&S<BCRSMatrix<FieldMatrix<T,m,n>,A> >::decompose));
1188 for (SolverIterator solverIt = solvers.begin(); solverIt != solvers.end(); ++solverIt)
1190 assert(solverIt->getInternalMatrix().N() == solverIt->getInternalMatrix().M());
1191 maxlength = std::max(maxlength, solverIt->getInternalMatrix().N());
1199 template<
class M,
class X,
class Y>
1200 template<
class RowToDomain,
class Solvers,
class SubDomains>
1201 std::size_t SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>::assembleLocalProblems(
const RowToDomain& rowToDomain,
1202 const matrix_type& mat,
1204 const SubDomains& subDomains,
1208 typedef typename SubDomains::const_iterator DomainIterator;
1209 typedef typename Solvers::iterator SolverIterator;
1210 std::size_t maxlength = 0;
1213 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end(); ++domain)
1214 maxlength=std::max(maxlength, domain->size());
1217 SolverIterator solver=solvers.begin();
1218 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end();
1219 ++domain, ++solver) {
1220 solver->setSubMatrix(mat, *domain);
1221 maxlength=std::max(maxlength, domain->size());
1230 template<
class M,
class X,
class TM,
class TD,
class TA>
1236 template<
class M,
class X,
class TM,
class TD,
class TA>
1237 template<
bool forward>
1240 typedef slu_vector solver_vector;
1241 typedef typename IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::solver_iterator iterator;
1242 typedef typename IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::domain_iterator
1245 OverlappingAssigner<TD> assigner(maxlength, mat, b, x);
1253 Adder adder(v, x, assigner, relax);
1257 std::for_each(domain->begin(), domain->end(), assigner);
1258 assigner.resetIndexForNextDomain();
1262 sdsolver.setSubMatrix(mat, *domain);
1264 sdsolver.apply(assigner.lhs(), assigner.rhs());
1266 solver->apply(assigner.lhs(), assigner.rhs());
1271 std::for_each(domain->begin(), domain->end(), adder);
1272 assigner.resetIndexForNextDomain();
1277 assigner.deallocate();
1280 template<
class K,
int n,
class Al,
class X,
class Y>
1281 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,
false>
1283 const X& b_, Y& x_) :
1290 maxlength_(maxlength)
1293 template<
class K,
int n,
class Al,
class X,
class Y>
1295 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,
false>
1302 template<
class K,
int n,
class Al,
class X,
class Y>
1304 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,
false>
1305 ::resetIndexForNextDomain()
1310 template<
class K,
int n,
class Al,
class X,
class Y>
1312 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,
false>
1318 template<
class K,
int n,
class Al,
class X,
class Y>
1320 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,
false>
1326 template<
class K,
int n,
class Al,
class X,
class Y>
1328 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,
false>
1329 ::relaxResult(field_type relax)
1334 template<
class K,
int n,
class Al,
class X,
class Y>
1336 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,
false>
1337 ::operator()(
const size_type& domainIndex)
1342 for(size_type j=0; j<n; ++j, ++i) {
1343 assert(i<maxlength_);
1344 rhs()[i]=(*b)[domainIndex][j];
1351 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1352 block_type tmp(0.0);
1353 (*col).mv((*x)[col.index()], tmp);
1355 for(size_type j=0; j<n; ++j, ++i) {
1356 assert(i<maxlength_);
1362 for(size_type j=0; j<n; ++j, ++i) {
1363 assert(i<maxlength_);
1364 rhs()[i]=(*b)[domainIndex][j];
1370 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1371 for(size_type k=0; k<n; ++k) {
1372 rhs()[i]-=(*col)[j][k] * (*x)[col.index()][k];
1379 template<
class K,
int n,
class Al,
class X,
class Y>
1381 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,
false>
1382 ::assignResult(block_type& res)
1385 for(size_type j=0; j<n; ++j, ++i) {
1386 assert(i<maxlength_);
1391#if HAVE_SUPERLU || HAVE_UMFPACK
1393 template<
template<
class>
class S,
int n,
int m,
typename T,
typename A>
1394 OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,
true>
1395 ::OverlappingAssignerHelper(std::size_t maxlength,
1397 const range_type& b_,
1401 x(&x_), i(0), maxlength_(maxlength)
1403 rhs_ =
new field_type[maxlength];
1404 lhs_ =
new field_type[maxlength];
1408 template<
template<
class>
class S,
int n,
int m,
typename T,
typename A>
1409 void OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,
true>::deallocate()
1415 template<
template<
class>
class S,
int n,
int m,
typename T,
typename A>
1416 void OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,
true>::operator()(
const size_type& domainIndex)
1421 for(size_type j=0; j<n; ++j, ++i) {
1422 assert(i<maxlength_);
1423 rhs_[i]=(*b)[domainIndex][j];
1431 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1433 (*col).mv((*x)[col.index()], tmp);
1435 for(size_type j=0; j<n; ++j, ++i) {
1436 assert(i<maxlength_);
1444 template<
template<
class>
class S,
int n,
int m,
typename T,
typename A>
1445 void OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,
true>::relaxResult(field_type relax)
1447 for(size_type j=i+n; i<j; ++i) {
1448 assert(i<maxlength_);
1454 template<
template<
class>
class S,
int n,
int m,
typename T,
typename A>
1455 void OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,
true>::assignResult(block_type& res)
1458 for(size_type j=0; j<n; ++j, ++i) {
1459 assert(i<maxlength_);
1464 template<
template<
class>
class S,
int n,
int m,
typename T,
typename A>
1465 void OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,
true>::resetIndexForNextDomain()
1470 template<
template<
class>
class S,
int n,
int m,
typename T,
typename A>
1471 typename OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,
true>::field_type*
1472 OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,
true>::lhs()
1477 template<
template<
class>
class S,
int n,
int m,
typename T,
typename A>
1478 typename OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,
true>::field_type*
1479 OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,
true>::rhs()
1486 template<
class M,
class X,
class Y>
1487 OverlappingAssignerILUBase<M,X,Y>::OverlappingAssignerILUBase(std::size_t maxlength,
1495 rhs_=
new Y(maxlength);
1496 lhs_ =
new X(maxlength);
1499 template<
class M,
class X,
class Y>
1500 void OverlappingAssignerILUBase<M,X,Y>::deallocate()
1506 template<
class M,
class X,
class Y>
1507 void OverlappingAssignerILUBase<M,X,Y>::operator()(
const size_type& domainIndex)
1509 (*rhs_)[i]=(*b)[domainIndex];
1512 typedef typename matrix_type::ConstColIterator col_iterator;
1515 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1516 (*col).mmv((*x)[col.index()], (*rhs_)[i]);
1522 template<
class M,
class X,
class Y>
1523 void OverlappingAssignerILUBase<M,X,Y>::relaxResult(field_type relax)
1528 template<
class M,
class X,
class Y>
1529 void OverlappingAssignerILUBase<M,X,Y>::assignResult(block_type& res)
1534 template<
class M,
class X,
class Y>
1535 X& OverlappingAssignerILUBase<M,X,Y>::lhs()
1540 template<
class M,
class X,
class Y>
1541 Y& OverlappingAssignerILUBase<M,X,Y>::rhs()
1546 template<
class M,
class X,
class Y>
1547 void OverlappingAssignerILUBase<M,X,Y>::resetIndexForNextDomain()
1552 template<
typename S,
typename T,
typename A,
int n>
1555 OverlappingAssigner<S>& assigner_,
1557 : v(&v_), x(&x_), assigner(&assigner_), relax(relax_)
1560 template<
typename S,
typename T,
typename A,
int n>
1561 void AdditiveAdder<S,BlockVector<FieldVector<T,n>,A> >::operator()(
const size_type& domainIndex)
1564 assigner->assignResult((*v)[domainIndex]);
1568 template<
typename S,
typename T,
typename A,
int n>
1569 void AdditiveAdder<S,BlockVector<FieldVector<T,n>,A> >::axpy()
1576 template<
typename S,
typename T,
typename A,
int n>
1577 MultiplicativeAdder<S,BlockVector<FieldVector<T,n>,A> >
1578 ::MultiplicativeAdder(BlockVector<FieldVector<T,n>,A>& v_,
1579 BlockVector<FieldVector<T,n>,A>& x_,
1580 OverlappingAssigner<S>& assigner_,
const T& relax_)
1581 : x(&x_), assigner(&assigner_), relax(relax_)
1587 template<
typename S,
typename T,
typename A,
int n>
1588 void MultiplicativeAdder<S,BlockVector<FieldVector<T,n>,A> >::operator()(
const size_type& domainIndex)
1591 assigner->relaxResult(relax);
1592 assigner->assignResult((*x)[domainIndex]);
1596 template<
typename S,
typename T,
typename A,
int n>
1597 void MultiplicativeAdder<S,BlockVector<FieldVector<T,n>,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:413
row_type::ConstIterator ConstColIterator
Const iterator to the entries of a row.
Definition: bcrsmatrix.hh:690
A vector of blocks with memory management.
Definition: bvector.hh:253
Iterator begin()
begin iterator
Definition: densevector.hh:307
size_type size() const
size method
Definition: densevector.hh:296
Iterator end()
end iterator
Definition: densevector.hh:313
Iterator find(size_type i)
return iterator to given element or end()
Definition: densevector.hh:333
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:126
A dense n x m matrix.
Definition: fmatrix.hh:67
Index Set Interface base class.
Definition: indexidset.hh:76
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:79
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:26
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.
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
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 void pre(X &x, X &b)
Prepare the preconditioner.
Definition: overlappingschwarz.hh:849
virtual void post(X &x)
Postprocess the preconditioner.
Definition: overlappingschwarz.hh:867
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:795
TA allocator
The allocator to use.
Definition: overlappingschwarz.hh:787
SLList< size_type, typename TA::template rebind< size_type >::other > subdomain_list
The type for the row to subdomain mapping.
Definition: overlappingschwarz.hh:798
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:801
std::vector< slu, typename TA::template rebind< slu >::other > slu_vector
The vector type containing subdomain solvers.
Definition: overlappingschwarz.hh:807
matrix_type::size_type size_type
The return type of the size method.
Definition: overlappingschwarz.hh:784
@ category
The category the precondtioner is part of.
Definition: overlappingschwarz.hh:811
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:792
RealIterator< const B > const_iterator
iterator class for sequential access
Definition: basearray.hh:195
size_type N() const
number of blocks in the vector (are of size 1 here)
Definition: bvector.hh:217
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:1231
SeqOverlappingSchwarz(const matrix_type &mat, const subdomain_vector &subDomains, field_type relaxationFactor=1, bool onTheFly_=true)
Construct the overlapping Schwarz method.
Definition: overlappingschwarz.hh:1058
SeqOverlappingSchwarz(const matrix_type &mat, const rowtodomain_vector &rowToDomain, field_type relaxationFactor=1, bool onTheFly_=true)
Definition: overlappingschwarz.hh:1011
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: alignment.hh:10
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:569
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:602
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:666
Definition: overlappingschwarz.hh:1105
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:21
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.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentional unused function parameters with.
Definition: unused.hh:18