dune-istl
2.2.0
|
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