1 #ifndef DUNE_AMGSMOOTHER_HH
2 #define DUNE_AMGSMOOTHER_HH
9 #include<dune/common/propertymap.hh>
66 template<
class X,
class Y,
class C,
class T>
73 template<
class C,
class T>
86 typedef typename T::matrix_type Matrix;
136 template<
class T,
class C=SequentialInformation>
159 class ConstructionTraits;
164 template<
class M,
class X,
class Y,
int l>
186 template<
class M,
class X,
class Y,
int l>
206 template<
class M,
class X,
class Y,
int l>
228 template<
class M,
class X,
class Y>
246 template<
class M,
class X,
class Y>
272 template<
class M,
class X,
class Y>
293 template<
class M,
class X,
class Y,
class C>
310 template<
class X,
class Y,
class C,
class T>
323 SeqConstructionTraits::deconstruct(static_cast<T*>(&bp->preconditioner));
329 template<
class C,
class T>
342 SeqConstructionTraits::deconstruct(static_cast<T*>(&bp->preconditioner));
362 typedef typename Smoother::range_type
Range;
363 typedef typename Smoother::domain_type
Domain;
390 template<
class M,
class X,
class Y,
int l>
399 smoother.template apply<true>(v,d);
405 smoother.template apply<false>(v,d);
409 template<
class M,
class X,
class Y,
class C,
int l>
418 smoother.template apply<true>(v,d);
424 smoother.template apply<false>(v,d);
428 template<
class M,
class X,
class Y,
class C,
int l>
437 smoother.template apply<true>(v,d);
443 smoother.template apply<false>(v,d);
450 template<
class M,
class X,
class MO,
class MS,
class A>
451 class SeqOverlappingSchwarz;
453 class MultiplicativeSchwarzMode;
457 template<
class M,
class X,
class MS,
class TA>
467 smoother.template apply<true>(v,d);
473 smoother.template apply<false>(v,d);
491 bool onthefly_=
false)
496 template<
class M,
class X,
class TM,
class TS,
class TA>
502 template<
class M,
class X,
class TM,
class TS,
class TA>
517 Father::setMatrix(matrix);
519 std::vector<bool> visited(amap.
noVertices(),
false);
520 typedef IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap> VisitedMapType;
521 VisitedMapType visitedMap(visited.begin());
527 switch(Father::getArgs().overlap){
528 case SmootherArgs::vertex:
530 VertexAdder visitor(subdomains, amap);
531 createSubdomains(matrix, graph, amap, visitor, visitedMap);
534 case SmootherArgs::pairwise:
536 createPairDomains(graph);
539 case SmootherArgs::aggregate:
541 AggregateAdder<VisitedMapType> visitor(subdomains, amap, graph, visitedMap);
542 createSubdomains(matrix, graph, amap, visitor, visitedMap);
545 case SmootherArgs::none:
547 createSubdomains(matrix, graph, amap, visitor, visitedMap);
550 DUNE_THROW(NotImplemented,
"This overlapping scheme is not supported!");
555 Father::setMatrix(matrix);
561 iter!=amap.end(); ++iter)
564 std::vector<bool> visited(amap.noVertices(),
false);
565 typedef IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap> VisitedMapType;
566 VisitedMapType visitedMap(visited.begin());
572 switch(Father::getArgs().overlap){
573 case SmootherArgs::vertex:
575 VertexAdder visitor(subdomains, amap);
576 createSubdomains(matrix, graph, amap, visitor, visitedMap);
579 case SmootherArgs::aggregate:
581 DUNE_THROW(NotImplemented,
"Aggregate overlap is not supported yet");
588 case SmootherArgs::pairwise:
590 createPairDomains(graph);
593 case SmootherArgs::none:
595 createSubdomains(matrix, graph, amap, visitor, visitedMap);
609 : subdomains(subdomains_), max(-1), subdomain(-1), aggregates(aggregates_)
612 void operator()(
const T& edge)
615 subdomains[subdomain].insert(edge.target());
619 subdomain=aggregate_;
620 max = std::max(subdomain, aggregate_);
623 int noSubdomains()
const
636 void operator()(
const T& edge)
642 int noSubdomains()
const
649 struct AggregateAdder
653 : subdomains(subdomains_), subdomain(-1), aggregates(aggregates_),
654 adder(subdomains_, aggregates_), graph(graph_), visitedMap(visitedMap_)
657 void operator()(
const T& edge)
659 subdomains[subdomain].insert(edge.target());
664 assert(aggregates[edge.target()]!=aggregate);
666 aggregates.template breadthFirstSearch<true,false>(edge.target(), aggregate,
667 graph, vlist, adder, adder,
674 adder.setAggregate(aggregate_);
675 aggregate=aggregate_;
678 int noSubdomains()
const
697 typedef typename M::size_type size_type;
699 std::set<std::pair<size_type,size_type> > pairs;
701 for(VIter v=graph.
begin(), ve=graph.
end(); ve != v; ++v)
702 for(EIter e = v.begin(), ee=v.end(); ee!=e; ++e)
705 if(e.source()<e.target())
706 pairs.insert(std::make_pair(e.source(),e.target()));
708 pairs.insert(std::make_pair(e.target(),e.source()));
712 subdomains.resize(pairs.size());
713 Dune::dinfo <<std::endl<<
"Created "<<pairs.size()<<
" ("<<total<<
") pair domains"<<std::endl<<std::endl;
714 typedef typename std::set<std::pair<size_type,size_type> >::const_iterator SIter;
715 typename Vector::iterator subdomain=subdomains.begin();
717 for(SIter s=pairs.begin(), se =pairs.end(); se!=s; ++s)
719 subdomain->insert(s->first);
720 subdomain->insert(s->second);
723 std::size_t minsize=10000;
724 std::size_t maxsize=0;
726 for(
typename Vector::size_type i=0; i < subdomains.size(); ++i){
727 sum+=subdomains[i].size();
728 minsize=std::min(minsize, subdomains[i].size());
729 maxsize=std::max(maxsize, subdomains[i].size());
731 Dune::dinfo<<
"Subdomain size: min="<<minsize<<
" max="<<maxsize<<
" avg="<<(sum/subdomains.size())
732 <<
" no="<<subdomains.size()<<std::endl;
735 template<
class Visitor>
736 void createSubdomains(
const M& matrix,
const MatrixGraph<const M>& graph,
737 const AggregatesMap& amap, Visitor& overlapVisitor,
738 IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap>& visitedMap )
745 AggregateDescriptor maxAggregate=0;
747 for(std::size_t i=0; i < amap.noVertices(); ++i)
751 maxAggregate = std::max(maxAggregate, amap[i]);
753 subdomains.resize(maxAggregate+1+isolated);
756 for(
typename Vector::size_type i=0; i < subdomains.size(); ++i)
757 subdomains[i].clear();
762 VertexAdder aggregateVisitor(subdomains, amap);
764 for(VertexDescriptor i=0; i < amap.noVertices(); ++i)
765 if(!
get(visitedMap, i)){
766 AggregateDescriptor aggregate=amap[i];
770 subdomains.push_back(Subdomain());
771 aggregate=subdomains.size()-1;
773 overlapVisitor.setAggregate(aggregate);
774 aggregateVisitor.setAggregate(aggregate);
775 subdomains[aggregate].insert(i);
777 amap.template breadthFirstSearch<false,false>(i, aggregate, graph, vlist, aggregateVisitor,
778 overlapVisitor, visitedMap);
781 std::size_t minsize=10000;
782 std::size_t maxsize=0;
784 for(
typename Vector::size_type i=0; i < subdomains.size(); ++i){
785 sum+=subdomains[i].size();
786 minsize=std::min(minsize, subdomains[i].size());
787 maxsize=std::max(maxsize, subdomains[i].size());
789 Dune::dinfo<<
"Subdomain size: min="<<minsize<<
" max="<<maxsize<<
" avg="<<(sum/subdomains.size())
790 <<
" no="<<subdomains.size()<<
" isolated="<<isolated<<std::endl;
799 template<
class M,
class X,
class TM,
class TS,
class TA>