dune-istl  2.2.0
smoother.hh
Go to the documentation of this file.
00001 #ifndef DUNE_AMGSMOOTHER_HH
00002 #define DUNE_AMGSMOOTHER_HH
00003 
00004 #include<dune/istl/paamg/construction.hh>
00005 #include<dune/istl/paamg/aggregates.hh>
00006 #include<dune/istl/preconditioners.hh>
00007 #include<dune/istl/schwarz.hh>
00008 #include<dune/istl/novlpschwarz.hh>
00009 #include<dune/common/propertymap.hh>
00010 
00011 namespace Dune
00012 {
00013   namespace Amg
00014   {
00015     
00031     template<class T>
00032     struct DefaultSmootherArgs
00033     {
00037       typedef T RelaxationFactor;
00038       
00042       int iterations;
00046       RelaxationFactor relaxationFactor;
00047       
00051       DefaultSmootherArgs()
00052         : iterations(1), relaxationFactor(1.0)
00053       {}
00054     };
00055     
00059     template<class T>
00060     struct SmootherTraits
00061     {
00062       typedef DefaultSmootherArgs<typename T::matrix_type::field_type> Arguments;
00063       
00064     };
00065     
00066     template<class X, class Y, class C, class T>
00067     struct SmootherTraits<BlockPreconditioner<X,Y,C,T> >
00068     {
00069       typedef DefaultSmootherArgs<typename T::matrix_type::field_type> Arguments;
00070       
00071     };
00072 
00073    template<class C, class T>
00074    struct SmootherTraits<NonoverlappingBlockPreconditioner<C,T> >
00075    {
00076      typedef DefaultSmootherArgs<typename T::matrix_type::field_type> Arguments;
00077      
00078    };
00079     
00083     template<class T>
00084     class DefaultConstructionArgs
00085     {      
00086       typedef typename T::matrix_type Matrix;
00087       
00088       typedef typename SmootherTraits<T>::Arguments SmootherArgs;
00089       
00090       typedef Dune::Amg::AggregatesMap<typename MatrixGraph<Matrix>::VertexDescriptor> AggregatesMap;
00091       
00092     public:
00093       virtual ~DefaultConstructionArgs()
00094       {}
00095       
00096       void setMatrix(const Matrix& matrix)
00097       {
00098         matrix_=&matrix;
00099       }
00100       virtual void setMatrix(const Matrix& matrix, const AggregatesMap& amap)
00101       {
00102         setMatrix(matrix);
00103       }
00104       
00105      
00106       const Matrix& getMatrix() const
00107       {
00108         return *matrix_;
00109       }
00110       
00111       void setArgs(const SmootherArgs& args)
00112       {
00113         args_=&args;
00114       }
00115 
00116       template<class T1>
00117       void setComm(T1& comm)
00118       {}
00119       
00120       const SmootherArgs getArgs() const
00121       {
00122         return *args_;
00123       }
00124       
00125     protected:
00126       const Matrix* matrix_;
00127     private:
00128       const SmootherArgs* args_;
00129     };
00130 
00131     template<class T>
00132     struct ConstructionArgs
00133       : public DefaultConstructionArgs<T>
00134     {};
00135 
00136     template<class T, class C=SequentialInformation>
00137     class DefaultParallelConstructionArgs 
00138       : public ConstructionArgs<T>
00139     {
00140     public:
00141       virtual ~DefaultParallelConstructionArgs()
00142       {}
00143       
00144       void setComm(const C& comm)
00145       {
00146         comm_ = &comm;
00147       }
00148 
00149       const C& getComm() const
00150       {
00151         return *comm_;
00152       }
00153     private:
00154       const C* comm_;
00155     };
00156     
00157 
00158     template<class T>
00159     class ConstructionTraits;
00160 
00164     template<class M, class X, class Y, int l>
00165     struct ConstructionTraits<SeqSSOR<M,X,Y,l> >
00166     {
00167       typedef DefaultConstructionArgs<SeqSSOR<M,X,Y,l> > Arguments;
00168       
00169       static inline SeqSSOR<M,X,Y,l>* construct(Arguments& args)
00170       {
00171         return new SeqSSOR<M,X,Y,l>(args.getMatrix(), args.getArgs().iterations,
00172                                     args.getArgs().relaxationFactor);
00173       }
00174       
00175       static inline void deconstruct(SeqSSOR<M,X,Y,l>* ssor)
00176       {
00177         delete ssor;
00178       }
00179       
00180     };
00181 
00182     
00186     template<class M, class X, class Y, int l>
00187     struct ConstructionTraits<SeqSOR<M,X,Y,l> >
00188     {
00189       typedef DefaultConstructionArgs<SeqSOR<M,X,Y,l> > Arguments;
00190       
00191       static inline SeqSOR<M,X,Y,l>* construct(Arguments& args)
00192       {
00193         return new SeqSOR<M,X,Y,l>(args.getMatrix(), args.getArgs().iterations,
00194                                     args.getArgs().relaxationFactor);
00195       }
00196       
00197       static inline void deconstruct(SeqSOR<M,X,Y,l>* sor)
00198       {
00199         delete sor;
00200       }
00201       
00202     };
00206     template<class M, class X, class Y, int l>
00207     struct ConstructionTraits<SeqJac<M,X,Y,l> >
00208     {
00209       typedef DefaultConstructionArgs<SeqJac<M,X,Y,l> > Arguments;
00210       
00211       static inline SeqJac<M,X,Y,l>* construct(Arguments& args)
00212       {
00213         return new SeqJac<M,X,Y,l>(args.getMatrix(), args.getArgs().iterations,
00214                                  args.getArgs().relaxationFactor);
00215       }
00216       
00217       static void deconstruct(SeqJac<M,X,Y,l>* jac)
00218       {
00219         delete jac;
00220       }
00221       
00222     };
00223 
00224     
00228     template<class M, class X, class Y>
00229     struct ConstructionTraits<SeqILU0<M,X,Y> >
00230     {
00231       typedef DefaultConstructionArgs<SeqILU0<M,X,Y> > Arguments;
00232       
00233       static inline SeqILU0<M,X,Y>* construct(Arguments& args)
00234       {
00235         return new SeqILU0<M,X,Y>(args.getMatrix(),
00236                                  args.getArgs().relaxationFactor);
00237       }
00238       
00239       static void deconstruct(SeqILU0<M,X,Y>* ilu)
00240       {
00241         delete ilu;
00242       }
00243       
00244     };
00245 
00246     template<class M, class X, class Y>
00247     class ConstructionArgs<SeqILUn<M,X,Y> >
00248       : public DefaultConstructionArgs<SeqILUn<M,X,Y> >
00249     {
00250     public:
00251       ConstructionArgs(int n=1)
00252         : n_(n)
00253       {}
00254       
00255       void setN(int n)
00256       {
00257         n_ = n;
00258       }
00259       int getN()
00260       {
00261         return n_;
00262       }
00263       
00264     private:
00265       int n_;
00266     };
00267     
00268     
00272     template<class M, class X, class Y>
00273     struct ConstructionTraits<SeqILUn<M,X,Y> >
00274     {
00275       typedef ConstructionArgs<SeqILUn<M,X,Y> > Arguments;
00276       
00277       static inline SeqILUn<M,X,Y>* construct(Arguments& args)
00278       {
00279         return new SeqILUn<M,X,Y>(args.getMatrix(), args.getN(),
00280                                  args.getArgs().relaxationFactor);
00281       }
00282       
00283       static void deconstruct(SeqILUn<M,X,Y>* ilu)
00284       {
00285         delete ilu;
00286       }
00287       
00288     };
00289     
00293     template<class M, class X, class Y, class C>
00294     struct ConstructionTraits<ParSSOR<M,X,Y,C> >
00295     {
00296       typedef DefaultParallelConstructionArgs<M,C> Arguments;
00297       
00298       static inline ParSSOR<M,X,Y,C>* construct(Arguments& args)
00299       {
00300         return new ParSSOR<M,X,Y,C>(args.getMatrix(), args.getArgs().iterations,
00301                                     args.getArgs().relaxationFactor,
00302                                     args.getComm());
00303       }      
00304       static inline void deconstruct(ParSSOR<M,X,Y,C>* ssor)
00305       {
00306         delete ssor;
00307       }
00308     };
00309 
00310     template<class X, class Y, class C, class T>
00311     struct ConstructionTraits<BlockPreconditioner<X,Y,C,T> >
00312     {
00313       typedef DefaultParallelConstructionArgs<T,C> Arguments;
00314       typedef ConstructionTraits<T> SeqConstructionTraits;
00315       static inline BlockPreconditioner<X,Y,C,T>* construct(Arguments& args)
00316       {
00317         return new BlockPreconditioner<X,Y,C,T>(*SeqConstructionTraits::construct(args),
00318                                                 args.getComm());
00319       }
00320 
00321       static inline void deconstruct(BlockPreconditioner<X,Y,C,T>* bp)
00322       {
00323         SeqConstructionTraits::deconstruct(static_cast<T*>(&bp->preconditioner));
00324         delete bp;
00325       }
00326       
00327     };
00328 
00329     template<class C, class T>
00330     struct ConstructionTraits<NonoverlappingBlockPreconditioner<C,T> >
00331     {
00332       typedef DefaultParallelConstructionArgs<T,C> Arguments;
00333       typedef ConstructionTraits<T> SeqConstructionTraits;
00334       static inline NonoverlappingBlockPreconditioner<C,T>* construct(Arguments& args)
00335       {
00336         return new NonoverlappingBlockPreconditioner<C,T>(*SeqConstructionTraits::construct(args),
00337                                                 args.getComm());
00338       }
00339 
00340       static inline void deconstruct(NonoverlappingBlockPreconditioner<C,T>* bp)
00341       {
00342         SeqConstructionTraits::deconstruct(static_cast<T*>(&bp->preconditioner));
00343         delete bp;
00344       }
00345       
00346     };
00347 
00358     template<class T>
00359     struct SmootherApplier
00360     {
00361       typedef T Smoother;
00362       typedef typename Smoother::range_type Range;
00363       typedef typename Smoother::domain_type Domain;
00364       
00372       static void preSmooth(Smoother& smoother, Domain& v, const Range& d)
00373       {
00374         smoother.apply(v,d);
00375       }
00376 
00384       static void postSmooth(Smoother& smoother, Domain& v, const Range& d)
00385       {
00386         smoother.apply(v,d);
00387       }
00388     };
00389     
00390     template<class M, class X, class Y, int l>
00391     struct SmootherApplier<SeqSOR<M,X,Y,l> >
00392     {
00393       typedef SeqSOR<M,X,Y,l> Smoother;
00394       typedef typename Smoother::range_type Range;
00395       typedef typename Smoother::domain_type Domain;
00396       
00397       static void preSmooth(Smoother& smoother, Domain& v, Range& d)
00398       {
00399         smoother.template apply<true>(v,d);
00400       }
00401 
00402        
00403       static void postSmooth(Smoother& smoother, Domain& v, Range& d)
00404       {
00405         smoother.template apply<false>(v,d);
00406       }
00407     };
00408 
00409     template<class M, class X, class Y, class C, int l>
00410     struct SmootherApplier<BlockPreconditioner<X,Y,C,SeqSOR<M,X,Y,l> > >
00411     {
00412       typedef BlockPreconditioner<X,Y,C,SeqSOR<M,X,Y,l> > Smoother;
00413       typedef typename Smoother::range_type Range;
00414       typedef typename Smoother::domain_type Domain;
00415       
00416       static void preSmooth(Smoother& smoother, Domain& v, Range& d)
00417       {
00418         smoother.template apply<true>(v,d);
00419       }
00420 
00421        
00422       static void postSmooth(Smoother& smoother, Domain& v, Range& d)
00423       {
00424         smoother.template apply<false>(v,d);
00425       }
00426     };
00427 
00428     template<class M, class X, class Y, class C, int l>
00429     struct SmootherApplier<NonoverlappingBlockPreconditioner<C,SeqSOR<M,X,Y,l> > >
00430     {
00431       typedef NonoverlappingBlockPreconditioner<C,SeqSOR<M,X,Y,l> > Smoother;
00432       typedef typename Smoother::range_type Range;
00433       typedef typename Smoother::domain_type Domain;
00434       
00435       static void preSmooth(Smoother& smoother, Domain& v, Range& d)
00436       {
00437         smoother.template apply<true>(v,d);
00438       }
00439 
00440        
00441       static void postSmooth(Smoother& smoother, Domain& v, Range& d)
00442       {
00443         smoother.template apply<false>(v,d);
00444       }
00445     };
00446 
00447   }// end namespace Amg
00448 
00449     // forward declarations
00450   template<class M, class X, class MO, class MS, class A>
00451     class SeqOverlappingSchwarz;
00452 
00453     class MultiplicativeSchwarzMode;
00454 
00455   namespace Amg 
00456   {
00457     template<class M, class X, class MS, class TA>
00458     struct SmootherApplier<SeqOverlappingSchwarz<M,X,MultiplicativeSchwarzMode,
00459                                                  MS,TA> >
00460     {
00461       typedef SeqOverlappingSchwarz<M,X,MultiplicativeSchwarzMode,MS,TA> Smoother;
00462       typedef typename Smoother::range_type Range;
00463       typedef typename Smoother::domain_type Domain;
00464       
00465       static void preSmooth(Smoother& smoother, Domain& v, const Range& d)
00466       {
00467         smoother.template apply<true>(v,d);
00468       }
00469 
00470        
00471       static void postSmooth(Smoother& smoother, Domain& v, const Range& d)
00472       {
00473         smoother.template apply<false>(v,d);
00474       
00475        }
00476     };
00477     
00478     //    template<class M, class X, class TM, class TA>
00479     //    class SeqOverlappingSchwarz;
00480 
00481     template<class T>
00482     struct SeqOverlappingSchwarzSmootherArgs
00483       : public DefaultSmootherArgs<T>
00484     {
00485       enum Overlap {vertex, aggregate, pairwise, none};
00486       
00487       Overlap overlap;
00488       bool onthefly;
00489       
00490       SeqOverlappingSchwarzSmootherArgs(Overlap overlap_=vertex,
00491                                         bool onthefly_=false)
00492         : overlap(overlap_), onthefly(onthefly_)
00493       {}
00494     };  
00495       
00496     template<class M, class X, class TM, class TS, class TA>
00497     struct SmootherTraits<SeqOverlappingSchwarz<M,X,TM,TS,TA> >
00498     {
00499       typedef  SeqOverlappingSchwarzSmootherArgs<typename M::field_type> Arguments;
00500     };
00501     
00502     template<class M, class X, class TM, class TS, class TA>
00503     class ConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> >
00504       : public DefaultConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> >
00505     {
00506       typedef DefaultConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> > Father;
00507 
00508     public:
00509       typedef typename MatrixGraph<M>::VertexDescriptor VertexDescriptor;
00510       typedef Dune::Amg::AggregatesMap<VertexDescriptor> AggregatesMap;
00511       typedef typename AggregatesMap::AggregateDescriptor AggregateDescriptor;
00512       typedef typename SeqOverlappingSchwarz<M,X,TM,TS,TA>::subdomain_vector Vector;
00513       typedef typename Vector::value_type Subdomain;
00514       
00515       virtual void setMatrix(const M& matrix, const AggregatesMap& amap)
00516       {
00517         Father::setMatrix(matrix);
00518         
00519         std::vector<bool> visited(amap.noVertices(), false);
00520         typedef IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap> VisitedMapType;
00521         VisitedMapType visitedMap(visited.begin());
00522         
00523         MatrixGraph<const M> graph(matrix);
00524         
00525         typedef SeqOverlappingSchwarzSmootherArgs<typename M::field_type> SmootherArgs;
00526         
00527         switch(Father::getArgs().overlap){
00528         case SmootherArgs::vertex:
00529           {  
00530           VertexAdder visitor(subdomains, amap);
00531           createSubdomains(matrix, graph, amap, visitor,  visitedMap);
00532           }
00533           break;
00534         case SmootherArgs::pairwise:
00535           {
00536             createPairDomains(graph);
00537           }
00538           break;
00539         case SmootherArgs::aggregate:
00540           {
00541           AggregateAdder<VisitedMapType> visitor(subdomains, amap, graph, visitedMap);
00542           createSubdomains(matrix, graph, amap, visitor, visitedMap);
00543           }
00544           break;
00545         case SmootherArgs::none:
00546           NoneAdder visitor;
00547           createSubdomains(matrix, graph, amap, visitor, visitedMap);
00548           break;
00549         default:
00550           DUNE_THROW(NotImplemented, "This overlapping scheme is not supported!");
00551         }
00552       }
00553       void setMatrix(const M& matrix)
00554       {
00555         Father::setMatrix(matrix);
00556         
00557         /* Create aggregates map where each aggregate is just one vertex. */
00558         AggregatesMap amap(matrix.N());
00559         VertexDescriptor v=0;
00560         for(typename AggregatesMap::iterator iter=amap.begin();
00561             iter!=amap.end(); ++iter)
00562           *iter=v++;
00563         
00564         std::vector<bool> visited(amap.noVertices(), false);
00565         typedef IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap> VisitedMapType;
00566         VisitedMapType visitedMap(visited.begin());
00567         
00568         MatrixGraph<const M> graph(matrix);
00569         
00570         typedef SeqOverlappingSchwarzSmootherArgs<typename M::field_type> SmootherArgs;
00571         
00572         switch(Father::getArgs().overlap){
00573         case SmootherArgs::vertex:
00574           {  
00575           VertexAdder visitor(subdomains, amap);
00576           createSubdomains(matrix, graph, amap, visitor,  visitedMap);
00577           }
00578           break;
00579         case SmootherArgs::aggregate:
00580           {
00581             DUNE_THROW(NotImplemented, "Aggregate overlap is not supported yet");
00582             /*
00583           AggregateAdder<VisitedMapType> visitor(subdomains, amap, graph, visitedMap);
00584           createSubdomains(matrix, graph, amap, visitor, visitedMap);
00585             */
00586           }
00587           break;
00588         case SmootherArgs::pairwise:
00589           {
00590             createPairDomains(graph);
00591           }
00592           break;
00593         case SmootherArgs::none:
00594           NoneAdder visitor;
00595           createSubdomains(matrix, graph, amap, visitor, visitedMap);
00596           
00597         }
00598       }
00599       
00600       const Vector& getSubDomains()
00601       {
00602         return subdomains;
00603       }
00604       
00605     private:
00606       struct VertexAdder
00607       {
00608         VertexAdder(Vector& subdomains_, const AggregatesMap& aggregates_)
00609           : subdomains(subdomains_), max(-1), subdomain(-1), aggregates(aggregates_)
00610         {}
00611         template<class T>
00612         void operator()(const T& edge)
00613         {
00614           if(aggregates[edge.target()]!=AggregatesMap::ISOLATED)
00615             subdomains[subdomain].insert(edge.target());
00616         }
00617         int setAggregate(const AggregateDescriptor& aggregate_)
00618         {
00619           subdomain=aggregate_;
00620           max = std::max(subdomain, aggregate_);
00621           return subdomain;
00622         }
00623         int noSubdomains() const
00624         {
00625           return max+1;
00626         }
00627       private:
00628         Vector& subdomains;
00629         AggregateDescriptor max;
00630         AggregateDescriptor subdomain;
00631         const AggregatesMap& aggregates;
00632       };
00633       struct NoneAdder
00634       {
00635         template<class T>
00636         void operator()(const T& edge)
00637         {}
00638         int setAggregate(const AggregateDescriptor& aggregate_)
00639         {
00640         return -1;
00641         }
00642         int noSubdomains() const
00643         {
00644           return -1;
00645         }
00646       };
00647 
00648       template<class VM>
00649       struct AggregateAdder
00650       {
00651         AggregateAdder(Vector& subdomains_, const AggregatesMap& aggregates_, 
00652                        const MatrixGraph<const M>& graph_, VM& visitedMap_)
00653           : subdomains(subdomains_), subdomain(-1), aggregates(aggregates_),
00654             adder(subdomains_, aggregates_), graph(graph_), visitedMap(visitedMap_)
00655         {}
00656         template<class T>
00657         void operator()(const T& edge)
00658         {
00659           subdomains[subdomain].insert(edge.target());
00660           // If we (the neighbouring vertex of the aggregate)
00661           // are not isolated, add the aggregate we belong to 
00662           // to the same subdomain using the OneOverlapAdder
00663           if(aggregates[edge.target()]!=AggregatesMap::ISOLATED){
00664             assert(aggregates[edge.target()]!=aggregate);
00665             typename AggregatesMap::VertexList vlist;       
00666             aggregates.template breadthFirstSearch<true,false>(edge.target(), aggregate,
00667                                                                graph, vlist, adder, adder,
00668                                                          visitedMap);
00669           }
00670         }
00671         
00672         int setAggregate(const AggregateDescriptor& aggregate_)
00673         {
00674           adder.setAggregate(aggregate_);
00675           aggregate=aggregate_;
00676           return ++subdomain;
00677         }
00678         int noSubdomains() const
00679         {
00680           return subdomain+1;
00681         }
00682         
00683       private:
00684         AggregateDescriptor aggregate;
00685         Vector& subdomains;
00686         int subdomain;
00687         const AggregatesMap& aggregates;
00688         VertexAdder adder;
00689         const MatrixGraph<const M>& graph;
00690         VM& visitedMap;
00691       };
00692       
00693       void createPairDomains(const MatrixGraph<const M>& graph)
00694       {
00695         typedef typename MatrixGraph<const M>::ConstVertexIterator VIter;
00696         typedef typename MatrixGraph<const M>::ConstEdgeIterator EIter;
00697         typedef typename M::size_type size_type;
00698         
00699         std::set<std::pair<size_type,size_type> > pairs;
00700         int total=0;
00701         for(VIter v=graph.begin(), ve=graph.end(); ve != v; ++v)
00702           for(EIter e = v.begin(), ee=v.end(); ee!=e; ++e)
00703             {
00704               ++total;
00705             if(e.source()<e.target())
00706               pairs.insert(std::make_pair(e.source(),e.target()));
00707             else
00708               pairs.insert(std::make_pair(e.target(),e.source()));
00709             }
00710         
00711         
00712         subdomains.resize(pairs.size());
00713         Dune::dinfo <<std::endl<< "Created "<<pairs.size()<<" ("<<total<<") pair domains"<<std::endl<<std::endl;
00714         typedef typename std::set<std::pair<size_type,size_type> >::const_iterator SIter;
00715         typename Vector::iterator subdomain=subdomains.begin();
00716         
00717         for(SIter s=pairs.begin(), se =pairs.end(); se!=s; ++s)
00718           {
00719             subdomain->insert(s->first);
00720             subdomain->insert(s->second);
00721             ++subdomain;
00722           }
00723         std::size_t minsize=10000;
00724         std::size_t maxsize=0;
00725         int sum=0;
00726         for(typename Vector::size_type i=0; i < subdomains.size(); ++i){
00727           sum+=subdomains[i].size();
00728           minsize=std::min(minsize, subdomains[i].size());
00729           maxsize=std::max(maxsize, subdomains[i].size());
00730         }
00731         Dune::dinfo<<"Subdomain size: min="<<minsize<<" max="<<maxsize<<" avg="<<(sum/subdomains.size())
00732                  <<" no="<<subdomains.size()<<std::endl;
00733       }
00734       
00735       template<class Visitor>
00736       void createSubdomains(const M& matrix, const MatrixGraph<const M>& graph, 
00737                             const AggregatesMap& amap, Visitor& overlapVisitor, 
00738                             IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap>& visitedMap )
00739       {
00740         // count  number ag aggregates. We asume that the
00741         // aggregates are numbered consecutively from 0 exept
00742         // for the isolated ones. All isolated vertices form 
00743         // one aggregate, here.
00744         int isolated=0;
00745         AggregateDescriptor maxAggregate=0;
00746         
00747         for(std::size_t i=0; i < amap.noVertices(); ++i)
00748           if(amap[i]==AggregatesMap::ISOLATED)
00749             isolated++;
00750           else
00751             maxAggregate = std::max(maxAggregate, amap[i]);
00752         
00753         subdomains.resize(maxAggregate+1+isolated);
00754 
00755         // reset the subdomains
00756         for(typename Vector::size_type i=0; i < subdomains.size(); ++i)
00757           subdomains[i].clear();
00758         
00759         // Create the subdomains from the aggregates mapping.
00760         // For each aggregate we mark all entries and the 
00761         // neighbouring vertices as belonging to the same subdomain
00762         VertexAdder aggregateVisitor(subdomains, amap);
00763                 
00764         for(VertexDescriptor i=0; i < amap.noVertices(); ++i)
00765           if(!get(visitedMap, i)){
00766             AggregateDescriptor aggregate=amap[i];
00767 
00768             if(amap[i]==AggregatesMap::ISOLATED){
00769               // isolated vertex gets its own aggregate
00770               subdomains.push_back(Subdomain());
00771               aggregate=subdomains.size()-1;
00772             }
00773             overlapVisitor.setAggregate(aggregate);
00774             aggregateVisitor.setAggregate(aggregate);
00775             subdomains[aggregate].insert(i);
00776             typename AggregatesMap::VertexList vlist;       
00777             amap.template breadthFirstSearch<false,false>(i, aggregate, graph, vlist, aggregateVisitor, 
00778                             overlapVisitor, visitedMap);
00779           }
00780         
00781         std::size_t minsize=10000;
00782         std::size_t maxsize=0;
00783         int sum=0;
00784         for(typename Vector::size_type i=0; i < subdomains.size(); ++i){
00785           sum+=subdomains[i].size();
00786           minsize=std::min(minsize, subdomains[i].size());
00787           maxsize=std::max(maxsize, subdomains[i].size());
00788         }
00789         Dune::dinfo<<"Subdomain size: min="<<minsize<<" max="<<maxsize<<" avg="<<(sum/subdomains.size())
00790                  <<" no="<<subdomains.size()<<" isolated="<<isolated<<std::endl;
00791         
00792           
00793           
00794       }
00795       Vector subdomains;
00796     };
00797     
00798 
00799     template<class M, class X, class TM, class TS, class TA>
00800     struct ConstructionTraits<SeqOverlappingSchwarz<M,X,TM,TS,TA> >
00801     {
00802       typedef ConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> > Arguments;
00803       
00804       static inline SeqOverlappingSchwarz<M,X,TM,TS,TA>* construct(Arguments& args)
00805       {
00806         return new SeqOverlappingSchwarz<M,X,TM,TS,TA>(args.getMatrix(), 
00807                                                       args.getSubDomains(),
00808                                                       args.getArgs().relaxationFactor,
00809                                                       args.getArgs().onthefly);
00810       }
00811       
00812       static void deconstruct(SeqOverlappingSchwarz<M,X,TM,TS,TA>* schwarz)
00813       {
00814         delete schwarz;
00815       }
00816     };
00817 
00818 
00819   } // namespace Amg
00820 } // namespace Dune
00821 
00822 
00823 
00824 #endif