3 #ifndef DUNE_AMGSMOOTHER_HH
4 #define DUNE_AMGSMOOTHER_HH
11 #include <dune/common/propertymap.hh>
12 #include <dune/common/unused.hh>
69 template<
class X,
class Y,
class C,
class T>
76 template<
class C,
class T>
89 typedef typename T::matrix_type Matrix;
105 DUNE_UNUSED_PARAMETER(amap);
123 DUNE_UNUSED_PARAMETER(comm);
148 template<
class T,
class C=SequentialInformation>
171 class ConstructionTraits;
176 template<
class M,
class X,
class Y,
int l>
198 template<
class M,
class X,
class Y,
int l>
218 template<
class M,
class X,
class Y,
int l>
240 template<
class M,
class X,
class Y>
258 template<
class M,
class X,
class Y>
284 template<
class M,
class X,
class Y>
305 template<
class M,
class X,
class Y,
class C>
322 template<
class X,
class Y,
class C,
class T>
341 template<
class C,
class T>
374 typedef typename Smoother::range_type
Range;
375 typedef typename Smoother::domain_type
Domain;
407 template<
typename LevelContext>
408 void presmooth(LevelContext& levelContext,
size_t steps)
410 for(std::size_t i=0; i < steps; ++i) {
413 ::preSmooth(*levelContext.smoother, *levelContext.lhs,
416 *levelContext.update += *levelContext.lhs;
419 levelContext.matrix->applyscaleadd(-1, *levelContext.lhs, *levelContext.rhs);
420 levelContext.pinfo->project(*levelContext.rhs);
429 template<
typename LevelContext>
432 for(std::size_t i=0; i < steps; ++i) {
434 levelContext.matrix->applyscaleadd(-1, *levelContext.lhs,
437 levelContext.pinfo->project(*levelContext.rhs);
439 ::postSmooth(*levelContext.smoother, *levelContext.lhs, *levelContext.rhs);
441 *levelContext.update += *levelContext.lhs;
445 template<
class M,
class X,
class Y,
int l>
454 smoother.template apply<true>(v,d);
460 smoother.template apply<false>(v,d);
464 template<
class M,
class X,
class Y,
class C,
int l>
473 smoother.template apply<true>(v,d);
479 smoother.template apply<false>(v,d);
483 template<
class M,
class X,
class Y,
class C,
int l>
492 smoother.template apply<true>(v,d);
498 smoother.template apply<false>(v,d);
505 template<
class M,
class X,
class MO,
class MS,
class A>
506 class SeqOverlappingSchwarz;
508 struct MultiplicativeSchwarzMode;
512 template<
class M,
class X,
class MS,
class TA>
522 smoother.template apply<true>(v,d);
528 smoother.template apply<false>(v,d);
546 bool onthefly_=
false)
551 template<
class M,
class X,
class TM,
class TS,
class TA>
557 template<
class M,
class X,
class TM,
class TS,
class TA>
574 std::vector<bool> visited(amap.
noVertices(),
false);
575 typedef IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap> VisitedMapType;
576 VisitedMapType visitedMap(visited.begin());
582 switch(Father::getArgs().overlap) {
583 case SmootherArgs::vertex :
585 VertexAdder visitor(subdomains, amap);
586 createSubdomains(matrix, graph, amap, visitor, visitedMap);
589 case SmootherArgs::pairwise :
591 createPairDomains(graph);
594 case SmootherArgs::aggregate :
596 AggregateAdder<VisitedMapType> visitor(subdomains, amap, graph, visitedMap);
597 createSubdomains(matrix, graph, amap, visitor, visitedMap);
600 case SmootherArgs::none :
602 createSubdomains(matrix, graph, amap, visitor, visitedMap);
605 DUNE_THROW(NotImplemented,
"This overlapping scheme is not supported!");
616 iter!=amap.end(); ++iter)
619 std::vector<bool> visited(amap.noVertices(),
false);
620 typedef IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap> VisitedMapType;
621 VisitedMapType visitedMap(visited.begin());
627 switch(Father::getArgs().overlap) {
628 case SmootherArgs::vertex :
630 VertexAdder visitor(subdomains, amap);
631 createSubdomains(matrix, graph, amap, visitor, visitedMap);
634 case SmootherArgs::aggregate :
636 DUNE_THROW(NotImplemented,
"Aggregate overlap is not supported yet");
643 case SmootherArgs::pairwise :
645 createPairDomains(graph);
648 case SmootherArgs::none :
650 createSubdomains(matrix, graph, amap, visitor, visitedMap);
664 : subdomains(subdomains_), max(-1), subdomain(-1), aggregates(aggregates_)
667 void operator()(
const T& edge)
670 subdomains[subdomain].insert(edge.target());
674 subdomain=aggregate_;
675 max = std::max(subdomain, aggregate_);
678 int noSubdomains()
const
691 void operator()(
const T& edge)
697 int noSubdomains()
const
704 struct AggregateAdder
708 : subdomains(subdomains_), subdomain(-1), aggregates(aggregates_),
709 adder(subdomains_, aggregates_), graph(graph_), visitedMap(visitedMap_)
712 void operator()(
const T& edge)
714 subdomains[subdomain].insert(edge.target());
719 assert(aggregates[edge.target()]!=aggregate);
721 aggregates.template breadthFirstSearch<true,false>(edge.target(), aggregate,
722 graph, vlist, adder, adder,
729 adder.setAggregate(aggregate_);
730 aggregate=aggregate_;
733 int noSubdomains()
const
752 typedef typename M::size_type size_type;
754 std::set<std::pair<size_type,size_type> > pairs;
756 for(VIter v=graph.
begin(), ve=graph.
end(); ve != v; ++v)
757 for(EIter e = v.begin(), ee=v.end(); ee!=e; ++e)
760 if(e.source()<e.target())
761 pairs.insert(std::make_pair(e.source(),e.target()));
763 pairs.insert(std::make_pair(e.target(),e.source()));
767 subdomains.resize(pairs.size());
768 Dune::dinfo <<std::endl<<
"Created "<<pairs.size()<<
" ("<<total<<
") pair domains"<<std::endl<<std::endl;
769 typedef typename std::set<std::pair<size_type,size_type> >::const_iterator SIter;
770 typename Vector::iterator subdomain=subdomains.begin();
772 for(SIter s=pairs.begin(), se =pairs.end(); se!=s; ++s)
774 subdomain->insert(s->first);
775 subdomain->insert(s->second);
778 std::size_t minsize=10000;
779 std::size_t maxsize=0;
781 for(
typename Vector::size_type i=0; i < subdomains.size(); ++i) {
782 sum+=subdomains[i].size();
783 minsize=std::min(minsize, subdomains[i].size());
784 maxsize=std::max(maxsize, subdomains[i].size());
786 Dune::dinfo<<
"Subdomain size: min="<<minsize<<
" max="<<maxsize<<
" avg="<<(sum/subdomains.size())
787 <<
" no="<<subdomains.size()<<std::endl;
790 template<
class Visitor>
791 void createSubdomains(
const M& matrix,
const MatrixGraph<const M>& graph,
792 const AggregatesMap& amap, Visitor& overlapVisitor,
793 IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap>& visitedMap )
800 AggregateDescriptor maxAggregate=0;
802 for(std::size_t i=0; i < amap.noVertices(); ++i)
806 maxAggregate = std::max(maxAggregate, amap[i]);
808 subdomains.resize(maxAggregate+1+isolated);
811 for(
typename Vector::size_type i=0; i < subdomains.size(); ++i)
812 subdomains[i].clear();
817 VertexAdder aggregateVisitor(subdomains, amap);
819 for(VertexDescriptor i=0; i < amap.noVertices(); ++i)
820 if(!
get(visitedMap, i)) {
821 AggregateDescriptor aggregate=amap[i];
825 subdomains.push_back(Subdomain());
826 aggregate=subdomains.size()-1;
828 overlapVisitor.setAggregate(aggregate);
829 aggregateVisitor.setAggregate(aggregate);
830 subdomains[aggregate].insert(i);
832 amap.template breadthFirstSearch<false,false>(i, aggregate, graph, vlist, aggregateVisitor,
833 overlapVisitor, visitedMap);
836 std::size_t minsize=10000;
837 std::size_t maxsize=0;
839 for(
typename Vector::size_type i=0; i < subdomains.size(); ++i) {
840 sum+=subdomains[i].size();
841 minsize=std::min(minsize, subdomains[i].size());
842 maxsize=std::max(maxsize, subdomains[i].size());
844 Dune::dinfo<<
"Subdomain size: min="<<minsize<<
" max="<<maxsize<<
" avg="<<(sum/subdomains.size())
845 <<
" no="<<subdomains.size()<<
" isolated="<<isolated<<std::endl;
854 template<
class M,
class X,
class TM,
class TS,
class TA>