3#ifndef DUNE_AMGSMOOTHER_HH
4#define DUNE_AMGSMOOTHER_HH
9#include <dune/istl/schwarz.hh>
10#include <dune/istl/novlpschwarz.hh>
11#include <dune/common/propertymap.hh>
69 template<
class X,
class Y,
class C,
class T>
76 template<
class C,
class T>
79 typedef DefaultSmootherArgs<typename T::matrix_type::field_type> Arguments;
89 typedef typename T::matrix_type Matrix;
99 void setMatrix(
const Matrix& matrix)
103 virtual void setMatrix(
const Matrix& matrix,
const AggregatesMap& amap)
110 const Matrix& getMatrix()
const
121 void setComm(T1& comm)
126 const SequentialInformation& getComm()
137 const Matrix* matrix_;
140 SequentialInformation comm_;
144 struct ConstructionArgs
148 template<
class T,
class C=SequentialInformation>
149 class DefaultParallelConstructionArgs
150 :
public ConstructionArgs<T>
153 virtual ~DefaultParallelConstructionArgs()
156 void setComm(
const C& comm)
161 const C& getComm()
const
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>
259 class ConstructionArgs<
SeqILUn<M,X,Y> >
263 ConstructionArgs(
int n=1)
284 template<
class M,
class X,
class Y>
287 typedef ConstructionArgs<SeqILUn<M,X,Y> >
Arguments;
292 args.getArgs().relaxationFactor);
303 template<
class M,
class X,
class Y>
304 class ConstructionArgs<
SeqILU<M,X,Y> >
308 ConstructionArgs(
int n=0)
330 template<
class M,
class X,
class Y>
333 typedef ConstructionArgs<SeqILU<M,X,Y> >
Arguments;
338 args.getArgs().relaxationFactor);
351 template<
class M,
class X,
class Y,
class C>
354 typedef DefaultParallelConstructionArgs<M,C>
Arguments;
359 args.getArgs().relaxationFactor,
368 template<
class X,
class Y,
class C,
class T>
371 typedef DefaultParallelConstructionArgs<T,C>
Arguments;
381 SeqConstructionTraits::deconstruct(
static_cast<T*
>(&bp->preconditioner));
387 template<
class C,
class T>
388 struct ConstructionTraits<NonoverlappingBlockPreconditioner<C,T> >
390 typedef DefaultParallelConstructionArgs<T,C>
Arguments;
391 typedef ConstructionTraits<T> SeqConstructionTraits;
394 return new NonoverlappingBlockPreconditioner<C,T>(*SeqConstructionTraits::construct(args),
398 static inline void deconstruct(NonoverlappingBlockPreconditioner<C,T>* bp)
400 SeqConstructionTraits::deconstruct(
static_cast<T*
>(&bp->preconditioner));
420 typedef typename Smoother::range_type Range;
421 typedef typename Smoother::domain_type Domain;
430 static void preSmooth(Smoother& smoother, Domain& v,
const Range& d)
442 static void postSmooth(Smoother& smoother, Domain& v,
const Range& d)
453 template<
typename LevelContext>
454 void presmooth(LevelContext& levelContext,
size_t steps)
456 for(std::size_t i=0; i < steps; ++i) {
459 ::preSmooth(*levelContext.smoother, *levelContext.lhs,
462 *levelContext.update += *levelContext.lhs;
465 levelContext.matrix->applyscaleadd(-1, *levelContext.lhs, *levelContext.rhs);
466 levelContext.pinfo->project(*levelContext.rhs);
475 template<
typename LevelContext>
478 for(std::size_t i=0; i < steps; ++i) {
480 levelContext.matrix->applyscaleadd(-1, *levelContext.lhs,
483 levelContext.pinfo->project(*levelContext.rhs);
485 ::postSmooth(*levelContext.smoother, *levelContext.lhs, *levelContext.rhs);
487 *levelContext.update += *levelContext.lhs;
491 template<
class M,
class X,
class Y,
int l>
492 struct SmootherApplier<
SeqSOR<M,X,Y,l> >
495 typedef typename Smoother::range_type Range;
496 typedef typename Smoother::domain_type Domain;
498 static void preSmooth(Smoother& smoother, Domain& v, Range& d)
500 smoother.template apply<true>(v,d);
504 static void postSmooth(Smoother& smoother, Domain& v, Range& d)
506 smoother.template apply<false>(v,d);
510 template<
class M,
class X,
class Y,
class C,
int l>
511 struct SmootherApplier<BlockPreconditioner<X,Y,C,SeqSOR<M,X,Y,l> > >
513 typedef BlockPreconditioner<X,Y,C,SeqSOR<M,X,Y,l> > Smoother;
514 typedef typename Smoother::range_type Range;
515 typedef typename Smoother::domain_type Domain;
517 static void preSmooth(Smoother& smoother, Domain& v, Range& d)
519 smoother.template apply<true>(v,d);
523 static void postSmooth(Smoother& smoother, Domain& v, Range& d)
525 smoother.template apply<false>(v,d);
529 template<
class M,
class X,
class Y,
class C,
int l>
530 struct SmootherApplier<NonoverlappingBlockPreconditioner<C,SeqSOR<M,X,Y,l> > >
532 typedef NonoverlappingBlockPreconditioner<C,SeqSOR<M,X,Y,l> > Smoother;
533 typedef typename Smoother::range_type Range;
534 typedef typename Smoother::domain_type Domain;
536 static void preSmooth(Smoother& smoother, Domain& v, Range& d)
538 smoother.template apply<true>(v,d);
542 static void postSmooth(Smoother& smoother, Domain& v, Range& d)
544 smoother.template apply<false>(v,d);
551 template<
class M,
class X,
class MO,
class MS,
class A>
552 class SeqOverlappingSchwarz;
554 struct MultiplicativeSchwarzMode;
558 template<
class M,
class X,
class MS,
class TA>
559 struct SmootherApplier<SeqOverlappingSchwarz<M,X,MultiplicativeSchwarzMode,
562 typedef SeqOverlappingSchwarz<M,X,MultiplicativeSchwarzMode,MS,TA> Smoother;
563 typedef typename Smoother::range_type Range;
564 typedef typename Smoother::domain_type Domain;
566 static void preSmooth(Smoother& smoother, Domain& v,
const Range& d)
568 smoother.template apply<true>(v,d);
572 static void postSmooth(Smoother& smoother, Domain& v,
const Range& d)
574 smoother.template apply<false>(v,d);
583 struct SeqOverlappingSchwarzSmootherArgs
584 :
public DefaultSmootherArgs<T>
591 SeqOverlappingSchwarzSmootherArgs(Overlap overlap_=
vertex,
592 bool onthefly_=
false)
593 : overlap(overlap_), onthefly(onthefly_)
597 template<
class M,
class X,
class TM,
class TS,
class TA>
598 struct SmootherTraits<SeqOverlappingSchwarz<M,X,TM,TS,TA> >
600 typedef SeqOverlappingSchwarzSmootherArgs<typename M::field_type> Arguments;
603 template<
class M,
class X,
class TM,
class TS,
class TA>
604 class ConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> >
605 :
public DefaultConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> >
607 typedef DefaultConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> > Father;
614 typedef typename Vector::value_type Subdomain;
616 virtual void setMatrix(
const M& matrix,
const AggregatesMap& amap)
618 Father::setMatrix(matrix);
620 std::vector<bool> visited(amap.noVertices(),
false);
621 typedef IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap> VisitedMapType;
622 VisitedMapType visitedMap(visited.begin());
624 MatrixGraph<const M> graph(matrix);
626 typedef SeqOverlappingSchwarzSmootherArgs<typename M::field_type> SmootherArgs;
628 switch(Father::getArgs().overlap) {
631 VertexAdder visitor(subdomains, amap);
632 createSubdomains(matrix, graph, amap, visitor, visitedMap);
635 case SmootherArgs::pairwise :
637 createPairDomains(graph);
640 case SmootherArgs::aggregate :
642 AggregateAdder<VisitedMapType> visitor(subdomains, amap, graph, visitedMap);
643 createSubdomains(matrix, graph, amap, visitor, visitedMap);
648 createSubdomains(matrix, graph, amap, visitor, visitedMap);
651 DUNE_THROW(NotImplemented,
"This overlapping scheme is not supported!");
654 void setMatrix(
const M& matrix)
656 Father::setMatrix(matrix);
659 AggregatesMap amap(matrix.N());
660 VertexDescriptor v=0;
661 for(
typename AggregatesMap::iterator iter=amap.begin();
662 iter!=amap.end(); ++iter)
665 std::vector<bool> visited(amap.noVertices(),
false);
666 typedef IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap> VisitedMapType;
667 VisitedMapType visitedMap(visited.begin());
669 MatrixGraph<const M> graph(matrix);
671 typedef SeqOverlappingSchwarzSmootherArgs<typename M::field_type> SmootherArgs;
673 switch(Father::getArgs().overlap) {
676 VertexAdder visitor(subdomains, amap);
677 createSubdomains(matrix, graph, amap, visitor, visitedMap);
680 case SmootherArgs::aggregate :
682 DUNE_THROW(NotImplemented,
"Aggregate overlap is not supported yet");
689 case SmootherArgs::pairwise :
691 createPairDomains(graph);
696 createSubdomains(matrix, graph, amap, visitor, visitedMap);
701 const Vector& getSubDomains()
709 VertexAdder(Vector& subdomains_,
const AggregatesMap& aggregates_)
710 : subdomains(subdomains_), max(-1), subdomain(-1), aggregates(aggregates_)
713 void operator()(
const T& edge)
716 subdomains[subdomain].insert(edge.target());
718 int setAggregate(
const AggregateDescriptor& aggregate_)
720 subdomain=aggregate_;
721 max = std::max(subdomain, aggregate_);
724 int noSubdomains()
const
730 AggregateDescriptor max;
731 AggregateDescriptor subdomain;
732 const AggregatesMap& aggregates;
737 void operator()(
const T& edge)
739 int setAggregate(
const AggregateDescriptor& aggregate_)
743 int noSubdomains()
const
750 struct AggregateAdder
752 AggregateAdder(Vector& subdomains_,
const AggregatesMap& aggregates_,
753 const MatrixGraph<const M>& graph_, VM& visitedMap_)
754 : subdomains(subdomains_), subdomain(-1), aggregates(aggregates_),
755 adder(subdomains_, aggregates_), graph(graph_), visitedMap(visitedMap_)
758 void operator()(
const T& edge)
760 subdomains[subdomain].insert(edge.target());
765 assert(aggregates[edge.target()]!=aggregate);
767 aggregates.template breadthFirstSearch<true,false>(edge.target(), aggregate,
768 graph, vlist, adder, adder,
773 int setAggregate(
const AggregateDescriptor& aggregate_)
775 adder.setAggregate(aggregate_);
776 aggregate=aggregate_;
779 int noSubdomains()
const
785 AggregateDescriptor aggregate;
788 const AggregatesMap& aggregates;
790 const MatrixGraph<const M>& graph;
794 void createPairDomains(
const MatrixGraph<const M>& graph)
798 typedef typename M::size_type size_type;
800 std::set<std::pair<size_type,size_type> > pairs;
802 for(VIter v=graph.begin(), ve=graph.end(); ve != v; ++v)
803 for(EIter e = v.begin(), ee=v.end(); ee!=e; ++e)
806 if(e.source()<e.target())
807 pairs.insert(std::make_pair(e.source(),e.target()));
809 pairs.insert(std::make_pair(e.target(),e.source()));
813 subdomains.resize(pairs.size());
814 Dune::dinfo <<std::endl<<
"Created "<<pairs.size()<<
" ("<<total<<
") pair domains"<<std::endl<<std::endl;
815 typedef typename std::set<std::pair<size_type,size_type> >::const_iterator SIter;
816 typename Vector::iterator subdomain=subdomains.begin();
818 for(SIter s=pairs.begin(), se =pairs.end(); se!=s; ++s)
820 subdomain->insert(s->first);
821 subdomain->insert(s->second);
824 std::size_t minsize=10000;
825 std::size_t maxsize=0;
827 for(
typename Vector::size_type i=0; i < subdomains.size(); ++i) {
828 sum+=subdomains[i].size();
829 minsize=std::min(minsize, subdomains[i].size());
830 maxsize=std::max(maxsize, subdomains[i].size());
832 Dune::dinfo<<
"Subdomain size: min="<<minsize<<
" max="<<maxsize<<
" avg="<<(sum/subdomains.size())
833 <<
" no="<<subdomains.size()<<std::endl;
836 template<
class Visitor>
837 void createSubdomains(
const M& matrix,
const MatrixGraph<const M>& graph,
838 const AggregatesMap& amap, Visitor& overlapVisitor,
839 IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap>& visitedMap )
846 AggregateDescriptor maxAggregate=0;
848 for(std::size_t i=0; i < amap.noVertices(); ++i)
852 maxAggregate = std::max(maxAggregate, amap[i]);
854 subdomains.resize(maxAggregate+1+isolated);
857 for(
typename Vector::size_type i=0; i < subdomains.size(); ++i)
858 subdomains[i].clear();
863 VertexAdder aggregateVisitor(subdomains, amap);
865 for(VertexDescriptor i=0; i < amap.noVertices(); ++i)
866 if(!get(visitedMap, i)) {
867 AggregateDescriptor aggregate=amap[i];
871 subdomains.push_back(Subdomain());
872 aggregate=subdomains.size()-1;
874 overlapVisitor.setAggregate(aggregate);
875 aggregateVisitor.setAggregate(aggregate);
876 subdomains[aggregate].insert(i);
878 amap.template breadthFirstSearch<false,false>(i, aggregate, graph, vlist, aggregateVisitor,
879 overlapVisitor, visitedMap);
882 std::size_t minsize=10000;
883 std::size_t maxsize=0;
885 for(
typename Vector::size_type i=0; i < subdomains.size(); ++i) {
886 sum+=subdomains[i].size();
887 minsize=std::min(minsize, subdomains[i].size());
888 maxsize=std::max(maxsize, subdomains[i].size());
890 Dune::dinfo<<
"Subdomain size: min="<<minsize<<
" max="<<maxsize<<
" avg="<<(sum/subdomains.size())
891 <<
" no="<<subdomains.size()<<
" isolated="<<isolated<<std::endl;
900 template<
class M,
class X,
class TM,
class TS,
class TA>
901 struct ConstructionTraits<SeqOverlappingSchwarz<M,X,TM,TS,TA> >
903 typedef ConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> >
Arguments;
907 return new SeqOverlappingSchwarz<M,X,TM,TS,TA>(args.getMatrix(),
908 args.getSubDomains(),
909 args.getArgs().relaxationFactor,
910 args.getArgs().onthefly);
913 static void deconstruct(SeqOverlappingSchwarz<M,X,TM,TS,TA>* schwarz)
Provides classes for the Coloring process of AMG.
Class providing information about the mapping of the vertices onto aggregates.
Definition: aggregates.hh:543
Traits class for generically constructing non default constructable types.
Definition: construction.hh:38
Construction Arguments for the default smoothers.
Definition: smoother.hh:88
VertexIteratorT< const MatrixGraph< Matrix > > ConstVertexIterator
The constant vertex iterator type.
Definition: graph.hh:306
M::size_type VertexDescriptor
The vertex descriptor.
Definition: graph.hh:71
EdgeIteratorT< const MatrixGraph< Matrix > > ConstEdgeIterator
The constant edge iterator type.
Definition: graph.hh:296
Block parallel preconditioner.
Definition: schwarz.hh:269
Nonoverlapping parallel preconditioner.
Definition: novlpschwarz.hh:266
Sequential ILU0 preconditioner.
Definition: preconditioners.hh:652
Sequential ILU preconditioner.
Definition: preconditioners.hh:504
Sequential ILU(n) preconditioner.
Definition: preconditioners.hh:739
The sequential jacobian preconditioner.
Definition: preconditioners.hh:416
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:784
Sequential SOR preconditioner.
Definition: preconditioners.hh:225
Sequential SSOR preconditioner.
Definition: preconditioners.hh:135
Helper classes for the construction of classes without empty constructor.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:25
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
constexpr GeometryType none(unsigned int dim)
Returns a GeometryType representing a singular of dimension dim.
Definition: type.hh:714
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:727
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition: smoother.hh:454
V AggregateDescriptor
The aggregate descriptor type.
Definition: aggregates.hh:563
static const V ISOLATED
Identifier of isolated vertices.
Definition: aggregates.hh:554
DefaultSmootherArgs()
Default constructor.
Definition: smoother.hh:54
const void * Arguments
A type holding all the arguments needed to call the constructor.
Definition: construction.hh:44
static void postSmooth(Smoother &smoother, Domain &v, const Range &d)
apply post smoothing in forward direction
Definition: smoother.hh:442
SLList< VertexDescriptor, Allocator > VertexList
The type of a single linked list of vertex descriptors.
Definition: aggregates.hh:575
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition: smoother.hh:476
static T * construct(Arguments &args)
Construct an object with the specified arguments.
Definition: construction.hh:52
RelaxationFactor relaxationFactor
The relaxation factor to use.
Definition: smoother.hh:49
static void preSmooth(Smoother &smoother, Domain &v, const Range &d)
apply pre smoothing in forward direction
Definition: smoother.hh:430
static void deconstruct(T *t)
Destroys an object.
Definition: construction.hh:61
int iterations
The numbe of iterations to perform.
Definition: smoother.hh:45
T RelaxationFactor
The type of the relaxation factor.
Definition: smoother.hh:40
DInfoType dinfo(std::cout)
Stream for informative output.
Definition: stdstreams.hh:138
PartitionSet<... > Overlap
Type of PartitionSet for the overlap partition.
Definition: partitionset.hh:248
Dune namespace.
Definition: alignedallocator.hh:10
Define general preconditioner interface.
The default class for the smoother arguments.
Definition: smoother.hh:36
Helper class for applying the smoothers.
Definition: smoother.hh:418
Traits class for getting the attribute class of a smoother.
Definition: smoother.hh:64
Definition of the DUNE_UNUSED macro for the case that config.h is not available.