4 #ifndef DUNE_AMGHIERARCHY_HH
5 #define DUNE_AMGHIERARCHY_HH
16 #include<dune/common/stdstreams.hh>
17 #include<dune/common/timer.hh>
18 #include<dune/common/tuples.hh>
19 #include<dune/common/bigunsignedint.hh>
21 #include<dune/common/parallel/indexset.hh>
64 template<
typename T,
typename A=std::allocator<T> >
73 template<
typename T1,
typename T2>
109 typedef typename A::template rebind<Element>::other
Allocator;
144 template<
class C,
class T1>
146 :
public BidirectionalIteratorFacade<LevelIterator<C,T1>,T1,T1&>
150 friend class
LevelIterator<const typename remove_const<C>::type,
151 const typename remove_const<T1>::type >;
166 : element_(other.element_)
172 : element_(other.element_)
181 return element_ == other.element_;
190 return element_ == other.element_;
196 return *(element_->element_);
202 element_ = element_->coarser_;
208 element_ = element_->finer_;
217 return element_->redistributed_;
226 assert(element_->redistributed_);
227 return *element_->redistributed_;
231 element_->redistributed_ = t;
236 element_->redistributed_ = 0;
278 std::size_t
levels()
const;
289 Element* nonAllocated_;
302 template<
class M,
class PI,
class A=std::allocator<M> >
310 typedef typename MatrixOperator::matrix_type
Matrix;
328 typedef typename Allocator::template rebind<AggregatesMap*>::other
AAllocator;
337 typedef typename Allocator::template rebind<RedistributeInfoType>::other
RILAllocator;
358 template<
typename O,
typename T>
359 void build(
const T& criterion);
375 template<
class V,
class TA>
383 template<
class S,
class TA>
391 std::size_t
levels()
const;
468 typename MatrixOperator::field_type prolongDamp_;
473 template<
class Matrix,
bool pr
int>
480 static void stats(
const Matrix& matrix)
484 template<
class Matrix>
485 struct MatrixStats<
Matrix,true>
494 min=std::numeric_limits<size_type>::max();
501 min=std::min(min, row.size());
502 max=std::max(max, row.size());
513 static void stats(
const Matrix& matrix)
515 calc c= for_each(matrix.begin(), matrix.end(), calc());
516 dinfo<<
"Matrix row: min="<<c.min<<
" max="<<c.max
517 <<
" average="<<
static_cast<double>(c.sum)/matrix.N()
557 template<
typename M,
typename C1>
562 int nparts, C1& criterion)
564 DUNE_THROW(NotImplemented,
"Redistribution does not make sense in sequential code!");
568 template<
typename M,
typename C,
typename C1>
571 int nparts, C1& criterion)
574 #ifdef AMG_REPART_ON_COMM_GRAPH
578 criterion.debugLevel()>1);
587 MatrixGraph graph(origMatrix);
588 PropertiesGraph pgraph(graph);
592 if(origComm.communicator().rank()==0)
593 std::cout<<
"Original matrix"<<std::endl;
594 origComm.communicator().barrier();
598 newComm, ri.getInterface(),
599 criterion.debugLevel()>1);
600 #endif // if else AMG_REPART
602 if(origComm.communicator().rank()==0 && criterion.debugLevel()>1)
603 std::cout<<
"Repartitioning took "<<time.elapsed()<<
" seconds."<<std::endl;
608 ri.checkInterface(origComm.indexSet(), newComm->indexSet(), origComm.communicator());
614 if(origComm.communicator().rank()==0)
615 std::cout<<
"Original matrix"<<std::endl;
616 origComm.communicator().barrier();
617 if(newComm->communicator().size()>0)
619 origComm.communicator().barrier();
622 if(origComm.communicator().rank()==0 && criterion.debugLevel()>1)
623 std::cout<<
"Redistributing matrix took "<<time.elapsed()<<
" seconds."<<std::endl;
624 return existentOnRedist;
637 template<
class M,
class IS,
class A>
643 dune_static_assert((static_cast<int>(MatrixOperator::category) ==
645 static_cast<int>(MatrixOperator::category) ==
647 static_cast<int>(MatrixOperator::category) ==
649 "MatrixOperator must be of category sequential or overlapping or nonoverlapping");
650 if (static_cast<int>(MatrixOperator::category) != static_cast<int>(pinfo.getSolverCategory()))
651 DUNE_THROW(
ISTLError,
"MatrixOperator and ParallelInformation must belong to the same category!");
655 template<
class M,
class IS,
class A>
656 template<
typename O,
typename T>
659 prolongDamp_ = criterion.getProlongationDampingFactor();
660 typedef O OverlapFlags;
666 typedef bigunsignedint<sizeof(int)*8*noints> BIGINT;
668 MatIterator mlevel = matrices_.finest();
669 MatrixStats<typename M::matrix_type,MINIMAL_DEBUG_LEVEL<=INFO_DEBUG_LEVEL>::stats(mlevel->getmat());
671 PInfoIterator infoLevel = parallelInformation_.finest();
673 finenonzeros = infoLevel->communicator().sum(finenonzeros);
674 BIGINT allnonzeros = finenonzeros;
680 BIGINT unknowns = mlevel->getmat().N();
682 unknowns = infoLevel->communicator().sum(unknowns);
683 double dunknowns=unknowns.todouble();
684 infoLevel->buildGlobalLookup(mlevel->getmat().N());
687 for(; level < criterion.maxLevel(); ++level, ++mlevel){
688 assert(matrices_.levels()==redistributes_.size());
689 rank = infoLevel->communicator().rank();
690 if(rank==0 && criterion.debugLevel()>1)
691 std::cout<<
"Level "<<level<<
" has "<<dunknowns<<
" unknowns, "<<dunknowns/infoLevel->communicator().size()
692 <<
" unknowns per proc (procs="<<infoLevel->communicator().size()<<
")"<<std::endl;
704 && dunknowns < 30*infoLevel->communicator().size()))
705 && infoLevel->communicator().size()>1 &&
706 dunknowns/infoLevel->communicator().size() <= criterion.coarsenTarget())
711 std::size_t nodomains = (std::size_t)std::ceil(dunknowns/(criterion.minAggregateSize()
712 *criterion.coarsenTarget()));
713 if( nodomains<=criterion.minAggregateSize()/2 ||
714 dunknowns <= criterion.coarsenTarget() )
717 bool existentOnNextLevel =
719 redistComm, redistributes_.back(), nodomains,
721 BIGINT unknowns = redistMat->N();
722 unknowns = infoLevel->communicator().sum(unknowns);
723 dunknowns= unknowns.todouble();
724 if(redistComm->communicator().rank()==0 && criterion.debugLevel()>1)
725 std::cout<<
"Level "<<level<<
" (redistributed) has "<<dunknowns<<
" unknowns, "<<dunknowns/redistComm->communicator().size()
726 <<
" unknowns per proc (procs="<<redistComm->communicator().size()<<
")"<<std::endl;
727 MatrixArgs args(*redistMat, *redistComm);
729 assert(mlevel.isRedistributed());
730 infoLevel.addRedistributed(redistComm);
731 infoLevel->freeGlobalLookup();
733 if(!existentOnNextLevel)
738 matrix = &(mlevel.getRedistributed());
739 info = &(infoLevel.getRedistributed());
740 info->buildGlobalLookup(matrix->getmat().N());
743 rank = info->communicator().rank();
744 if(dunknowns <= criterion.coarsenTarget())
750 typedef typename GraphCreator::MatrixGraph
MatrixGraph;
751 typedef typename GraphCreator::GraphTuple GraphTuple;
755 std::vector<bool> excluded(matrix->getmat().N(),
false);
757 GraphTuple graphs = GraphCreator::create(*matrix, excluded, *info, OverlapFlags());
761 aggregatesMaps_.push_back(aggregatesMap);
765 int noAggregates, isoAggregates, oneAggregates, skippedAggregates;
767 tie(noAggregates, isoAggregates, oneAggregates, skippedAggregates) =
768 aggregatesMap->
buildAggregates(matrix->getmat(), *(get<1>(graphs)), criterion, level==0);
770 if(rank==0 && criterion.debugLevel()>2)
771 std::cout<<
" Have built "<<noAggregates<<
" aggregates totally ("<<isoAggregates<<
" isolated aggregates, "<<
772 oneAggregates<<
" aggregates of one vertex, and skipped "<<
773 skippedAggregates<<
" aggregates)."<<std::endl;
777 int start, end, overlapStart, overlapEnd;
778 int procs=info->communicator().rank();
779 int n = UNKNOWNS/procs;
780 int bigger = UNKNOWNS%procs;
785 end = (rank+1)*(n+1);
787 start = bigger + rank * n;
788 end = bigger + (rank + 1) * n;
793 overlapStart = start - 1;
795 overlapStart = start;
798 overlapEnd = end + 1;
802 assert((UNKNOWNS)*(overlapEnd-overlapStart)==aggregatesMap->
noVertices());
803 for(
int j=0; j< UNKNOWNS; ++j)
804 for(
int i=0; i < UNKNOWNS; ++i)
806 if(i>=overlapStart && i<overlapEnd)
808 int no = (j/2)*((UNKNOWNS)/2)+i/2;
809 (*aggregatesMap)[j*(overlapEnd-overlapStart)+i-overlapStart]=no;
814 if(criterion.debugLevel()>1 && info->communicator().rank()==0)
815 std::cout<<
"aggregating finished."<<std::endl;
817 BIGINT gnoAggregates=noAggregates;
818 gnoAggregates = info->communicator().sum(gnoAggregates);
819 double dgnoAggregates = gnoAggregates.todouble();
821 BIGINT gnoAggregates=((UNKNOWNS)/2)*((UNKNOWNS)/2);
824 if(criterion.debugLevel()>2 && rank==0)
825 std::cout <<
"Building "<<dgnoAggregates<<
" aggregates took "<<watch.elapsed()<<
" seconds."<<std::endl;
827 if(dgnoAggregates==0 || dunknowns/dgnoAggregates<criterion.minCoarsenRate())
832 std::cerr <<
"Stopped coarsening because of rate breakdown "<<dunknowns<<
"/"<<dgnoAggregates
833 <<
"="<<dunknowns/dgnoAggregates<<
"<"
834 <<criterion.minCoarsenRate()<<std::endl;
836 std::cerr<<
"Could not build any aggregates. Probably no connected nodes."<<std::endl;
838 aggregatesMap->
free();
839 delete aggregatesMap;
840 aggregatesMaps_.pop_back();
842 if(criterion.accumulate() && mlevel.isRedistributed() && info->communicator().size()>1){
846 delete &(mlevel.getRedistributed().getmat());
847 mlevel.deleteRedistributed();
848 delete &(infoLevel.getRedistributed());
849 infoLevel.deleteRedistributed();
850 redistributes_.back().resetSetup();
855 unknowns = noAggregates;
856 dunknowns = dgnoAggregates;
858 CommunicationArgs commargs(info->communicator(),info->getSolverCategory());
859 parallelInformation_.addCoarser(commargs);
863 typename PropertyMapTypeSelector<VertexVisitedTag,PropertiesGraph>::Type visitedMap =
874 GraphCreator::free(graphs);
876 if(criterion.debugLevel()>2){
878 std::cout<<
"Coarsening of index sets took "<<watch.elapsed()<<
" seconds."<<std::endl;
883 infoLevel->buildGlobalLookup(aggregates);
886 infoLevel->globalLookup());
889 if(criterion.debugLevel()>2){
891 std::cout<<
"Communicating global aggregate numbers took "<<watch.elapsed()<<
" seconds."<<std::endl;
895 std::vector<bool>& visited=excluded;
897 typedef std::vector<bool>::iterator Iterator;
898 typedef IteratorPropertyMap<Iterator, IdentityMap> VisitedMap2;
899 Iterator end = visited.end();
900 for(Iterator iter= visited.begin(); iter != end; ++iter)
903 VisitedMap2 visitedMap2(visited.begin(), Dune::IdentityMap());
905 typename MatrixOperator::matrix_type* coarseMatrix;
907 coarseMatrix = productBuilder.
build(matrix->getmat(), *(get<0>(graphs)), visitedMap2,
913 info->freeGlobalLookup();
915 delete get<0>(graphs);
916 productBuilder.
calculate(matrix->getmat(), *aggregatesMap, *coarseMatrix, *infoLevel, OverlapFlags());
918 if(criterion.debugLevel()>2){
920 std::cout<<
"Calculation of Galerkin product took "<<watch.elapsed()<<
" seconds."<<std::endl;
924 allnonzeros = allnonzeros + infoLevel->communicator().sum(nonzeros);
925 MatrixArgs args(*coarseMatrix, *infoLevel);
927 matrices_.addCoarser(args);
932 infoLevel->freeGlobalLookup();
936 aggregatesMaps_.push_back(aggregatesMap);
938 if(criterion.debugLevel()>0){
939 if(level==criterion.maxLevel()){
940 BIGINT unknowns = mlevel->getmat().N();
941 unknowns = infoLevel->communicator().sum(unknowns);
942 double dunknowns = unknowns.todouble();
943 if(rank==0 && criterion.debugLevel()>1){
944 std::cout<<
"Level "<<level<<
" has "<<dunknowns<<
" unknowns, "<<dunknowns/infoLevel->communicator().size()
945 <<
" unknowns per proc (procs="<<infoLevel->communicator().size()<<
")"<<std::endl;
950 if(criterion.accumulate() && !redistributes_.back().isSetup() &&
951 infoLevel->communicator().size()>1){
952 #if HAVE_MPI && !HAVE_PARMETIS
954 infoLevel->communicator().rank()==0)
955 std::cerr<<
"Successive accumulation of data on coarse levels only works with ParMETIS installed."
956 <<
" Fell back to accumulation to one domain on coarsest level"<<std::endl;
965 redistComm, redistributes_.back(), nodomains,criterion);
966 MatrixArgs args(*redistMat, *redistComm);
967 BIGINT unknowns = redistMat->N();
968 unknowns = infoLevel->communicator().sum(unknowns);
970 if(redistComm->communicator().rank()==0 && criterion.debugLevel()>1){
971 double dunknowns= unknowns.todouble();
972 std::cout<<
"Level "<<level<<
" redistributed has "<<dunknowns<<
" unknowns, "<<dunknowns/redistComm->communicator().size()
973 <<
" unknowns per proc (procs="<<redistComm->communicator().size()<<
")"<<std::endl;
976 infoLevel.addRedistributed(redistComm);
977 infoLevel->freeGlobalLookup();
980 int levels = matrices_.levels();
981 maxlevels_ = parallelInformation_.finest()->communicator().max(levels);
982 assert(matrices_.levels()==redistributes_.size());
983 if(hasCoarsest() && rank==0 && criterion.debugLevel()>1)
984 std::cout<<
"operator complexity: "<<allnonzeros.todouble()/finenonzeros.todouble()<<std::endl;
988 template<
class M,
class IS,
class A>
995 template<
class M,
class IS,
class A>
999 return parallelInformation_;
1002 template<
class M,
class IS,
class A>
1005 int levels=aggregatesMaps().size();
1006 int maxlevels=parallelInformation_.finest()->communicator().max(levels);
1007 std::size_t size=(*(aggregatesMaps().begin()))->noVertices();
1009 std::vector<std::size_t> tmp;
1010 std::vector<std::size_t> *coarse, *fine;
1027 if(levels==maxlevels){
1028 const AggregatesMap& map = *(*(++aggregatesMaps().rbegin()));
1033 m=std::max(*iter,m);
1035 coarse->resize(m+1);
1037 srand((
unsigned)std::clock());
1038 std::set<size_t> used;
1039 for(
typename std::vector<std::size_t>::iterator iter=coarse->begin(); iter != coarse->end();
1042 std::pair<std::set<std::size_t>::iterator,
bool> ibpair
1043 = used.insert(static_cast<std::size_t>((((
double)rand())/(RAND_MAX+1.0)))*coarse->size());
1045 while(!ibpair.second)
1046 ibpair = used.insert(static_cast<std::size_t>((((
double)rand())/(RAND_MAX+1.0))*coarse->size()));
1047 *iter=*(ibpair.first);
1055 for(
typename AggregatesMapList::const_reverse_iterator aggregates=++aggregatesMaps().rbegin();
1056 aggregates != aggregatesMaps().rend(); ++aggregates,--levels){
1058 fine->resize((*aggregates)->noVertices());
1059 fine->assign(fine->size(), 0);
1061 ::prolongate(*(*aggregates), *coarse, *fine, static_cast<std::size_t>(1), *pinfo);
1063 std::swap(coarse, fine);
1067 assert(coarse==&data);
1070 template<
class M,
class IS,
class A>
1074 return aggregatesMaps_;
1076 template<
class M,
class IS,
class A>
1080 return redistributes_;
1083 template<
class M,
class IS,
class A>
1086 typedef typename AggregatesMapList::reverse_iterator AggregatesMapIterator;
1090 AggregatesMapIterator amap = aggregatesMaps_.rbegin();
1091 InfoIterator info = parallelInformation_.coarsest();
1092 for(Iterator level=matrices_.coarsest(), finest=matrices_.finest(); level != finest; --level, --info, ++amap){
1095 delete &level->getmat();
1096 if(level.isRedistributed())
1097 delete &(level.getRedistributed().getmat());
1102 template<
class M,
class IS,
class A>
1103 template<
class V,
class TA>
1106 assert(hierarchy.levels()==1);
1108 typedef typename RedistributeInfoList::const_iterator RIter;
1109 RIter redist = redistributes_.begin();
1111 Iterator matrix = matrices_.finest(), coarsest = matrices_.coarsest();
1113 if(redist->isSetup())
1114 hierarchy.addRedistributedOnCoarsest(matrix.getRedistributed().getmat().N());
1115 Dune::dvverb<<
"Level "<<level<<
" has "<<matrices_.finest()->getmat().N()<<
" unknowns!"<<std::endl;
1117 while(matrix != coarsest){
1118 ++matrix; ++level; ++redist;
1119 Dune::dvverb<<
"Level "<<level<<
" has "<<matrix->getmat().N()<<
" unknowns!"<<std::endl;
1121 hierarchy.addCoarser(matrix->getmat().N());
1122 if(redist->isSetup())
1123 hierarchy.addRedistributedOnCoarsest(matrix.getRedistributed().getmat().N());
1129 template<
class M,
class IS,
class A>
1130 template<
class S,
class TA>
1134 assert(smoothers.
levels()==0);
1137 typedef typename AggregatesMapList::const_iterator AggregatesIterator;
1140 cargs.setArgs(sargs);
1141 PinfoIterator pinfo = parallelInformation_.finest();
1142 AggregatesIterator aggregates = aggregatesMaps_.begin();
1144 for(MatrixIterator matrix = matrices_.finest(), coarsest = matrices_.coarsest();
1145 matrix != coarsest; ++matrix, ++pinfo, ++aggregates, ++level){
1146 cargs.setMatrix(matrix->getmat(), **aggregates);
1147 cargs.setComm(*pinfo);
1150 if(maxlevels()>levels()){
1152 cargs.setMatrix(matrices_.coarsest()->getmat(), **aggregates);
1153 cargs.setComm(*pinfo);
1159 template<
class M,
class IS,
class A>
1163 typedef typename AggregatesMapList::iterator AggregatesMapIterator;
1167 AggregatesMapIterator amap = aggregatesMaps_.begin();
1169 InfoIterator info = parallelInformation_.finest();
1170 typename RedistributeInfoList::iterator riIter = redistributes_.begin();
1171 Iterator level = matrices_.finest(), coarsest=matrices_.coarsest();
1172 if(level.isRedistributed()){
1173 info->buildGlobalLookup(info->indexSet().size());
1175 const_cast<Matrix&>(level.getRedistributed().getmat()),
1176 *info,info.getRedistributed(), *riIter);
1177 info->freeGlobalLookup();
1180 for(; level!=coarsest; ++amap){
1181 const Matrix& fine = (level.isRedistributed()?level.getRedistributed():*level).getmat();
1185 productBuilder.
calculate(fine, *(*amap), const_cast<Matrix&>(level->getmat()), *info, copyFlags);
1186 if(level.isRedistributed()){
1187 info->buildGlobalLookup(info->indexSet().size());
1189 const_cast<Matrix&>(level.getRedistributed().getmat()), *info,
1190 info.getRedistributed(), *riIter);
1191 info->freeGlobalLookup();
1196 template<
class M,
class IS,
class A>
1199 return matrices_.levels();
1202 template<
class M,
class IS,
class A>
1208 template<
class M,
class IS,
class A>
1211 return levels()==maxlevels() &&
1212 (!matrices_.coarsest().isRedistributed() ||matrices_.coarsest()->getmat().N()>0);
1215 template<
class M,
class IS,
class A>
1221 template<
class T,
class A>
1223 : finest_(0), coarsest_(0), nonAllocated_(0), allocator_(), levels_(0)
1226 template<
class T,
class A>
1230 finest_ = allocator_.allocate(1,0);
1231 finest_->element_ = &first;
1232 finest_->redistributed_ = 0;
1233 nonAllocated_ = finest_;
1234 coarsest_ = finest_;
1235 coarsest_->coarser_ = coarsest_->finer_ = 0;
1239 template<
class T,
class A>
1243 Element* current = coarsest_;
1244 coarsest_ = coarsest_->finer_;
1245 if(current != nonAllocated_){
1246 if(current->redistributed_)
1250 allocator_.deallocate(current, 1);
1255 template<
class T,
class A>
1261 template<
class T,
class A>
1267 template<
class T,
class A>
1272 coarsest_ = allocator_.allocate(1,0);
1274 finest_ = coarsest_;
1275 coarsest_->finer_ = 0;
1277 coarsest_->coarser_ = allocator_.allocate(1,0);
1278 coarsest_->coarser_->finer_ = coarsest_;
1279 coarsest_ = coarsest_->coarser_;
1282 coarsest_->redistributed_ = 0;
1283 coarsest_->coarser_=0;
1288 template<
class T,
class A>
1293 finest_ = allocator_.allocate(1,0);
1295 coarsest_ = finest_;
1296 coarsest_->coarser_ = coarsest_->finer_ = 0;
1298 finest_->finer_ = allocator_.allocate(1,0);
1299 finest_->finer_->coarser_ = finest_;
1300 finest_ = finest_->finer_;
1307 template<
class T,
class A>
1313 template<
class T,
class A>
1319 template<
class T,
class A>
1325 template<
class T,
class A>