3 #ifndef DUNE_OVERLAPPINGSCHWARZ_HH
4 #define DUNE_OVERLAPPINGSCHWARZ_HH
10 #include<dune/common/dynmatrix.hh>
11 #include<dune/common/sllist.hh>
32 template<
class M,
class X,
class TM,
class TD,
class TA>
38 template<
class I,
class S,
class D>
47 typedef typename AtomInitializer::Matrix
Matrix;
48 typedef typename Matrix::const_iterator
Iter;
49 typedef typename Matrix::row_type::const_iterator
CIter;
75 typedef std::map<size_type,size_type> Map;
76 typedef typename Map::iterator iterator;
77 typedef typename Map::const_iterator const_iterator;
83 const_iterator find(
size_type grow)
const;
89 const_iterator begin()
const;
93 const_iterator end()
const;
96 std::map<size_type,size_type> map_;
101 typedef typename InitializerList::iterator InitIterator;
102 typedef typename IndexSet::const_iterator IndexIteratur;
105 mutable std::vector<IndexMap> indexMaps;
132 template<
class M,
class X,
class Y>
136 template<
class K,
int n,
class Al,
class X,
class Y>
154 void apply (DynamicVector<field_type>& v, DynamicVector<field_type>& d)
156 assert(v.size() > 0);
157 assert(v.size() == d.size());
158 assert(
A.rows() <= v.size());
159 assert(
A.cols() <= v.size());
160 size_t sz =
A.rows();
164 v.resize(v.capacity());
165 d.resize(d.capacity());
176 void setSubMatrix(
const M& BCRS, S& rowset)
178 size_t sz = rowset.size();
180 typename DynamicMatrix<K>::RowIterator rIt =
A.begin();
181 typedef typename S::const_iterator SIter;
183 for(SIter rowIdx = rowset.begin(), rowEnd=rowset.end();
184 rowIdx!= rowEnd; ++rowIdx, r++)
187 for(SIter colIdx = rowset.begin(), colEnd=rowset.end();
188 colIdx!= colEnd; ++colIdx, c++)
190 if (BCRS[*rowIdx].find(*colIdx) == BCRS[*rowIdx].
end())
192 for (
size_t i=0; i<n; i++)
194 for (
size_t j=0; j<n; j++)
196 A[r*n+i][c*n+j] = BCRS[*rowIdx][*colIdx][i][j];
212 template<
class K,
int n,
class Al,
class X,
class Y>
241 void resetIndexForNextDomain();
248 DynamicVector<K> & lhs();
255 DynamicVector<K> & rhs();
268 void operator()(
const size_type& domainIndex);
287 DynamicVector<field_type> * rhs_;
290 DynamicVector<field_type> * lhs_;
298 std::size_t maxlength_;
302 template<
typename T,
typename A,
int n,
int m>
331 void resetIndexForNextDomain();
355 void operator()(
const size_type& domain);
382 std::size_t maxlength_;
387 template<
class M,
class X,
class Y>
470 template<
class M,
class X,
class Y>
489 template<
class M,
class X,
class Y>
507 template<
typename S,
typename T>
512 template<
typename S,
typename T,
typename A,
int n>
518 void operator()(
const size_type& domain);
528 template<
typename S,
typename T>
533 template<
typename S,
typename T,
typename A,
int n>
539 void operator()(
const size_type& domain);
557 template<
typename T,
class X,
class S>
561 template<
class X,
class S>
567 template<
class X,
class S>
573 template<
class X,
class S>
590 template<
typename T1,
typename T2,
bool forward>
618 template<
typename T1,
typename T2>
662 sm.template apply<true>(v, b);
666 template<
class M,
class X,
class TD,
class TA>
674 sm.template apply<true>(v, b);
675 sm.template apply<false>(v, b);
684 template<
class K,
int n,
class Al,
class X,
class Y>
688 template<
class RowToDomain,
class Solvers,
class SubDomains>
689 static std::size_t assembleLocalProblems(
const RowToDomain& rowToDomain,
const matrix_type&
mat,
690 Solvers& solvers,
const SubDomains& domains,
699 template<
class RowToDomain,
class Solvers,
class SubDomains>
700 static std::size_t assembleLocalProblems(
const RowToDomain& rowToDomain,
const matrix_type&
mat,
701 Solvers& solvers,
const SubDomains& domains,
706 template<
class M,
class X,
class Y>
710 template<
class RowToDomain,
class Solvers,
class SubDomains>
712 Solvers& solvers,
const SubDomains& domains,
716 template<
class M,
class X,
class Y>
721 template<
class M,
class X,
class Y>
777 typedef std::set<size_type, std::less<size_type>,
778 typename TA::template rebind<size_type>::other>
782 typedef std::vector<subdomain_type, typename TA::template rebind<subdomain_type>::other>
subdomain_vector;
785 typedef SLList<size_type, typename TA::template rebind<size_type>::other>
subdomain_list;
788 typedef std::vector<subdomain_list, typename TA::template rebind<subdomain_list>::other >
rowtodomain_vector;
794 typedef std::vector<slu, typename TA::template rebind<slu>::other>
slu_vector;
815 field_type relaxationFactor=1,
bool onTheFly_=
true);
829 field_type relaxationFactor=1,
bool onTheFly_=
true);
836 virtual void pre (X& x, X& b) {}
843 virtual void apply (X& v,
const X& d);
845 template<
bool forward>
846 void apply(X& v,
const X& d);
854 Dune::dverb<<
" avg nnz over subdomain is "<<nnz<<std::endl;
863 typename M::size_type maxlength;
871 template<
class I,
class S,
class D>
875 : initializers(&il), indices(&idx), indexMaps(il.size()), domains(domains_)
880 template<
class I,
class S,
class D>
883 typedef typename IndexSet::value_type::const_iterator iterator;
884 for(iterator domain=(*indices)[row.index()].begin(); domain != (*indices)[row.index()].end(); ++domain){
885 (*initializers)[*domain].addRowNnz(row, domains[*domain]);
886 indexMaps[*domain].insert(row.index());
890 template<
class I,
class S,
class D>
893 std::for_each(initializers->begin(), initializers->end(),
894 std::mem_fun_ref(&AtomInitializer::allocateMatrixStorage));
895 std::for_each(initializers->begin(), initializers->end(),
896 std::mem_fun_ref(&AtomInitializer::allocateMarker));
899 template<
class I,
class S,
class D>
902 typedef typename IndexSet::value_type::const_iterator iterator;
903 for(iterator domain=(*indices)[row.index()].begin(); domain != (*indices)[row.index()].end(); ++domain){
904 typename std::map<size_type,size_type>::const_iterator v = indexMaps[*domain].find(col.index());
905 if(v!= indexMaps[*domain].end()){
906 (*initializers)[*domain].countEntries(indexMaps[*domain].find(col.index())->second);
911 template<
class I,
class S,
class D>
914 std::for_each(initializers->begin(), initializers->end(),
915 std::mem_fun_ref(&AtomInitializer::calcColstart));
918 template<
class I,
class S,
class D>
921 typedef typename IndexSet::value_type::const_iterator iterator;
922 for(iterator domain=(*indices)[row.index()].begin(); domain!= (*indices)[row.index()].end(); ++domain){
923 typename std::map<size_type,size_type>::const_iterator v = indexMaps[*domain].find(col.index());
924 if(v!= indexMaps[*domain].end()){
925 assert(indexMaps[*domain].end()!=indexMaps[*domain].find(row.index()));
926 (*initializers)[*domain].copyValue(col, indexMaps[*domain].find(row.index())->second,
932 template<
class I,
class S,
class D>
936 indexMaps.swap(std::vector<IndexMap>(indexMaps));
937 std::for_each(initializers->begin(), initializers->end(),
938 std::mem_fun_ref(&AtomInitializer::createMatrix));
941 template<
class I,
class S,
class D>
946 template<
class I,
class S,
class D>
949 assert(map_.find(grow)==map_.end());
950 map_.insert(std::make_pair(grow,
row++));
953 template<
class I,
class S,
class D>
954 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator
957 return map_.find(grow);
960 template<
class I,
class S,
class D>
961 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator
964 return map_.find(grow);
967 template<
class I,
class S,
class D>
968 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator
974 template<
class I,
class S,
class D>
975 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator
981 template<
class I,
class S,
class D>
982 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator
988 template<
class I,
class S,
class D>
989 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator
995 template<
class M,
class X,
class TM,
class TD,
class TA>
998 :
mat(mat_), relax(relaxationFactor), onTheFly(fly)
1000 typedef typename rowtodomain_vector::const_iterator RowDomainIterator;
1001 typedef typename subdomain_list::const_iterator DomainIterator;
1002 #ifdef DUNE_ISTL_WITH_CHECKING
1003 assert(rowToDomain.size()==mat.N());
1004 assert(rowToDomain.size()==mat.M());
1006 for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter)
1007 assert(iter->size()>0);
1012 for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter)
1013 for(DomainIterator d=iter->begin(); d != iter->end(); ++d)
1014 domains=std::max(domains, *d);
1017 solvers.resize(domains);
1018 subDomains.resize(domains);
1022 for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter, ++
row)
1023 for(DomainIterator d=iter->begin(); d != iter->end(); ++d)
1024 subDomains[*d].insert(row);
1026 #ifdef DUNE_ISTL_WITH_CHECKING
1028 typedef typename subdomain_vector::const_iterator iterator;
1029 for(iterator iter=subDomains.begin(); iter != subDomains.end(); ++iter){
1030 typedef typename subdomain_type::const_iterator entry_iterator;
1031 Dune::dvverb<<
"domain "<<i++<<
":";
1032 for(entry_iterator entry = iter->begin(); entry != iter->end(); ++entry){
1033 Dune::dvverb<<
" "<<*entry;
1035 Dune::dvverb<<std::endl;
1042 template<
class M,
class X,
class TM,
class TD,
class TA>
1047 :
mat(mat_), solvers(sd.size()), subDomains(sd), relax(relaxationFactor),
1050 typedef typename subdomain_vector::const_iterator DomainIterator;
1052 #ifdef DUNE_ISTL_WITH_CHECKING
1055 for(DomainIterator d=sd.begin(); d != sd.end(); ++d,++i){
1057 assert(d->size()>0);
1058 typedef typename DomainIterator::value_type::const_iterator entry_iterator;
1059 Dune::dvverb<<
"domain "<<i<<
":";
1060 for(entry_iterator entry = d->begin(); entry != d->end(); ++entry){
1061 Dune::dvverb<<
" "<<*entry;
1063 Dune::dvverb<<std::endl;
1073 for(DomainIterator domain=sd.begin(); domain != sd.end(); ++domain, ++domainId){
1074 typedef typename subdomain_type::const_iterator iterator;
1075 for(iterator
row=domain->begin();
row != domain->end(); ++
row)
1076 rowToDomain[*
row].push_back(domainId);
1092 template<
typename T,
typename A,
int n,
int m>
1095 template<
class Domain>
1096 static int size(
const Domain & d)
1103 template<
class K,
int n,
class Al,
class X,
class Y>
1104 template<
class RowToDomain,
class Solvers,
class SubDomains>
1106 SeqOverlappingSchwarzAssembler< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >::
1107 assembleLocalProblems(
const RowToDomain& rowToDomain,
1110 const SubDomains& subDomains,
1113 typedef typename SubDomains::const_iterator DomainIterator;
1114 std::size_t maxlength = 0;
1118 for(DomainIterator domain=subDomains.begin();domain!=subDomains.end();++domain)
1119 maxlength=std::max(maxlength, domain->size());
1127 template<
class RowToDomain,
class Solvers,
class SubDomains>
1131 const SubDomains& subDomains,
1134 typedef typename std::vector<SuperMatrixInitializer<matrix_type> >::iterator InitializerIterator;
1135 typedef typename SubDomains::const_iterator DomainIterator;
1136 typedef typename Solvers::iterator SolverIterator;
1137 std::size_t maxlength = 0;
1140 for(DomainIterator domain=subDomains.begin();domain!=subDomains.end();++domain)
1141 maxlength=std::max(maxlength, domain->size());
1142 maxlength*=mat[0].begin()->N();
1145 DomainIterator domain=subDomains.begin();
1148 std::vector<SuperMatrixInitializer<matrix_type> > initializers(subDomains.size());
1150 SolverIterator solver=solvers.begin();
1151 for(InitializerIterator initializer=initializers.begin(); initializer!=initializers.end();
1152 ++initializer, ++solver, ++domain){
1161 RowToDomain, SubDomains> Initializer;
1163 Initializer initializer(initializers, rowToDomain, subDomains);
1165 if(solvers.size()==1)
1166 assert(solvers[0].mat==mat);
1173 for(SolverIterator solver=solvers.begin(); solver!=solvers.end(); ++solver){
1174 assert(solver->mat.N()==solver->mat.M());
1175 maxlength=std::max(maxlength, solver->mat.N());
1184 template<
class M,
class X,
class Y>
1185 template<
class RowToDomain,
class Solvers,
class SubDomains>
1189 const SubDomains& subDomains,
1192 typedef typename SubDomains::const_iterator DomainIterator;
1193 typedef typename Solvers::iterator SolverIterator;
1194 std::size_t maxlength = 0;
1197 for(DomainIterator domain=subDomains.begin();domain!=subDomains.end();++domain)
1198 maxlength=std::max(maxlength, domain->size());
1201 SolverIterator solver=solvers.begin();
1202 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end();
1203 ++domain, ++solver){
1204 solver->setSubMatrix(mat, *domain);
1205 maxlength=std::max(maxlength, domain->size());
1214 template<
class M,
class X,
class TM,
class TD,
class TA>
1220 template<
class M,
class X,
class TM,
class TD,
class TA>
1221 template<
bool forward>
1224 typedef typename X::block_type block;
1225 typedef slu_vector solver_vector;
1238 Adder adder(v, x, assigner, relax);
1244 std::for_each(domain->begin(), domain->end(), assigner);
1245 assigner.resetIndexForNextDomain();
1249 sdsolver.setSubMatrix(
mat, *domain);
1251 sdsolver.apply(assigner.lhs(), assigner.rhs());
1254 solver->apply(assigner.lhs(), assigner.rhs());
1260 std::for_each(domain->begin(), domain->end(), adder);
1261 assigner.resetIndexForNextDomain();
1267 assigner.deallocate();
1270 template<
class K,
int n,
class Al,
class X,
class Y>
1271 OverlappingAssigner< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
1273 const X& b_, Y& x_) :
1275 rhs_( new DynamicVector<
field_type>(maxlength, 42) ),
1276 lhs_( new DynamicVector<
field_type>(maxlength, -42) ),
1280 maxlength_(maxlength)
1283 template<
class K,
int n,
class Al,
class X,
class Y>
1292 template<
class K,
int n,
class Al,
class X,
class Y>
1295 ::resetIndexForNextDomain()
1300 template<
class K,
int n,
class Al,
class X,
class Y>
1308 template<
class K,
int n,
class Al,
class X,
class Y>
1316 template<
class K,
int n,
class Al,
class X,
class Y>
1324 template<
class K,
int n,
class Al,
class X,
class Y>
1333 assert(i<maxlength_);
1334 rhs()[i]=(*b)[domainIndex][j];
1341 for(col_iterator
col=(*
mat)[domainIndex].begin();
col!=(*mat)[domainIndex].end(); ++
col){
1343 (*col).mv((*x)[
col.index()], tmp);
1346 assert(i<maxlength_);
1353 assert(i<maxlength_);
1354 rhs()[i]=(*b)[domainIndex][j];
1360 for(col_iterator
col=(*
mat)[domainIndex].begin();
col!=(*mat)[domainIndex].end(); ++
col){
1362 rhs()[i]-=(*col)[j][k] * (*x)[
col.index()][k];
1369 template<
class K,
int n,
class Al,
class X,
class Y>
1376 assert(i<maxlength_);
1383 template<
typename T,
typename A,
int n,
int m>
1391 x(&x_), i(0), maxlength_(maxlength)
1398 template<
typename T,
typename A,
int n,
int m>
1405 template<
typename T,
typename A,
int n,
int m>
1412 assert(i<maxlength_);
1413 rhs_[i]=(*b)[domainIndex][j];
1421 for(col_iterator
col=(*
mat)[domainIndex].begin();
col!=(*mat)[domainIndex].end(); ++
col){
1423 (*col).mv((*x)[
col.index()], tmp);
1426 assert(i<maxlength_);
1433 template<
typename T,
typename A,
int n,
int m>
1437 assert(i<maxlength_);
1443 template<
typename T,
typename A,
int n,
int m>
1448 assert(i<maxlength_);
1453 template<
typename T,
typename A,
int n,
int m>
1459 template<
typename T,
typename A,
int n,
int m>
1466 template<
typename T,
typename A,
int n,
int m>
1475 template<
class M,
class X,
class Y>
1484 rhs_=
new Y(maxlength);
1485 lhs_ =
new X(maxlength);
1488 template<
class M,
class X,
class Y>
1495 template<
class M,
class X,
class Y>
1498 (*rhs_)[i]=(*b)[domainIndex];
1501 typedef typename matrix_type::ConstColIterator col_iterator;
1504 for(col_iterator
col=(*
mat)[domainIndex].begin();
col!=(*mat)[domainIndex].end(); ++
col){
1505 (*col).mmv((*x)[
col.index()], (*rhs_)[i]);
1511 template<
class M,
class X,
class Y>
1517 template<
class M,
class X,
class Y>
1523 template<
class M,
class X,
class Y>
1529 template<
class M,
class X,
class Y>
1535 template<
class M,
class X,
class Y>
1541 template<
typename S,
typename T,
typename A,
int n>
1546 : v(&v_), x(&x_), assigner(&assigner_), relax(relax_)
1549 template<
typename S,
typename T,
typename A,
int n>
1553 assigner->assignResult((*v)[domainIndex]);
1557 template<
typename S,
typename T,
typename A,
int n>
1565 template<
typename S,
typename T,
typename A,
int n>
1570 : x(&x_), assigner(&assigner_), relax(relax_)
1574 template<
typename S,
typename T,
typename A,
int n>
1578 assigner->relaxResult(relax);
1579 assigner->assignResult((*x)[domainIndex]);
1583 template<
typename S,
typename T,
typename A,
int n>