3 #ifndef DUNE_AMGSMOOTHER_HH
4 #define DUNE_AMGSMOOTHER_HH
11 #include <dune/common/propertymap.hh>
12 #include <dune/common/ftraits.hh>
69 template<
class X,
class Y>
76 template<
class X,
class Y,
class C,
class T>
81 template<
class C,
class T>
92 typedef typename T::matrix_type Matrix;
148 template<
class T,
class C=SequentialInformation>
170 template<
class X,
class Y>
181 template <
class... Args>
205 const SmootherArgs* args_;
212 struct ConstructionTraits;
217 template<
class M,
class X,
class Y,
int l>
224 return std::make_shared<SeqSSOR<M,X,Y,l>>
233 template<
class M,
class X,
class Y,
int l>
240 return std::make_shared<SeqSOR<M,X,Y,l>>
249 template<
class M,
class X,
class Y,
int l>
256 return std::make_shared<SeqJac<M,X,Y,l>>
264 template<
class X,
class Y>
271 return std::make_shared<Richardson<X,Y>>
277 template<
class M,
class X,
class Y>
304 template<
class M,
class X,
class Y>
311 return std::make_shared<SeqILU<M,X,Y>>
319 template<
class M,
class X,
class Y,
class C>
326 return std::make_shared<ParSSOR<M,X,Y,C>>
332 template<
class X,
class Y,
class C,
class T>
339 auto seqPrec = SeqConstructionTraits::construct(args);
340 return std::make_shared<BlockPreconditioner<X,Y,C,T>> (seqPrec, args.
getComm());
344 template<
class C,
class T>
351 auto seqPrec = SeqConstructionTraits::construct(args);
352 return std::make_shared<NonoverlappingBlockPreconditioner<C,T>> (seqPrec, args.
getComm());
370 typedef typename Smoother::range_type
Range;
371 typedef typename Smoother::domain_type
Domain;
403 template<
typename LevelContext>
404 void presmooth(LevelContext& levelContext,
size_t steps)
406 for(std::size_t i=0; i < steps; ++i) {
409 ::preSmooth(*levelContext.smoother, *levelContext.lhs,
412 *levelContext.update += *levelContext.lhs;
415 levelContext.matrix->applyscaleadd(-1, *levelContext.lhs, *levelContext.rhs);
416 levelContext.pinfo->project(*levelContext.rhs);
425 template<
typename LevelContext>
428 for(std::size_t i=0; i < steps; ++i) {
430 levelContext.matrix->applyscaleadd(-1, *levelContext.lhs,
433 levelContext.pinfo->project(*levelContext.rhs);
435 ::postSmooth(*levelContext.smoother, *levelContext.lhs, *levelContext.rhs);
437 *levelContext.update += *levelContext.lhs;
441 template<
class M,
class X,
class Y,
int l>
450 smoother.template apply<true>(v,d);
456 smoother.template apply<false>(v,d);
460 template<
class M,
class X,
class Y,
class C,
int l>
469 smoother.template apply<true>(v,d);
475 smoother.template apply<false>(v,d);
479 template<
class M,
class X,
class Y,
class C,
int l>
488 smoother.template apply<true>(v,d);
494 smoother.template apply<false>(v,d);
501 template<
class M,
class X,
class MO,
class MS,
class A>
502 class SeqOverlappingSchwarz;
504 struct MultiplicativeSchwarzMode;
508 template<
class M,
class X,
class MS,
class TA>
518 smoother.template apply<true>(v,d);
524 smoother.template apply<false>(v,d);
542 bool onthefly_=
false)
547 template<
class M,
class X,
class TM,
class TS,
class TA>
553 template<
class M,
class X,
class TM,
class TS,
class TA>
568 Father::setMatrix(matrix);
570 std::vector<bool> visited(amap.
noVertices(),
false);
571 typedef IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap> VisitedMapType;
572 VisitedMapType visitedMap(visited.begin());
578 switch(Father::getArgs().overlap) {
579 case SmootherArgs::vertex :
581 VertexAdder visitor(subdomains, amap);
582 createSubdomains(matrix, graph, amap, visitor, visitedMap);
585 case SmootherArgs::pairwise :
587 createPairDomains(graph);
590 case SmootherArgs::aggregate :
592 AggregateAdder<VisitedMapType> visitor(subdomains, amap, graph, visitedMap);
593 createSubdomains(matrix, graph, amap, visitor, visitedMap);
596 case SmootherArgs::none :
598 createSubdomains(matrix, graph, amap, visitor, visitedMap);
601 DUNE_THROW(NotImplemented,
"This overlapping scheme is not supported!");
606 Father::setMatrix(matrix);
612 iter!=amap.
end(); ++iter)
615 std::vector<bool> visited(amap.
noVertices(),
false);
616 typedef IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap> VisitedMapType;
617 VisitedMapType visitedMap(visited.begin());
623 switch(Father::getArgs().overlap) {
624 case SmootherArgs::vertex :
626 VertexAdder visitor(subdomains, amap);
627 createSubdomains(matrix, graph, amap, visitor, visitedMap);
630 case SmootherArgs::aggregate :
632 DUNE_THROW(NotImplemented,
"Aggregate overlap is not supported yet");
639 case SmootherArgs::pairwise :
641 createPairDomains(graph);
644 case SmootherArgs::none :
646 createSubdomains(matrix, graph, amap, visitor, visitedMap);
660 : subdomains(subdomains_), max(-1), subdomain(-1), aggregates(aggregates_)
666 subdomains[subdomain].insert(edge.target());
670 subdomain=aggregate_;
671 max = std::max(subdomain, aggregate_);
700 struct AggregateAdder
704 : subdomains(subdomains_), subdomain(-1), aggregates(aggregates_),
705 adder(subdomains_, aggregates_), graph(graph_), visitedMap(visitedMap_)
710 subdomains[subdomain].insert(edge.target());
715 assert(aggregates[edge.target()]!=aggregate);
717 aggregates.template breadthFirstSearch<true,false>(edge.target(), aggregate,
718 graph, vlist, adder, adder,
725 adder.setAggregate(aggregate_);
726 aggregate=aggregate_;
748 typedef typename M::size_type size_type;
750 std::set<std::pair<size_type,size_type> > pairs;
752 for(VIter v=graph.
begin(), ve=graph.
end(); ve != v; ++v)
753 for(EIter e = v.begin(), ee=v.end(); ee!=e; ++e)
756 if(e.source()<e.target())
757 pairs.insert(std::make_pair(e.source(),e.target()));
759 pairs.insert(std::make_pair(e.target(),e.source()));
763 subdomains.resize(pairs.size());
764 Dune::dinfo <<std::endl<<
"Created "<<pairs.size()<<
" ("<<total<<
") pair domains"<<std::endl<<std::endl;
765 typedef typename std::set<std::pair<size_type,size_type> >::const_iterator SIter;
766 typename Vector::iterator subdomain=subdomains.
begin();
768 for(SIter s=pairs.begin(), se =pairs.end(); se!=s; ++s)
770 subdomain->insert(s->first);
771 subdomain->insert(s->second);
774 std::size_t minsize=10000;
775 std::size_t maxsize=0;
777 for(
typename Vector::size_type i=0; i < subdomains.size(); ++i) {
778 sum+=subdomains[i].size();
779 minsize=std::min(minsize, subdomains[i].size());
780 maxsize=std::max(maxsize, subdomains[i].size());
782 Dune::dinfo<<
"Subdomain size: min="<<minsize<<
" max="<<maxsize<<
" avg="<<(sum/subdomains.size())
783 <<
" no="<<subdomains.size()<<std::endl;
786 template<
class Visitor>
787 void createSubdomains(
const M& matrix,
const MatrixGraph<const M>& graph,
788 const AggregatesMap& amap, Visitor& overlapVisitor,
789 IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap>& visitedMap )
796 AggregateDescriptor maxAggregate=0;
798 for(std::size_t i=0; i < amap.noVertices(); ++i)
802 maxAggregate = std::max(maxAggregate, amap[i]);
804 subdomains.resize(maxAggregate+1+isolated);
807 for(
typename Vector::size_type i=0; i < subdomains.size(); ++i)
808 subdomains[i].clear();
813 VertexAdder aggregateVisitor(subdomains, amap);
815 for(VertexDescriptor i=0; i < amap.noVertices(); ++i)
816 if(!
get(visitedMap, i)) {
817 AggregateDescriptor aggregate=amap[i];
821 subdomains.push_back(Subdomain());
822 aggregate=subdomains.size()-1;
824 overlapVisitor.setAggregate(aggregate);
825 aggregateVisitor.setAggregate(aggregate);
826 subdomains[aggregate].insert(i);
828 amap.template breadthFirstSearch<false,false>(i, aggregate, graph, vlist, aggregateVisitor,
829 overlapVisitor, visitedMap);
832 std::size_t minsize=10000;
833 std::size_t maxsize=0;
835 for(
typename Vector::size_type i=0; i < subdomains.size(); ++i) {
836 sum+=subdomains[i].size();
837 minsize=std::min(minsize, subdomains[i].size());
838 maxsize=std::max(maxsize, subdomains[i].size());
840 Dune::dinfo<<
"Subdomain size: min="<<minsize<<
" max="<<maxsize<<
" avg="<<(sum/subdomains.size())
841 <<
" no="<<subdomains.size()<<
" isolated="<<isolated<<std::endl;
850 template<
class M,
class X,
class TM,
class TS,
class TA>
857 return std::make_shared<SeqOverlappingSchwarz<M,X,TM,TS,TA>>
Provides classes for the Coloring process of AMG.
Helper classes for the construction of classes without empty constructor.
Define general preconditioner interface.
const Matrix & getMatrix() const
Definition: smoother.hh:112
DefaultSmootherArgs< typename X::field_type > Arguments
Definition: smoother.hh:72
static std::shared_ptr< SeqILU< M, X, Y > > construct(Arguments &args)
Definition: smoother.hh:309
static void postSmooth(Smoother &smoother, Domain &v, Range &d)
Definition: smoother.hh:492
void operator()(const T &edge)
Definition: smoother.hh:687
const SequentialInformation & getComm()
Definition: smoother.hh:194
static void preSmooth(Smoother &smoother, Domain &v, const Range &d)
Definition: smoother.hh:516
ConstructionArgs< SeqILU< M, X, Y > > Arguments
Definition: smoother.hh:307
Smoother::range_type Range
Definition: smoother.hh:513
int setAggregate(const AggregateDescriptor &aggregate_)
Definition: smoother.hh:689
DefaultSmootherArgs< typename T::matrix_type::field_type > Arguments
Definition: smoother.hh:65
int setAggregate(const AggregateDescriptor &aggregate_)
Definition: smoother.hh:668
ConstructionTraits< T > SeqConstructionTraits
Definition: smoother.hh:348
void setArgs(const SmootherArgs &args)
Definition: smoother.hh:117
ConstructionArgs< SeqOverlappingSchwarz< M, X, TM, TS, TA > > Arguments
Definition: smoother.hh:853
NonoverlappingBlockPreconditioner< C, SeqSOR< M, X, Y, l > > Smoother
Definition: smoother.hh:482
static void postSmooth(Smoother &smoother, Domain &v, Range &d)
Definition: smoother.hh:454
AggregateAdder(Vector &subdomains_, const AggregatesMap &aggregates_, const MatrixGraph< const M > &graph_, VM &visitedMap_)
Definition: smoother.hh:702
Overlap
Definition: smoother.hh:536
void setComm(const C &comm)
Definition: smoother.hh:156
DefaultConstructionArgs< Richardson< X, Y > > Arguments
Definition: smoother.hh:267
virtual ~DefaultConstructionArgs()
Definition: smoother.hh:178
SeqOverlappingSchwarzSmootherArgs< typename M::field_type > Arguments
Definition: smoother.hh:550
Vector::value_type Subdomain
Definition: smoother.hh:564
static void preSmooth(Smoother &smoother, Domain &v, Range &d)
Definition: smoother.hh:486
int getN()
Definition: smoother.hh:291
Smoother::range_type Range
Definition: smoother.hh:483
SeqSOR< M, X, Y, l > Smoother
Definition: smoother.hh:444
void operator()(const T &edge)
Definition: smoother.hh:708
static std::shared_ptr< SeqSSOR< M, X, Y, l > > construct(Arguments &args)
Definition: smoother.hh:222
AggregatesMap::AggregateDescriptor AggregateDescriptor
Definition: smoother.hh:562
static std::shared_ptr< Richardson< X, Y > > construct(Arguments &args)
Definition: smoother.hh:269
bool onthefly
Definition: smoother.hh:539
DefaultConstructionArgs< SeqSOR< M, X, Y, l > > Arguments
Definition: smoother.hh:236
virtual void setMatrix(const M &matrix, const AggregatesMap &amap)
Definition: smoother.hh:566
void setMatrix(const Args &...)
Definition: smoother.hh:182
static std::shared_ptr< SeqJac< M, X, Y, l > > construct(Arguments &args)
Definition: smoother.hh:254
DefaultConstructionArgs< SeqJac< M, X, Y, l > > Arguments
Definition: smoother.hh:252
DefaultConstructionArgs< SeqSSOR< M, X, Y, l > > Arguments
Definition: smoother.hh:220
static std::shared_ptr< NonoverlappingBlockPreconditioner< C, T > > construct(Arguments &args)
Definition: smoother.hh:349
const SmootherArgs getArgs() const
Definition: smoother.hh:131
VertexAdder(Vector &subdomains_, const AggregatesMap &aggregates_)
Definition: smoother.hh:659
static std::shared_ptr< SeqSOR< M, X, Y, l > > construct(Arguments &args)
Definition: smoother.hh:238
SeqOverlappingSchwarz< M, X, TM, TS, TA >::subdomain_vector Vector
Definition: smoother.hh:563
static void preSmooth(Smoother &smoother, Domain &v, Range &d)
Definition: smoother.hh:467
const_iterator begin() const
Definition: aggregates.hh:723
void setMatrix(const Matrix &matrix)
Definition: smoother.hh:102
T Smoother
Definition: smoother.hh:369
static void preSmooth(Smoother &smoother, Domain &v, Range &d)
Definition: smoother.hh:448
BlockPreconditioner< X, Y, C, SeqSOR< M, X, Y, l > > Smoother
Definition: smoother.hh:463
static std::shared_ptr< BlockPreconditioner< X, Y, C, T > > construct(Arguments &args)
Definition: smoother.hh:337
int noSubdomains() const
Definition: smoother.hh:674
AggregateDescriptor * iterator
Definition: aggregates.hh:733
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition: smoother.hh:404
Smoother::domain_type Domain
Definition: smoother.hh:514
SeqOverlappingSchwarzSmootherArgs(Overlap overlap_=vertex, bool onthefly_=false)
Definition: smoother.hh:541
DefaultParallelConstructionArgs< T, C > Arguments
Definition: smoother.hh:335
DefaultParallelConstructionArgs< T, C > Arguments
Definition: smoother.hh:347
const Matrix * matrix_
Definition: smoother.hh:137
const_iterator end() const
Definition: aggregates.hh:728
V AggregateDescriptor
The aggregate descriptor type.
Definition: aggregates.hh:578
static const V ISOLATED
Identifier of isolated vertices.
Definition: aggregates.hh:569
const SequentialInformation & getComm()
Definition: smoother.hh:126
DefaultSmootherArgs()
Default constructor.
Definition: smoother.hh:54
void setMatrix(const M &matrix)
Definition: smoother.hh:604
void setComm([[maybe_unused]] T1 &comm)
Definition: smoother.hh:191
std::size_t noVertices() const
Get the number of vertices.
void operator()(const T &edge)
Definition: smoother.hh:663
static void postSmooth(Smoother &smoother, Domain &v, const Range &d)
apply post smoothing in forward direction
Definition: smoother.hh:392
FieldTraits< T >::real_type RelaxationFactor
The type of the relaxation factor.
Definition: smoother.hh:40
Smoother::domain_type Domain
Definition: smoother.hh:484
static std::shared_ptr< ParSSOR< M, X, Y, C > > construct(Arguments &args)
Definition: smoother.hh:324
virtual void setMatrix(const Matrix &matrix, [[maybe_unused]] const AggregatesMap &amap)
Definition: smoother.hh:106
Smoother::domain_type Domain
Definition: smoother.hh:371
void setN(int n)
Definition: smoother.hh:286
SLList< VertexDescriptor, Allocator > VertexList
The type of a single linked list of vertex descriptors.
Definition: aggregates.hh:590
SeqOverlappingSchwarz< M, X, MultiplicativeSchwarzMode, MS, TA > Smoother
Definition: smoother.hh:512
MatrixGraph< M >::VertexDescriptor VertexDescriptor
Definition: smoother.hh:560
int noSubdomains() const
Definition: smoother.hh:729
const Vector & getSubDomains()
Definition: smoother.hh:651
const C & getComm() const
Definition: smoother.hh:161
int noSubdomains() const
Definition: smoother.hh:693
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition: smoother.hh:426
Dune::Amg::AggregatesMap< VertexDescriptor > AggregatesMap
Definition: smoother.hh:561
const SmootherArgs getArgs() const
Definition: smoother.hh:199
static void postSmooth(Smoother &smoother, Domain &v, const Range &d)
Definition: smoother.hh:522
DefaultParallelConstructionArgs< M, C > Arguments
Definition: smoother.hh:322
RelaxationFactor relaxationFactor
The relaxation factor to use.
Definition: smoother.hh:49
virtual ~DefaultParallelConstructionArgs()
Definition: smoother.hh:153
static std::shared_ptr< SeqOverlappingSchwarz< M, X, TM, TS, TA > > construct(Arguments &args)
Definition: smoother.hh:855
Smoother::domain_type Domain
Definition: smoother.hh:446
int setAggregate(const AggregateDescriptor &aggregate_)
Definition: smoother.hh:723
Smoother::range_type Range
Definition: smoother.hh:445
virtual ~DefaultConstructionArgs()
Definition: smoother.hh:99
static void preSmooth(Smoother &smoother, Domain &v, const Range &d)
apply pre smoothing in forward direction
Definition: smoother.hh:380
static void postSmooth(Smoother &smoother, Domain &v, Range &d)
Definition: smoother.hh:473
ConstructionTraits< T > SeqConstructionTraits
Definition: smoother.hh:336
ConstructionArgs(int n=0)
Definition: smoother.hh:282
Smoother::range_type Range
Definition: smoother.hh:464
void setComm([[maybe_unused]] T1 &comm)
Definition: smoother.hh:123
void setArgs(const SmootherArgs &args)
Definition: smoother.hh:185
int iterations
The numbe of iterations to perform.
Definition: smoother.hh:45
Smoother::range_type Range
Definition: smoother.hh:370
Overlap overlap
Definition: smoother.hh:538
Smoother::domain_type Domain
Definition: smoother.hh:465
@ aggregate
Definition: smoother.hh:536
@ none
Definition: smoother.hh:536
@ pairwise
Definition: smoother.hh:536
@ vertex
Definition: smoother.hh:536
Definition: allocator.hh:9
PropertyMapTypeSelector< Amg::VertexVisitedTag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > >::Type get([[maybe_unused]] const Amg::VertexVisitedTag &tag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > &graph)
Definition: dependency.hh:291
Sequential overlapping Schwarz preconditioner.
Definition: overlappingschwarz.hh:753
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
Traits class for generically constructing non default constructable types.
Definition: construction.hh:37
Nonoverlapping parallel preconditioner.
Definition: novlpschwarz.hh:277
P::range_type range_type
The range type of the preconditioner.
Definition: novlpschwarz.hh:285
P::domain_type domain_type
The domain type of the preconditioner.
Definition: novlpschwarz.hh:283
Tag that tells the Schwarz method to be multiplicative.
Definition: overlappingschwarz.hh:124
Class providing information about the mapping of the vertices onto aggregates.
Definition: aggregates.hh:558
VertexIterator end()
Get an iterator over the vertices.
M::size_type VertexDescriptor
The vertex descriptor.
Definition: graph.hh:71
VertexIterator begin()
Get an iterator over the vertices.
The default class for the smoother arguments.
Definition: smoother.hh:36
Traits class for getting the attribute class of a smoother.
Definition: smoother.hh:64
Construction Arguments for the default smoothers.
Definition: smoother.hh:91
Definition: smoother.hh:146
Definition: smoother.hh:151
Definition: smoother.hh:172
Definition: smoother.hh:280
Helper class for applying the smoothers.
Definition: smoother.hh:368
Definition: smoother.hh:535
Definition: smoother.hh:556
Sequential SSOR preconditioner.
Definition: preconditioners.hh:139
Sequential SOR preconditioner.
Definition: preconditioners.hh:259
X domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:264
Y range_type
The range type of the preconditioner.
Definition: preconditioners.hh:266
The sequential jacobian preconditioner.
Definition: preconditioners.hh:410
Sequential ILU preconditioner.
Definition: preconditioners.hh:530
Richardson preconditioner.
Definition: preconditioners.hh:711
A parallel SSOR preconditioner.
Definition: schwarz.hh:176
Block parallel preconditioner.
Definition: schwarz.hh:279
X domain_type
The domain type of the preconditioner.
Definition: schwarz.hh:286
Y range_type
The range type of the preconditioner.
Definition: schwarz.hh:291