dune-istl  2.2.0
hierarchy.hh
Go to the documentation of this file.
00001 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
00002 // vi: set et ts=8 sw=2 sts=2:
00003 // $Id: hierarchy.hh 1610 2012-05-29 16:01:01Z sander $
00004 #ifndef DUNE_AMGHIERARCHY_HH
00005 #define DUNE_AMGHIERARCHY_HH
00006 
00007 #include<list>
00008 #include<memory>
00009 #include<limits>
00010 #include<algorithm>
00011 #include"aggregates.hh"
00012 #include"graph.hh"
00013 #include"galerkin.hh"
00014 #include"renumberer.hh"
00015 #include"graphcreator.hh"
00016 #include<dune/common/stdstreams.hh>
00017 #include<dune/common/timer.hh>
00018 #include<dune/common/tuples.hh>
00019 #include<dune/common/bigunsignedint.hh>
00020 #include<dune/istl/bvector.hh>
00021 #include<dune/common/parallel/indexset.hh>
00022 #include<dune/istl/matrixutils.hh>
00023 #include<dune/istl/matrixredistribute.hh>
00024 #include<dune/istl/paamg/dependency.hh>
00025 #include<dune/istl/paamg/graph.hh>
00026 #include<dune/istl/paamg/indicescoarsener.hh>
00027 #include<dune/istl/paamg/globalaggregates.hh>
00028 #include<dune/istl/paamg/construction.hh>
00029 #include<dune/istl/paamg/smoother.hh>
00030 #include<dune/istl/paamg/transfer.hh>
00031 
00032 namespace Dune
00033 {
00034   namespace Amg
00035   {
00047     enum{ 
00055       MAX_PROCESSES = 72000};
00056 
00064     template<typename T, typename A=std::allocator<T> >
00065     class Hierarchy
00066     {
00067     public:
00071       typedef T MemberType;
00072 
00073       template<typename T1, typename T2>
00074       class LevelIterator;
00075 
00076     private:
00080       struct Element
00081       {
00082         friend class LevelIterator<Hierarchy<T,A>, T>;
00083         friend class LevelIterator<const Hierarchy<T,A>, const T>;
00084 
00086         Element* coarser_;
00087         
00089         Element* finer_;
00090         
00092         MemberType* element_;
00093         
00095         MemberType* redistributed_;
00096       };
00097     public:
00098 //       enum{
00099 //      /**
00100 //       * @brief If true only the method addCoarser will be usable
00101 //       * otherwise only the method addFiner will be usable.
00102 //       */
00103 //      coarsen = b
00104 //        };
00105       
00109       typedef typename A::template rebind<Element>::other Allocator;
00110 
00111       typedef typename ConstructionTraits<T>::Arguments Arguments;
00112       
00117       Hierarchy(MemberType& first);
00118       
00122       Hierarchy();
00123 
00128       void addCoarser(Arguments& args);
00129 
00130       void addRedistributedOnCoarsest(Arguments& args);
00131       
00136       void addFiner(Arguments& args);
00137 
00144       template<class C, class T1>
00145       class LevelIterator 
00146         : public BidirectionalIteratorFacade<LevelIterator<C,T1>,T1,T1&>
00147       {
00148         friend class LevelIterator<typename remove_const<C>::type,
00149                                    typename remove_const<T1>::type >;
00150         friend class LevelIterator<const typename remove_const<C>::type, 
00151                                    const typename remove_const<T1>::type >;
00152 
00153       public:
00155         LevelIterator()
00156           : element_(0)
00157         {}
00158         
00159         LevelIterator(Element* element)
00160           : element_(element)
00161         {}
00162         
00164         LevelIterator(const LevelIterator<typename remove_const<C>::type,
00165                       typename remove_const<T1>::type>& other)
00166           : element_(other.element_)
00167         {}
00168         
00170         LevelIterator(const LevelIterator<const typename remove_const<C>::type,
00171                       const typename remove_const<T1>::type>& other)
00172           : element_(other.element_)
00173         {}
00174 
00178         bool equals(const LevelIterator<typename remove_const<C>::type, 
00179                     typename remove_const<T1>::type>& other) const
00180         {
00181           return element_ == other.element_;
00182         }
00183         
00187         bool equals(const LevelIterator<const typename remove_const<C>::type, 
00188                     const typename remove_const<T1>::type>& other) const
00189         {
00190           return element_ == other.element_;
00191         }
00192 
00194         T1& dereference() const
00195         {
00196           return *(element_->element_);
00197         }
00198         
00200         void increment()
00201         {
00202           element_ = element_->coarser_;
00203         }
00204         
00206         void decrement()
00207         {
00208           element_ = element_->finer_;
00209         }
00210         
00215         bool isRedistributed() const
00216         {
00217           return element_->redistributed_;
00218         }
00219         
00224         T1& getRedistributed() const
00225         {
00226           assert(element_->redistributed_);
00227           return *element_->redistributed_;
00228         }
00229         void addRedistributed(T1* t)
00230         {
00231           element_->redistributed_ = t;
00232         }
00233 
00234         void deleteRedistributed()
00235         {
00236           element_->redistributed_ = 0;
00237         }
00238         
00239       private:
00240         Element* element_;
00241       };
00242       
00244       typedef LevelIterator<Hierarchy<T,A>,T> Iterator;
00245 
00247       typedef LevelIterator<const Hierarchy<T,A>, const T> ConstIterator;
00248 
00253       Iterator finest();
00254       
00259       Iterator coarsest();
00260       
00261       
00266       ConstIterator finest() const;
00267       
00272       ConstIterator coarsest() const;
00273 
00278       std::size_t levels() const;
00279       
00281       ~Hierarchy();
00282       
00283     private:
00285       Element* finest_;
00287       Element* coarsest_;
00289       Element* nonAllocated_;
00291       Allocator allocator_;
00293       int levels_;
00294     };
00295     
00302     template<class M, class PI, class A=std::allocator<M> >
00303     class MatrixHierarchy
00304     {
00305     public:
00307       typedef M MatrixOperator;
00308       
00310       typedef typename MatrixOperator::matrix_type Matrix;
00311 
00313       typedef PI ParallelInformation;
00314       
00316       typedef A Allocator;
00317 
00319       typedef Dune::Amg::AggregatesMap<typename MatrixGraph<Matrix>::VertexDescriptor> AggregatesMap;
00320 
00322       typedef Dune::Amg::Hierarchy<MatrixOperator,Allocator> ParallelMatrixHierarchy;
00323 
00325       typedef Dune::Amg::Hierarchy<ParallelInformation,Allocator> ParallelInformationHierarchy;
00326 
00328       typedef typename Allocator::template rebind<AggregatesMap*>::other AAllocator;
00329       
00331       typedef std::list<AggregatesMap*,AAllocator> AggregatesMapList;
00332 
00334       typedef RedistributeInformation<ParallelInformation> RedistributeInfoType;
00335       
00337       typedef typename Allocator::template rebind<RedistributeInfoType>::other RILAllocator;
00338       
00340       typedef std::list<RedistributeInfoType,RILAllocator> RedistributeInfoList;
00341 
00347       MatrixHierarchy(const MatrixOperator& fineMatrix,
00348                       const ParallelInformation& pinfo=ParallelInformation());
00349       
00350       
00351       ~MatrixHierarchy();
00352       
00358       template<typename O, typename T>
00359       void build(const T& criterion);
00360 
00368       template<class F>
00369       void recalculateGalerkin(const F& copyFlags);
00370 
00375       template<class V, class TA>
00376       void coarsenVector(Hierarchy<BlockVector<V,TA> >& hierarchy) const;
00377 
00383       template<class S, class TA>
00384       void coarsenSmoother(Hierarchy<S,TA>& smoothers, 
00385                            const typename SmootherTraits<S>::Arguments& args) const;
00386       
00391       std::size_t levels() const;
00392       
00397       std::size_t maxlevels() const;
00398       
00399       bool hasCoarsest() const;
00400       
00405       bool isBuilt() const;
00406 
00411       const ParallelMatrixHierarchy& matrices() const;
00412       
00417       const ParallelInformationHierarchy& parallelInformation() const;
00418       
00423       const AggregatesMapList& aggregatesMaps() const;
00424       
00430       const RedistributeInfoList& redistributeInformation() const;
00431       
00432 
00433       typename MatrixOperator::field_type getProlongationDampingFactor() const
00434       {
00435         return prolongDamp_;
00436       }
00437       
00448       void getCoarsestAggregatesOnFinest(std::vector<std::size_t>& data) const;
00449       
00450     private:
00451       typedef typename ConstructionTraits<MatrixOperator>::Arguments MatrixArgs;
00452       typedef typename ConstructionTraits<ParallelInformation>::Arguments CommunicationArgs;
00454       AggregatesMapList aggregatesMaps_;
00456       RedistributeInfoList redistributes_;
00458       ParallelMatrixHierarchy matrices_;
00460       ParallelInformationHierarchy parallelInformation_;
00461 
00463       bool built_;
00464 
00466       int maxlevels_;
00467       
00468       typename MatrixOperator::field_type prolongDamp_;
00469       
00473       template<class Matrix, bool print>
00474       struct MatrixStats
00475       {
00476         
00480         static void stats(const Matrix& matrix)
00481         {}
00482       };
00483       
00484       template<class Matrix>
00485       struct MatrixStats<Matrix,true>
00486       {
00487         struct calc
00488         {
00489           typedef typename Matrix::size_type size_type;
00490           typedef typename Matrix::row_type matrix_row;
00491           
00492           calc()
00493           {
00494             min=std::numeric_limits<size_type>::max();
00495             max=0;
00496             sum=0;
00497           }
00498           
00499           void operator()(const matrix_row& row)
00500           {
00501             min=std::min(min, row.size());
00502             max=std::max(max, row.size());
00503             sum += row.size();
00504           }
00505           
00506           size_type min;
00507           size_type max;
00508           size_type sum;
00509         };
00513         static void stats(const Matrix& matrix)
00514         {
00515           calc c= for_each(matrix.begin(), matrix.end(), calc());
00516           dinfo<<"Matrix row: min="<<c.min<<" max="<<c.max
00517                <<" average="<<static_cast<double>(c.sum)/matrix.N()
00518                <<std::endl;
00519         }
00520       };
00521     };
00522 
00526     template<class T>
00527     class CoarsenCriterion : public T
00528     {
00529     public:
00534       typedef T AggregationCriterion;
00535       
00546       CoarsenCriterion(int maxLevel=100, int coarsenTarget=1000, double minCoarsenRate=1.2,
00547                        double prolongDamp=1.6, AccumulationMode accumulate=successiveAccu)
00548         : AggregationCriterion(Dune::Amg::Parameters(maxLevel, coarsenTarget, minCoarsenRate, prolongDamp, accumulate))
00549       {}
00550 
00551       CoarsenCriterion(const Dune::Amg::Parameters& parms)
00552         : AggregationCriterion(parms)
00553       {}
00554       
00555     };
00556     
00557     template<typename M, typename C1>
00558     bool repartitionAndDistributeMatrix(const M& origMatrix, M& newMatrix, 
00559                                         SequentialInformation& origSequentialInformationomm, 
00560                                         SequentialInformation*& newComm, 
00561                                         RedistributeInformation<SequentialInformation>& ri,
00562                                         int nparts, C1& criterion)
00563     {
00564       DUNE_THROW(NotImplemented, "Redistribution does not make sense in sequential code!");
00565     }
00566     
00567     
00568     template<typename M, typename C, typename C1>
00569     bool repartitionAndDistributeMatrix(const M& origMatrix, M& newMatrix, C& origComm, C*& newComm, 
00570                                         RedistributeInformation<C>& ri,
00571                                         int nparts, C1& criterion)
00572   {
00573     Timer time;
00574 #ifdef AMG_REPART_ON_COMM_GRAPH
00575     // Done not repartition the matrix graph, but a graph of the communication scheme.
00576     bool existentOnRedist=Dune::commGraphRepartition(origMatrix, origComm, nparts, newComm, 
00577                                                      ri.getInterface(),
00578                                                      criterion.debugLevel()>1);
00579     
00580 #else
00581     typedef Dune::Amg::MatrixGraph<const M> MatrixGraph;
00582     typedef Dune::Amg::PropertiesGraph<MatrixGraph,
00583       VertexProperties,
00584       EdgeProperties,
00585       IdentityMap,
00586       IdentityMap> PropertiesGraph;
00587     MatrixGraph graph(origMatrix);
00588     PropertiesGraph pgraph(graph);
00589     buildDependency(pgraph, origMatrix, criterion, false);
00590     
00591 #ifdef DEBUG_REPART
00592     if(origComm.communicator().rank()==0)
00593       std::cout<<"Original matrix"<<std::endl;
00594     origComm.communicator().barrier();
00595   printGlobalSparseMatrix(origMatrix, origComm, std::cout);
00596 #endif
00597     bool existentOnRedist=Dune::graphRepartition(pgraph, origComm, nparts,
00598                                                  newComm, ri.getInterface(),
00599                                                  criterion.debugLevel()>1);
00600 #endif // if else AMG_REPART
00601     
00602     if(origComm.communicator().rank()==0  && criterion.debugLevel()>1)
00603       std::cout<<"Repartitioning took "<<time.elapsed()<<" seconds."<<std::endl;
00604 
00605     ri.setSetup();
00606 
00607 #ifdef DEBUG_REPART
00608     ri.checkInterface(origComm.indexSet(), newComm->indexSet(), origComm.communicator());
00609 #endif
00610 
00611     redistributeMatrix(const_cast<M&>(origMatrix), newMatrix, origComm, *newComm, ri);
00612 
00613 #ifdef DEBUG_REPART
00614     if(origComm.communicator().rank()==0)
00615       std::cout<<"Original matrix"<<std::endl;
00616     origComm.communicator().barrier();
00617     if(newComm->communicator().size()>0)
00618       printGlobalSparseMatrix(newMatrix, *newComm, std::cout);
00619     origComm.communicator().barrier();
00620 #endif
00621 
00622     if(origComm.communicator().rank()==0  && criterion.debugLevel()>1)
00623       std::cout<<"Redistributing matrix took "<<time.elapsed()<<" seconds."<<std::endl;
00624     return existentOnRedist;
00625     
00626   }
00627   
00628   template<typename M>
00629   bool repartitionAndDistributeMatrix(M& origMatrix, M& newMatrix, 
00630                                       SequentialInformation& origComm, 
00631                                       SequentialInformation& newComm, 
00632                                       RedistributeInformation<SequentialInformation>& ri)
00633   {
00634     return true;
00635   }
00636 
00637     template<class M, class IS, class A>
00638     MatrixHierarchy<M,IS,A>::MatrixHierarchy(const MatrixOperator& fineOperator,
00639                                              const ParallelInformation& pinfo)
00640       : matrices_(const_cast<MatrixOperator&>(fineOperator)),
00641         parallelInformation_(const_cast<ParallelInformation&>(pinfo))
00642     {
00643       dune_static_assert((static_cast<int>(MatrixOperator::category) == 
00644                           static_cast<int>(SolverCategory::sequential) ||
00645                           static_cast<int>(MatrixOperator::category) == 
00646                           static_cast<int>(SolverCategory::overlapping) ||
00647                           static_cast<int>(MatrixOperator::category) == 
00648                           static_cast<int>(SolverCategory::nonoverlapping)),
00649                          "MatrixOperator must be of category sequential or overlapping or nonoverlapping");
00650       if (static_cast<int>(MatrixOperator::category) != static_cast<int>(pinfo.getSolverCategory()))
00651         DUNE_THROW(ISTLError, "MatrixOperator and ParallelInformation must belong to the same category!");
00652       
00653     }
00654     
00655     template<class M, class IS, class A>
00656     template<typename O, typename T>
00657     void MatrixHierarchy<M,IS,A>::build(const T& criterion)
00658     {
00659       prolongDamp_ = criterion.getProlongationDampingFactor();
00660       typedef O OverlapFlags;
00661       typedef typename ParallelMatrixHierarchy::Iterator MatIterator;
00662       typedef typename ParallelInformationHierarchy::Iterator PInfoIterator;
00663 
00664       static const int noints=(Dune::Amg::MAX_PROCESSES/4096>0)?(Dune::Amg::MAX_PROCESSES/4096):1;
00665       
00666       typedef bigunsignedint<sizeof(int)*8*noints> BIGINT;      
00667       GalerkinProduct<ParallelInformation> productBuilder;
00668       MatIterator mlevel = matrices_.finest();
00669       MatrixStats<typename M::matrix_type,MINIMAL_DEBUG_LEVEL<=INFO_DEBUG_LEVEL>::stats(mlevel->getmat());
00670       
00671       PInfoIterator infoLevel = parallelInformation_.finest();
00672       BIGINT finenonzeros=countNonZeros(mlevel->getmat());
00673       finenonzeros = infoLevel->communicator().sum(finenonzeros);
00674       BIGINT allnonzeros = finenonzeros;
00675 
00676       
00677       int level = 0;
00678       int rank = 0;
00679 
00680       BIGINT unknowns = mlevel->getmat().N();
00681 
00682       unknowns = infoLevel->communicator().sum(unknowns);
00683       double dunknowns=unknowns.todouble();
00684       infoLevel->buildGlobalLookup(mlevel->getmat().N());
00685       redistributes_.push_back(RedistributeInfoType());
00686 
00687       for(; level < criterion.maxLevel(); ++level, ++mlevel){
00688         assert(matrices_.levels()==redistributes_.size());
00689         rank = infoLevel->communicator().rank();
00690         if(rank==0 && criterion.debugLevel()>1)
00691           std::cout<<"Level "<<level<<" has "<<dunknowns<<" unknowns, "<<dunknowns/infoLevel->communicator().size()
00692                    <<" unknowns per proc (procs="<<infoLevel->communicator().size()<<")"<<std::endl;
00693 
00694         MatrixOperator* matrix=&(*mlevel);
00695         ParallelInformation* info =&(*infoLevel);
00696 
00697         if((
00698 #if HAVE_PARMETIS
00699               criterion.accumulate()==successiveAccu
00700 #else
00701               false
00702 #endif
00703              || (criterion.accumulate()==atOnceAccu
00704                  && dunknowns < 30*infoLevel->communicator().size()))
00705            && infoLevel->communicator().size()>1 && 
00706            dunknowns/infoLevel->communicator().size() <= criterion.coarsenTarget())
00707           {
00708             // accumulate to fewer processors
00709             Matrix* redistMat= new Matrix();
00710             ParallelInformation* redistComm=0;
00711             std::size_t nodomains = (std::size_t)std::ceil(dunknowns/(criterion.minAggregateSize()
00712                                                                   *criterion.coarsenTarget()));
00713             if( nodomains<=criterion.minAggregateSize()/2 || 
00714                 dunknowns <= criterion.coarsenTarget() )
00715               nodomains=1;
00716 
00717             bool existentOnNextLevel = 
00718               repartitionAndDistributeMatrix(mlevel->getmat(), *redistMat, *infoLevel,
00719                                              redistComm, redistributes_.back(), nodomains,
00720                                              criterion);
00721             BIGINT unknowns = redistMat->N();
00722             unknowns = infoLevel->communicator().sum(unknowns);
00723             dunknowns= unknowns.todouble();
00724             if(redistComm->communicator().rank()==0 && criterion.debugLevel()>1)
00725               std::cout<<"Level "<<level<<" (redistributed) has "<<dunknowns<<" unknowns, "<<dunknowns/redistComm->communicator().size()
00726                        <<" unknowns per proc (procs="<<redistComm->communicator().size()<<")"<<std::endl;
00727             MatrixArgs args(*redistMat, *redistComm);
00728             mlevel.addRedistributed(ConstructionTraits<MatrixOperator>::construct(args));
00729             assert(mlevel.isRedistributed());
00730             infoLevel.addRedistributed(redistComm);
00731             infoLevel->freeGlobalLookup();
00732             
00733             if(!existentOnNextLevel)
00734               // We do not hold any data on the redistributed partitioning
00735               break;
00736             
00737             // Work on the redistributed Matrix from now on
00738             matrix = &(mlevel.getRedistributed());
00739             info = &(infoLevel.getRedistributed());
00740             info->buildGlobalLookup(matrix->getmat().N());
00741           }
00742         
00743         rank = info->communicator().rank();
00744         if(dunknowns <= criterion.coarsenTarget())
00745           // No further coarsening needed
00746            break;
00747 
00748         typedef PropertiesGraphCreator<MatrixOperator> GraphCreator;
00749         typedef typename GraphCreator::PropertiesGraph PropertiesGraph;
00750         typedef typename GraphCreator::MatrixGraph MatrixGraph;
00751         typedef typename GraphCreator::GraphTuple  GraphTuple;
00752 
00753         typedef typename PropertiesGraph::VertexDescriptor Vertex;
00754         
00755         std::vector<bool> excluded(matrix->getmat().N(), false);
00756 
00757         GraphTuple graphs = GraphCreator::create(*matrix, excluded, *info, OverlapFlags());
00758         
00759         AggregatesMap* aggregatesMap=new AggregatesMap(get<1>(graphs)->maxVertex()+1);
00760 
00761         aggregatesMaps_.push_back(aggregatesMap);
00762 
00763         Timer watch;
00764         watch.reset();
00765         int noAggregates, isoAggregates, oneAggregates, skippedAggregates;
00766         
00767         tie(noAggregates, isoAggregates, oneAggregates, skippedAggregates) =    
00768           aggregatesMap->buildAggregates(matrix->getmat(), *(get<1>(graphs)), criterion, level==0);
00769 
00770         if(rank==0 && criterion.debugLevel()>2)
00771           std::cout<<" Have built "<<noAggregates<<" aggregates totally ("<<isoAggregates<<" isolated aggregates, "<<
00772             oneAggregates<<" aggregates of one vertex,  and skipped "<<
00773             skippedAggregates<<" aggregates)."<<std::endl;
00774 #ifdef TEST_AGGLO
00775         {
00776         // calculate size of local matrix in the distributed direction
00777         int start, end, overlapStart, overlapEnd;
00778         int procs=info->communicator().rank();
00779         int n = UNKNOWNS/procs; // number of unknowns per process
00780         int bigger = UNKNOWNS%procs; // number of process with n+1 unknows
00781 
00782         // Compute owner region
00783         if(rank<bigger){
00784           start = rank*(n+1);
00785           end   = (rank+1)*(n+1);
00786         }else{
00787           start = bigger + rank * n;
00788           end   = bigger + (rank + 1) * n;
00789         }
00790   
00791         // Compute overlap region
00792         if(start>0)
00793           overlapStart = start - 1;
00794         else
00795           overlapStart = start;
00796         
00797         if(end<UNKNOWNS)
00798           overlapEnd = end + 1;
00799         else
00800           overlapEnd = end;
00801   
00802         assert((UNKNOWNS)*(overlapEnd-overlapStart)==aggregatesMap->noVertices());
00803         for(int j=0; j< UNKNOWNS; ++j)
00804           for(int i=0; i < UNKNOWNS; ++i)
00805             {
00806               if(i>=overlapStart && i<overlapEnd)
00807                 {
00808                   int no = (j/2)*((UNKNOWNS)/2)+i/2;
00809                   (*aggregatesMap)[j*(overlapEnd-overlapStart)+i-overlapStart]=no;
00810                 }
00811             }
00812         }
00813 #endif
00814         if(criterion.debugLevel()>1 && info->communicator().rank()==0)
00815           std::cout<<"aggregating finished."<<std::endl;
00816 
00817         BIGINT gnoAggregates=noAggregates;
00818         gnoAggregates = info->communicator().sum(gnoAggregates);
00819         double dgnoAggregates = gnoAggregates.todouble();
00820 #ifdef TEST_AGGLO
00821         BIGINT gnoAggregates=((UNKNOWNS)/2)*((UNKNOWNS)/2);
00822 #endif
00823 
00824         if(criterion.debugLevel()>2 && rank==0)
00825           std::cout << "Building "<<dgnoAggregates<<" aggregates took "<<watch.elapsed()<<" seconds."<<std::endl;       
00826 
00827         if(dgnoAggregates==0 || dunknowns/dgnoAggregates<criterion.minCoarsenRate())
00828     {
00829             if(rank==0)
00830         {
00831               if(dgnoAggregates>0)
00832                 std::cerr << "Stopped coarsening because of rate breakdown "<<dunknowns<<"/"<<dgnoAggregates
00833                           <<"="<<dunknowns/dgnoAggregates<<"<"
00834                   <<criterion.minCoarsenRate()<<std::endl;
00835               else
00836             std::cerr<< "Could not build any aggregates. Probably no connected nodes."<<std::endl;
00837         }
00838             aggregatesMap->free();
00839             delete aggregatesMap;
00840             aggregatesMaps_.pop_back();
00841 
00842             if(criterion.accumulate() && mlevel.isRedistributed() && info->communicator().size()>1){
00843               // coarse level matrix was already redistributed, but to more than 1 process
00844               // Therefore need to delete the redistribution. Further down it will
00845               // then be redistributed to 1 process
00846               delete &(mlevel.getRedistributed().getmat());
00847               mlevel.deleteRedistributed();
00848               delete &(infoLevel.getRedistributed());
00849               infoLevel.deleteRedistributed();
00850               redistributes_.back().resetSetup();
00851             }
00852                 
00853             break;
00854     }
00855         unknowns =  noAggregates;
00856     dunknowns = dgnoAggregates;
00857 
00858     CommunicationArgs commargs(info->communicator(),info->getSolverCategory());
00859     parallelInformation_.addCoarser(commargs);
00860 
00861         ++infoLevel; // parallel information on coarse level
00862                 
00863         typename PropertyMapTypeSelector<VertexVisitedTag,PropertiesGraph>::Type visitedMap = 
00864           get(VertexVisitedTag(), *(get<1>(graphs)));
00865   
00866         watch.reset();
00867         int aggregates = IndicesCoarsener<ParallelInformation,OverlapFlags>
00868 	  ::coarsen(*info,
00869                     *(get<1>(graphs)),
00870                     visitedMap,
00871                     *aggregatesMap,
00872                     *infoLevel);
00873         
00874         GraphCreator::free(graphs);
00875         
00876         if(criterion.debugLevel()>2){
00877           if(rank==0)
00878             std::cout<<"Coarsening of index sets took "<<watch.elapsed()<<" seconds."<<std::endl;
00879         }
00880         
00881         watch.reset();
00882 
00883         infoLevel->buildGlobalLookup(aggregates);
00884         AggregatesPublisher<Vertex,OverlapFlags,ParallelInformation>::publish(*aggregatesMap,
00885                                                                               *info,
00886                                                                               infoLevel->globalLookup());
00887         
00888                 
00889         if(criterion.debugLevel()>2){
00890           if(rank==0)
00891             std::cout<<"Communicating global aggregate numbers took "<<watch.elapsed()<<" seconds."<<std::endl;
00892         }
00893 
00894         watch.reset();
00895         std::vector<bool>& visited=excluded;
00896         
00897         typedef std::vector<bool>::iterator Iterator;
00898         typedef IteratorPropertyMap<Iterator, IdentityMap> VisitedMap2;
00899         Iterator end = visited.end();
00900         for(Iterator iter= visited.begin(); iter != end; ++iter)
00901           *iter=false;
00902         
00903         VisitedMap2 visitedMap2(visited.begin(), Dune::IdentityMap());
00904 
00905         typename MatrixOperator::matrix_type* coarseMatrix;
00906 
00907         coarseMatrix = productBuilder.build(matrix->getmat(), *(get<0>(graphs)), visitedMap2, 
00908                                             *info, 
00909                                             *aggregatesMap,
00910                                             aggregates,
00911                                             OverlapFlags());
00912         
00913         info->freeGlobalLookup();
00914         
00915         delete get<0>(graphs);
00916         productBuilder.calculate(matrix->getmat(), *aggregatesMap, *coarseMatrix, *infoLevel, OverlapFlags());
00917         
00918         if(criterion.debugLevel()>2){
00919           if(rank==0)
00920             std::cout<<"Calculation of Galerkin product took "<<watch.elapsed()<<" seconds."<<std::endl;
00921         }
00922         
00923         BIGINT nonzeros = countNonZeros(*coarseMatrix);
00924         allnonzeros = allnonzeros + infoLevel->communicator().sum(nonzeros);
00925         MatrixArgs args(*coarseMatrix, *infoLevel);
00926         
00927         matrices_.addCoarser(args);
00928         redistributes_.push_back(RedistributeInfoType());
00929       } // end level loop
00930       
00931 
00932       infoLevel->freeGlobalLookup();
00933       
00934       built_=true;
00935       AggregatesMap* aggregatesMap=new AggregatesMap(0);
00936       aggregatesMaps_.push_back(aggregatesMap);
00937 
00938       if(criterion.debugLevel()>0){
00939         if(level==criterion.maxLevel()){
00940           BIGINT unknowns = mlevel->getmat().N();
00941           unknowns = infoLevel->communicator().sum(unknowns);
00942           double dunknowns = unknowns.todouble();
00943           if(rank==0 && criterion.debugLevel()>1){
00944             std::cout<<"Level "<<level<<" has "<<dunknowns<<" unknowns, "<<dunknowns/infoLevel->communicator().size()
00945                      <<" unknowns per proc (procs="<<infoLevel->communicator().size()<<")"<<std::endl;
00946           }
00947         }
00948       }
00949 
00950       if(criterion.accumulate() && !redistributes_.back().isSetup() && 
00951          infoLevel->communicator().size()>1){ 
00952 #if HAVE_MPI && !HAVE_PARMETIS
00953         if(criterion.accumulate()==successiveAccu &&
00954            infoLevel->communicator().rank()==0)
00955           std::cerr<<"Successive accumulation of data on coarse levels only works with ParMETIS installed."
00956                    <<"  Fell back to accumulation to one domain on coarsest level"<<std::endl;
00957 #endif
00958           
00959         // accumulate to fewer processors
00960         Matrix* redistMat= new Matrix();
00961         ParallelInformation* redistComm=0;
00962         int nodomains = 1;
00963             
00964         repartitionAndDistributeMatrix(mlevel->getmat(), *redistMat, *infoLevel,
00965                                        redistComm, redistributes_.back(), nodomains,criterion);
00966         MatrixArgs args(*redistMat, *redistComm);
00967         BIGINT unknowns = redistMat->N();
00968         unknowns = infoLevel->communicator().sum(unknowns);
00969         
00970         if(redistComm->communicator().rank()==0 && criterion.debugLevel()>1){
00971           double dunknowns= unknowns.todouble();
00972           std::cout<<"Level "<<level<<" redistributed has "<<dunknowns<<" unknowns, "<<dunknowns/redistComm->communicator().size()
00973                      <<" unknowns per proc (procs="<<redistComm->communicator().size()<<")"<<std::endl;
00974         }
00975         mlevel.addRedistributed(ConstructionTraits<MatrixOperator>::construct(args));
00976         infoLevel.addRedistributed(redistComm);
00977         infoLevel->freeGlobalLookup();
00978       }
00979 
00980         int levels = matrices_.levels();
00981         maxlevels_ = parallelInformation_.finest()->communicator().max(levels);
00982         assert(matrices_.levels()==redistributes_.size());
00983         if(hasCoarsest() && rank==0 && criterion.debugLevel()>1)
00984           std::cout<<"operator complexity: "<<allnonzeros.todouble()/finenonzeros.todouble()<<std::endl;
00985 
00986     }
00987     
00988     template<class M, class IS, class A>
00989     const typename MatrixHierarchy<M,IS,A>::ParallelMatrixHierarchy& 
00990     MatrixHierarchy<M,IS,A>::matrices() const
00991     {
00992       return matrices_;
00993     }
00994     
00995     template<class M, class IS, class A>
00996     const typename MatrixHierarchy<M,IS,A>::ParallelInformationHierarchy& 
00997     MatrixHierarchy<M,IS,A>::parallelInformation() const
00998     {
00999       return parallelInformation_;
01000     }
01001     
01002     template<class M, class IS, class A>
01003     void MatrixHierarchy<M,IS,A>::getCoarsestAggregatesOnFinest(std::vector<std::size_t>& data) const
01004     {
01005       int levels=aggregatesMaps().size();
01006       int maxlevels=parallelInformation_.finest()->communicator().max(levels);
01007       std::size_t size=(*(aggregatesMaps().begin()))->noVertices();
01008       // We need an auxiliary vector for the consecutive prolongation.
01009       std::vector<std::size_t> tmp;
01010       std::vector<std::size_t> *coarse, *fine;
01011 
01012       // make sure the allocated space suffices.
01013       tmp.reserve(size);
01014       data.reserve(size);
01015 
01016       // Correctly assign coarse and fine for the first prolongation such that
01017       // we end up in data in the end.
01018       if(levels%2==0){
01019         coarse=&tmp;
01020         fine=&data;
01021       }else{
01022         coarse=&data;
01023         fine=&tmp;
01024       }
01025 
01026       // Number the unknowns on the coarsest level consecutively for each process.
01027       if(levels==maxlevels){
01028         const AggregatesMap& map = *(*(++aggregatesMaps().rbegin()));
01029         std::size_t m=0;
01030         
01031         for(typename AggregatesMap::const_iterator iter = map.begin(); iter != map.end(); ++iter)
01032           if(*iter< AggregatesMap::ISOLATED)
01033             m=std::max(*iter,m);
01034 
01035         coarse->resize(m+1);
01036         std::size_t i=0;
01037         srand((unsigned)std::clock());
01038         std::set<size_t> used;
01039         for(typename std::vector<std::size_t>::iterator iter=coarse->begin(); iter != coarse->end();
01040             ++iter, ++i)
01041           {
01042             std::pair<std::set<std::size_t>::iterator,bool> ibpair
01043               = used.insert(static_cast<std::size_t>((((double)rand())/(RAND_MAX+1.0)))*coarse->size());
01044             
01045             while(!ibpair.second)
01046               ibpair = used.insert(static_cast<std::size_t>((((double)rand())/(RAND_MAX+1.0))*coarse->size()));
01047             *iter=*(ibpair.first);
01048           }
01049       }
01050 
01051       typename ParallelInformationHierarchy::Iterator pinfo = parallelInformation().coarsest();
01052       --pinfo;
01053             
01054       // Now consecutively project the numbers to the finest level.
01055       for(typename AggregatesMapList::const_reverse_iterator aggregates=++aggregatesMaps().rbegin();
01056           aggregates != aggregatesMaps().rend(); ++aggregates,--levels){
01057         
01058         fine->resize((*aggregates)->noVertices());
01059         fine->assign(fine->size(), 0);
01060         Transfer<typename AggregatesMap::AggregateDescriptor, std::vector<std::size_t>, ParallelInformation>
01061           ::prolongate(*(*aggregates), *coarse, *fine, static_cast<std::size_t>(1), *pinfo);
01062         --pinfo;
01063         std::swap(coarse, fine);
01064       }
01065 
01066       // Assertion to check that we really projected to data on the last step.
01067       assert(coarse==&data);
01068     }
01069     
01070     template<class M, class IS, class A>
01071     const typename MatrixHierarchy<M,IS,A>::AggregatesMapList& 
01072     MatrixHierarchy<M,IS,A>::aggregatesMaps() const
01073     {
01074       return aggregatesMaps_;
01075     }
01076     template<class M, class IS, class A>
01077     const typename MatrixHierarchy<M,IS,A>::RedistributeInfoList& 
01078     MatrixHierarchy<M,IS,A>::redistributeInformation() const 
01079     {
01080       return redistributes_;
01081     }
01082     
01083     template<class M, class IS, class A>
01084     MatrixHierarchy<M,IS,A>::~MatrixHierarchy()
01085     {
01086       typedef typename AggregatesMapList::reverse_iterator AggregatesMapIterator;
01087       typedef typename ParallelMatrixHierarchy::Iterator Iterator;
01088       typedef typename ParallelInformationHierarchy::Iterator InfoIterator;
01089       
01090       AggregatesMapIterator amap = aggregatesMaps_.rbegin();
01091       InfoIterator info = parallelInformation_.coarsest();
01092       for(Iterator level=matrices_.coarsest(), finest=matrices_.finest(); level != finest;  --level, --info, ++amap){
01093         (*amap)->free();
01094         delete *amap;
01095         delete &level->getmat();
01096         if(level.isRedistributed())
01097           delete &(level.getRedistributed().getmat());
01098       }
01099       delete *amap;
01100     }
01101 
01102     template<class M, class IS, class A>
01103     template<class V, class TA>
01104     void MatrixHierarchy<M,IS,A>::coarsenVector(Hierarchy<BlockVector<V,TA> >& hierarchy) const
01105     {
01106       assert(hierarchy.levels()==1);
01107       typedef typename ParallelMatrixHierarchy::ConstIterator Iterator;
01108       typedef typename RedistributeInfoList::const_iterator RIter;
01109       RIter redist = redistributes_.begin();
01110 
01111       Iterator matrix = matrices_.finest(), coarsest = matrices_.coarsest();
01112       int level=0;
01113       if(redist->isSetup())
01114           hierarchy.addRedistributedOnCoarsest(matrix.getRedistributed().getmat().N());
01115       Dune::dvverb<<"Level "<<level<<" has "<<matrices_.finest()->getmat().N()<<" unknowns!"<<std::endl;
01116 
01117       while(matrix != coarsest){
01118         ++matrix; ++level; ++redist;
01119         Dune::dvverb<<"Level "<<level<<" has "<<matrix->getmat().N()<<" unknowns!"<<std::endl;
01120           
01121         hierarchy.addCoarser(matrix->getmat().N());
01122         if(redist->isSetup())
01123           hierarchy.addRedistributedOnCoarsest(matrix.getRedistributed().getmat().N());
01124         
01125       }
01126           
01127     }
01128     
01129     template<class M, class IS, class A>
01130     template<class S, class TA>
01131     void MatrixHierarchy<M,IS,A>::coarsenSmoother(Hierarchy<S,TA>& smoothers, 
01132                                                     const typename SmootherTraits<S>::Arguments& sargs) const
01133     {
01134       assert(smoothers.levels()==0);
01135       typedef typename ParallelMatrixHierarchy::ConstIterator MatrixIterator;
01136       typedef typename ParallelInformationHierarchy::ConstIterator PinfoIterator;
01137       typedef typename AggregatesMapList::const_iterator AggregatesIterator;
01138       
01139       typename ConstructionTraits<S>::Arguments cargs;
01140       cargs.setArgs(sargs);
01141       PinfoIterator pinfo = parallelInformation_.finest();
01142       AggregatesIterator aggregates = aggregatesMaps_.begin();
01143       int level=0;
01144       for(MatrixIterator matrix = matrices_.finest(), coarsest = matrices_.coarsest(); 
01145           matrix != coarsest; ++matrix, ++pinfo, ++aggregates, ++level){
01146         cargs.setMatrix(matrix->getmat(), **aggregates);
01147         cargs.setComm(*pinfo);
01148         smoothers.addCoarser(cargs);
01149       }
01150       if(maxlevels()>levels()){
01151         // This is not the globally coarsest level and therefore smoothing is needed
01152         cargs.setMatrix(matrices_.coarsest()->getmat(), **aggregates);
01153         cargs.setComm(*pinfo);
01154         smoothers.addCoarser(cargs);
01155         ++level;
01156       }
01157     }
01158     
01159     template<class M, class IS, class A>
01160     template<class F>
01161     void MatrixHierarchy<M,IS,A>::recalculateGalerkin(const F& copyFlags)
01162     {
01163       typedef typename AggregatesMapList::iterator AggregatesMapIterator;
01164       typedef typename ParallelMatrixHierarchy::Iterator Iterator;
01165       typedef typename ParallelInformationHierarchy::Iterator InfoIterator;
01166       
01167       AggregatesMapIterator amap = aggregatesMaps_.begin();
01168       BaseGalerkinProduct productBuilder;
01169       InfoIterator info = parallelInformation_.finest();
01170       typename RedistributeInfoList::iterator riIter = redistributes_.begin();
01171       Iterator level = matrices_.finest(), coarsest=matrices_.coarsest();
01172       if(level.isRedistributed()){
01173         info->buildGlobalLookup(info->indexSet().size());
01174         redistributeMatrixEntries(const_cast<Matrix&>(level->getmat()), 
01175                                   const_cast<Matrix&>(level.getRedistributed().getmat()),
01176                                   *info,info.getRedistributed(), *riIter);
01177         info->freeGlobalLookup();
01178       }
01179       
01180       for(; level!=coarsest; ++amap){
01181         const Matrix& fine = (level.isRedistributed()?level.getRedistributed():*level).getmat();
01182         ++level;
01183         ++info;
01184         ++riIter;
01185         productBuilder.calculate(fine, *(*amap), const_cast<Matrix&>(level->getmat()), *info, copyFlags);
01186         if(level.isRedistributed()){
01187         info->buildGlobalLookup(info->indexSet().size());
01188           redistributeMatrixEntries(const_cast<Matrix&>(level->getmat()), 
01189                                     const_cast<Matrix&>(level.getRedistributed().getmat()), *info,
01190                                     info.getRedistributed(), *riIter);
01191           info->freeGlobalLookup();
01192         }
01193       }
01194     }
01195 
01196     template<class M, class IS, class A>
01197     std::size_t MatrixHierarchy<M,IS,A>::levels() const
01198     {
01199       return matrices_.levels();
01200     }
01201 
01202     template<class M, class IS, class A>
01203     std::size_t MatrixHierarchy<M,IS,A>::maxlevels() const
01204     {
01205       return maxlevels_;
01206     }
01207 
01208     template<class M, class IS, class A>
01209     bool MatrixHierarchy<M,IS,A>::hasCoarsest() const
01210     {
01211       return levels()==maxlevels() && 
01212         (!matrices_.coarsest().isRedistributed() ||matrices_.coarsest()->getmat().N()>0);
01213     }
01214     
01215     template<class M, class IS, class A>
01216     bool MatrixHierarchy<M,IS,A>::isBuilt() const
01217     {
01218       return built_;
01219     }
01220     
01221     template<class T, class A>
01222     Hierarchy<T,A>::Hierarchy()
01223       : finest_(0), coarsest_(0), nonAllocated_(0), allocator_(), levels_(0)
01224     {}
01225 
01226     template<class T, class A>
01227     Hierarchy<T,A>::Hierarchy(MemberType& first)
01228       : allocator_()
01229     {
01230       finest_ = allocator_.allocate(1,0);
01231       finest_->element_ = &first;
01232       finest_->redistributed_ = 0;
01233       nonAllocated_ = finest_;
01234       coarsest_ = finest_;
01235       coarsest_->coarser_ = coarsest_->finer_ = 0;
01236       levels_ = 1;
01237     }
01238 
01239     template<class T, class A>
01240     Hierarchy<T,A>::~Hierarchy()
01241     {
01242       while(coarsest_){
01243         Element* current = coarsest_;
01244         coarsest_ = coarsest_->finer_;
01245         if(current != nonAllocated_){
01246           if(current->redistributed_)
01247             ConstructionTraits<T>::deconstruct(current->redistributed_);
01248           ConstructionTraits<T>::deconstruct(current->element_);
01249         }
01250         allocator_.deallocate(current, 1);
01251         //coarsest_->coarser_ = 0;
01252       }
01253     }
01254 
01255     template<class T, class A>
01256     std::size_t Hierarchy<T,A>::levels() const
01257     {
01258       return levels_;
01259     }
01260 
01261     template<class T, class A>
01262     void Hierarchy<T,A>::addRedistributedOnCoarsest(Arguments& args)
01263     {
01264       coarsest_->redistributed_ = ConstructionTraits<MemberType>::construct(args);
01265     }
01266     
01267     template<class T, class A>
01268     void Hierarchy<T,A>::addCoarser(Arguments& args)
01269     {
01270       if(!coarsest_){
01271         assert(!finest_);
01272         coarsest_ = allocator_.allocate(1,0);
01273         coarsest_->element_ = ConstructionTraits<MemberType>::construct(args);
01274         finest_ = coarsest_;
01275         coarsest_->finer_ = 0;
01276       }else{
01277         coarsest_->coarser_ = allocator_.allocate(1,0);
01278         coarsest_->coarser_->finer_ = coarsest_;
01279         coarsest_ = coarsest_->coarser_;
01280         coarsest_->element_ = ConstructionTraits<MemberType>::construct(args);
01281       }
01282       coarsest_->redistributed_ = 0;
01283       coarsest_->coarser_=0;
01284       ++levels_;
01285     }
01286    
01287     
01288     template<class T, class A>
01289     void Hierarchy<T,A>::addFiner(Arguments& args)
01290     {
01291       if(!finest_){
01292         assert(!coarsest_);
01293         finest_ = allocator_.allocate(1,0);
01294         finest_->element = ConstructionTraits<T>::construct(args);
01295         coarsest_ = finest_;
01296         coarsest_->coarser_ = coarsest_->finer_ = 0;
01297       }else{
01298         finest_->finer_ = allocator_.allocate(1,0);
01299         finest_->finer_->coarser_ = finest_;
01300         finest_ = finest_->finer_;
01301         finest_->finer = 0;
01302         finest_->element = ConstructionTraits<T>::construct(args);
01303       }
01304       ++levels_;
01305     }
01306 
01307     template<class T, class A>
01308     typename Hierarchy<T,A>::Iterator Hierarchy<T,A>::finest()
01309     {
01310       return Iterator(finest_);
01311     }
01312 
01313     template<class T, class A>
01314     typename Hierarchy<T,A>::Iterator Hierarchy<T,A>::coarsest()
01315     {
01316       return Iterator(coarsest_);
01317     }
01318 
01319     template<class T, class A>
01320     typename Hierarchy<T,A>::ConstIterator Hierarchy<T,A>::finest() const
01321     {
01322       return ConstIterator(finest_);
01323     }
01324 
01325     template<class T, class A>
01326     typename Hierarchy<T,A>::ConstIterator Hierarchy<T,A>::coarsest() const
01327     {
01328       return ConstIterator(coarsest_);
01329     }
01331   }// namespace Amg
01332 }// namespace Dune
01333 
01334 #endif