3 #ifndef DUNE_AMG_MATRIXHIERARCHY_HH
4 #define DUNE_AMG_MATRIXHIERARCHY_HH
15 #include <dune/common/parallel/indexset.hh>
57 template<
class M,
class PI,
class A=std::allocator<M> >
65 typedef typename MatrixOperator::matrix_type
Matrix;
83 typedef typename Allocator::template rebind<AggregatesMap*>::other
AAllocator;
92 typedef typename Allocator::template rebind<RedistributeInfoType>::other
RILAllocator;
103 std::shared_ptr<ParallelInformation> pinfo = std::make_shared<ParallelInformation>());
112 template<
typename O,
typename T>
113 void build(
const T& criterion);
129 template<
class V,
class BA,
class TA>
137 template<
class S,
class TA>
145 std::size_t
levels()
const;
226 template<
class Matrix,
bool pr
int>
233 static void stats(
const Matrix& matrix)
235 DUNE_UNUSED_PARAMETER(matrix);
239 template<
class Matrix>
240 struct MatrixStats<
Matrix,true>
249 min=std::numeric_limits<size_type>::max();
256 min=std::min(min, row.size());
257 max=std::max(max, row.size());
268 static void stats(
const Matrix& matrix)
270 calc c= for_each(matrix.begin(), matrix.end(), calc());
271 dinfo<<
"Matrix row: min="<<c.min<<
" max="<<c.max
272 <<
" average="<<
static_cast<double>(c.sum)/matrix.N()
312 template<
typename M,
typename C1>
314 std::shared_ptr<M> newMatrix,
316 std::shared_ptr<SequentialInformation>& newComm,
318 int nparts, C1& criterion)
320 DUNE_UNUSED_PARAMETER(origMatrix);
321 DUNE_UNUSED_PARAMETER(newMatrix);
322 DUNE_UNUSED_PARAMETER(origComm);
323 DUNE_UNUSED_PARAMETER(newComm);
324 DUNE_UNUSED_PARAMETER(ri);
325 DUNE_UNUSED_PARAMETER(nparts);
326 DUNE_UNUSED_PARAMETER(criterion);
327 DUNE_THROW(NotImplemented,
"Redistribution does not make sense in sequential code!");
331 template<
typename M,
typename C,
typename C1>
333 std::shared_ptr<M> newMatrix,
335 std::shared_ptr<C>& newComm,
337 int nparts, C1& criterion)
340 #ifdef AMG_REPART_ON_COMM_GRAPH
344 criterion.debugLevel()>1);
358 if(origComm.communicator().rank()==0)
359 std::cout<<
"Original matrix"<<std::endl;
360 origComm.communicator().barrier();
364 newComm, ri.getInterface(),
365 criterion.debugLevel()>1);
366 #endif // if else AMG_REPART
368 if(origComm.communicator().rank()==0 && criterion.debugLevel()>1)
369 std::cout<<
"Repartitioning took "<<time.elapsed()<<
" seconds."<<std::endl;
374 ri.checkInterface(origComm.indexSet(), newComm->indexSet(), origComm.communicator());
380 if(origComm.communicator().rank()==0)
381 std::cout<<
"Original matrix"<<std::endl;
382 origComm.communicator().barrier();
383 if(newComm->communicator().size()>0)
385 origComm.communicator().barrier();
388 if(origComm.communicator().rank()==0 && criterion.debugLevel()>1)
389 std::cout<<
"Redistributing matrix took "<<time.elapsed()<<
" seconds."<<std::endl;
390 return existentOnRedist;
394 template<
class M,
class IS,
class A>
396 std::shared_ptr<ParallelInformation> pinfo)
397 : matrices_(fineMatrix),
398 parallelInformation_(pinfo)
401 DUNE_THROW(
ISTLError,
"MatrixOperator and ParallelInformation must belong to the same category!");
404 template<
class M,
class IS,
class A>
405 template<
typename O,
typename T>
408 prolongDamp_ = criterion.getProlongationDampingFactor();
409 typedef O OverlapFlags;
415 typedef bigunsignedint<
sizeof(int)*8*noints> BIGINT;
417 MatIterator mlevel = matrices_.finest();
418 MatrixStats<typename M::matrix_type,MINIMAL_DEBUG_LEVEL<=INFO_DEBUG_LEVEL>::stats(mlevel->getmat());
420 PInfoIterator infoLevel = parallelInformation_.finest();
422 finenonzeros = infoLevel->communicator().sum(finenonzeros);
423 BIGINT allnonzeros = finenonzeros;
429 BIGINT unknowns = mlevel->getmat().N();
431 unknowns = infoLevel->communicator().sum(unknowns);
432 double dunknowns=unknowns.todouble();
433 infoLevel->buildGlobalLookup(mlevel->getmat().N());
436 for(; level < criterion.maxLevel(); ++level, ++mlevel) {
437 assert(matrices_.levels()==redistributes_.size());
438 rank = infoLevel->communicator().rank();
439 if(rank==0 && criterion.debugLevel()>1)
440 std::cout<<
"Level "<<level<<
" has "<<dunknowns<<
" unknowns, "<<dunknowns/infoLevel->communicator().size()
441 <<
" unknowns per proc (procs="<<infoLevel->communicator().size()<<
")"<<std::endl;
453 && dunknowns < 30*infoLevel->communicator().size()))
454 && infoLevel->communicator().size()>1 &&
455 dunknowns/infoLevel->communicator().size() <= criterion.coarsenTarget())
458 std::shared_ptr<Matrix> redistMat = std::make_shared<Matrix>();
459 std::shared_ptr<ParallelInformation> redistComm;
460 std::size_t nodomains = (std::size_t)std::ceil(dunknowns/(criterion.minAggregateSize()
461 *criterion.coarsenTarget()));
462 if( nodomains<=criterion.minAggregateSize()/2 ||
463 dunknowns <= criterion.coarsenTarget() )
466 bool existentOnNextLevel =
468 redistComm, redistributes_.back(), nodomains,
470 BIGINT unknownsRedist = redistMat->N();
471 unknownsRedist = infoLevel->communicator().sum(unknownsRedist);
472 dunknowns= unknownsRedist.todouble();
473 if(redistComm->communicator().rank()==0 && criterion.debugLevel()>1)
474 std::cout<<
"Level "<<level<<
" (redistributed) has "<<dunknowns<<
" unknowns, "<<dunknowns/redistComm->communicator().size()
475 <<
" unknowns per proc (procs="<<redistComm->communicator().size()<<
")"<<std::endl;
476 MatrixArgs args(redistMat, *redistComm);
478 assert(mlevel.isRedistributed());
479 infoLevel.addRedistributed(redistComm);
480 infoLevel->freeGlobalLookup();
482 if(!existentOnNextLevel)
487 matrix = &(mlevel.getRedistributed());
488 info = &(infoLevel.getRedistributed());
489 info->buildGlobalLookup(matrix->getmat().N());
492 rank = info->communicator().rank();
493 if(dunknowns <= criterion.coarsenTarget())
499 typedef typename GraphCreator::GraphTuple GraphTuple;
503 std::vector<bool> excluded(matrix->getmat().N(),
false);
505 GraphTuple graphs = GraphCreator::create(*matrix, excluded, *info, OverlapFlags());
509 aggregatesMaps_.push_back(aggregatesMap);
513 int noAggregates, isoAggregates, oneAggregates, skippedAggregates;
515 std::tie(noAggregates, isoAggregates, oneAggregates, skippedAggregates) =
516 aggregatesMap->
buildAggregates(matrix->getmat(), *(std::get<1>(graphs)), criterion, level==0);
518 if(rank==0 && criterion.debugLevel()>2)
519 std::cout<<
" Have built "<<noAggregates<<
" aggregates totally ("<<isoAggregates<<
" isolated aggregates, "<<
520 oneAggregates<<
" aggregates of one vertex, and skipped "<<
521 skippedAggregates<<
" aggregates)."<<std::endl;
525 int start, end, overlapStart, overlapEnd;
526 int procs=info->communicator().rank();
527 int n = UNKNOWNS/procs;
528 int bigger = UNKNOWNS%procs;
533 end = (rank+1)*(n+1);
535 start = bigger + rank * n;
536 end = bigger + (rank + 1) * n;
541 overlapStart = start - 1;
543 overlapStart = start;
546 overlapEnd = end + 1;
550 assert((UNKNOWNS)*(overlapEnd-overlapStart)==aggregatesMap->
noVertices());
551 for(
int j=0; j< UNKNOWNS; ++j)
552 for(
int i=0; i < UNKNOWNS; ++i)
554 if(i>=overlapStart && i<overlapEnd)
556 int no = (j/2)*((UNKNOWNS)/2)+i/2;
557 (*aggregatesMap)[j*(overlapEnd-overlapStart)+i-overlapStart]=no;
562 if(criterion.debugLevel()>1 && info->communicator().rank()==0)
563 std::cout<<
"aggregating finished."<<std::endl;
565 BIGINT gnoAggregates=noAggregates;
566 gnoAggregates = info->communicator().sum(gnoAggregates);
567 double dgnoAggregates = gnoAggregates.todouble();
569 BIGINT gnoAggregates=((UNKNOWNS)/2)*((UNKNOWNS)/2);
572 if(criterion.debugLevel()>2 && rank==0)
573 std::cout <<
"Building "<<dgnoAggregates<<
" aggregates took "<<watch.elapsed()<<
" seconds."<<std::endl;
575 if(dgnoAggregates==0 || dunknowns/dgnoAggregates<criterion.minCoarsenRate())
580 std::cerr <<
"Stopped coarsening because of rate breakdown "<<dunknowns<<
"/"<<dgnoAggregates
581 <<
"="<<dunknowns/dgnoAggregates<<
"<"
582 <<criterion.minCoarsenRate()<<std::endl;
584 std::cerr<<
"Could not build any aggregates. Probably no connected nodes."<<std::endl;
586 aggregatesMap->
free();
587 delete aggregatesMap;
588 aggregatesMaps_.pop_back();
590 if(criterion.accumulate() && mlevel.isRedistributed() && info->communicator().size()>1) {
594 delete &(mlevel.getRedistributed().getmat());
595 mlevel.deleteRedistributed();
596 delete &(infoLevel.getRedistributed());
597 infoLevel.deleteRedistributed();
598 redistributes_.back().resetSetup();
603 unknowns = noAggregates;
604 dunknowns = dgnoAggregates;
606 CommunicationArgs commargs(info->communicator(),info->category());
607 parallelInformation_.addCoarser(commargs);
611 typename PropertyMapTypeSelector<VertexVisitedTag,PropertiesGraph>::Type visitedMap =
617 *(std::get<1>(graphs)),
622 GraphCreator::free(graphs);
624 if(criterion.debugLevel()>2) {
626 std::cout<<
"Coarsening of index sets took "<<watch.elapsed()<<
" seconds."<<std::endl;
631 infoLevel->buildGlobalLookup(aggregates);
634 infoLevel->globalLookup());
637 if(criterion.debugLevel()>2) {
639 std::cout<<
"Communicating global aggregate numbers took "<<watch.elapsed()<<
" seconds."<<std::endl;
643 std::vector<bool>& visited=excluded;
645 typedef std::vector<bool>::iterator Iterator;
646 typedef IteratorPropertyMap<Iterator, IdentityMap> VisitedMap2;
647 Iterator end = visited.end();
648 for(Iterator iter= visited.begin(); iter != end; ++iter)
651 VisitedMap2 visitedMap2(visited.begin(), Dune::IdentityMap());
653 std::shared_ptr<typename MatrixOperator::matrix_type>
654 coarseMatrix(productBuilder.build(*(std::get<0>(graphs)), visitedMap2,
659 dverb<<
"Building of sparsity pattern took "<<watch.elapsed()<<std::endl;
661 info->freeGlobalLookup();
663 delete std::get<0>(graphs);
664 productBuilder.calculate(matrix->getmat(), *aggregatesMap, *coarseMatrix, *infoLevel, OverlapFlags());
666 if(criterion.debugLevel()>2) {
668 std::cout<<
"Calculation entries of Galerkin product took "<<watch.elapsed()<<
" seconds."<<std::endl;
672 allnonzeros = allnonzeros + infoLevel->communicator().sum(nonzeros);
673 MatrixArgs args(coarseMatrix, *infoLevel);
675 matrices_.addCoarser(args);
680 infoLevel->freeGlobalLookup();
684 aggregatesMaps_.push_back(aggregatesMap);
686 if(criterion.debugLevel()>0) {
687 if(level==criterion.maxLevel()) {
688 BIGINT unknownsLevel = mlevel->getmat().N();
689 unknownsLevel = infoLevel->communicator().sum(unknownsLevel);
690 double dunknownsLevel = unknownsLevel.todouble();
691 if(rank==0 && criterion.debugLevel()>1) {
692 std::cout<<
"Level "<<level<<
" has "<<dunknownsLevel<<
" unknowns, "<<dunknownsLevel/infoLevel->communicator().size()
693 <<
" unknowns per proc (procs="<<infoLevel->communicator().size()<<
")"<<std::endl;
698 if(criterion.accumulate() && !redistributes_.back().isSetup() &&
699 infoLevel->communicator().size()>1) {
700 #if HAVE_MPI && !HAVE_PARMETIS
702 infoLevel->communicator().rank()==0)
703 std::cerr<<
"Successive accumulation of data on coarse levels only works with ParMETIS installed."
704 <<
" Fell back to accumulation to one domain on coarsest level"<<std::endl;
708 std::shared_ptr<Matrix> redistMat = std::make_shared<Matrix>();
709 std::shared_ptr<ParallelInformation> redistComm;
713 redistComm, redistributes_.back(), nodomains,criterion);
714 MatrixArgs args(redistMat, *redistComm);
715 BIGINT unknownsRedist = redistMat->N();
716 unknownsRedist = infoLevel->communicator().sum(unknownsRedist);
718 if(redistComm->communicator().rank()==0 && criterion.debugLevel()>1) {
719 double dunknownsRedist = unknownsRedist.todouble();
720 std::cout<<
"Level "<<level<<
" redistributed has "<<dunknownsRedist<<
" unknowns, "<<dunknownsRedist/redistComm->communicator().size()
721 <<
" unknowns per proc (procs="<<redistComm->communicator().size()<<
")"<<std::endl;
724 infoLevel.addRedistributed(redistComm);
725 infoLevel->freeGlobalLookup();
728 int levels = matrices_.levels();
729 maxlevels_ = parallelInformation_.finest()->communicator().max(levels);
730 assert(matrices_.levels()==redistributes_.size());
731 if(hasCoarsest() && rank==0 && criterion.debugLevel()>1)
732 std::cout<<
"operator complexity: "<<allnonzeros.todouble()/finenonzeros.todouble()<<std::endl;
736 template<
class M,
class IS,
class A>
743 template<
class M,
class IS,
class A>
747 return parallelInformation_;
750 template<
class M,
class IS,
class A>
753 int levels=aggregatesMaps().size();
754 int maxlevels=parallelInformation_.finest()->communicator().max(levels);
755 std::size_t size=(*(aggregatesMaps().begin()))->noVertices();
757 std::vector<std::size_t> tmp;
758 std::vector<std::size_t> *coarse, *fine;
775 if(levels==maxlevels) {
776 const AggregatesMap& map = *(*(++aggregatesMaps().rbegin()));
785 srand((
unsigned)std::clock());
786 std::set<size_t> used;
787 for(
typename std::vector<std::size_t>::iterator iter=coarse->begin(); iter != coarse->end();
790 std::pair<std::set<std::size_t>::iterator,
bool> ibpair
791 = used.insert(
static_cast<std::size_t
>((((
double)rand())/(RAND_MAX+1.0)))*coarse->size());
793 while(!ibpair.second)
794 ibpair = used.insert(
static_cast<std::size_t
>((((
double)rand())/(RAND_MAX+1.0))*coarse->size()));
795 *iter=*(ibpair.first);
803 for(
typename AggregatesMapList::const_reverse_iterator aggregates=++aggregatesMaps().rbegin();
804 aggregates != aggregatesMaps().rend(); ++aggregates,--levels) {
806 fine->resize((*aggregates)->noVertices());
807 fine->assign(fine->size(), 0);
809 ::prolongateVector(*(*aggregates), *coarse, *fine,
static_cast<std::size_t
>(1), *pinfo);
811 std::swap(coarse, fine);
815 assert(coarse==&data);
818 template<
class M,
class IS,
class A>
822 return aggregatesMaps_;
824 template<
class M,
class IS,
class A>
828 return redistributes_;
831 template<
class M,
class IS,
class A>
834 typedef typename AggregatesMapList::reverse_iterator AggregatesMapIterator;
838 AggregatesMapIterator amap = aggregatesMaps_.rbegin();
839 InfoIterator info = parallelInformation_.coarsest();
840 for(Iterator level=matrices_.coarsest(), finest=matrices_.finest(); level != finest; --level, --info, ++amap) {
847 template<
class M,
class IS,
class A>
848 template<
class V,
class BA,
class TA>
851 assert(hierarchy.levels()==1);
853 typedef typename RedistributeInfoList::const_iterator RIter;
854 RIter redist = redistributes_.begin();
856 Iterator matrix = matrices_.finest(), coarsest = matrices_.coarsest();
858 if(redist->isSetup())
859 hierarchy.addRedistributedOnCoarsest(matrix.getRedistributed().getmat().N());
860 Dune::dvverb<<
"Level "<<level<<
" has "<<matrices_.finest()->getmat().N()<<
" unknowns!"<<std::endl;
862 while(matrix != coarsest) {
863 ++matrix; ++level; ++redist;
864 Dune::dvverb<<
"Level "<<level<<
" has "<<matrix->getmat().N()<<
" unknowns!"<<std::endl;
866 hierarchy.addCoarser(matrix->getmat().N());
867 if(redist->isSetup())
868 hierarchy.addRedistributedOnCoarsest(matrix.getRedistributed().getmat().N());
874 template<
class M,
class IS,
class A>
875 template<
class S,
class TA>
879 assert(smoothers.
levels()==0);
882 typedef typename AggregatesMapList::const_iterator AggregatesIterator;
885 cargs.setArgs(sargs);
886 PinfoIterator pinfo = parallelInformation_.finest();
887 AggregatesIterator aggregates = aggregatesMaps_.begin();
889 for(MatrixIterator matrix = matrices_.finest(), coarsest = matrices_.coarsest();
890 matrix != coarsest; ++matrix, ++pinfo, ++aggregates, ++level) {
891 cargs.setMatrix(matrix->getmat(), **aggregates);
892 cargs.setComm(*pinfo);
895 if(maxlevels()>levels()) {
897 cargs.setMatrix(matrices_.coarsest()->getmat(), **aggregates);
898 cargs.setComm(*pinfo);
904 template<
class M,
class IS,
class A>
908 typedef typename AggregatesMapList::iterator AggregatesMapIterator;
912 AggregatesMapIterator amap = aggregatesMaps_.begin();
914 InfoIterator info = parallelInformation_.finest();
915 typename RedistributeInfoList::iterator riIter = redistributes_.begin();
916 Iterator level = matrices_.finest(), coarsest=matrices_.coarsest();
917 if(level.isRedistributed()) {
918 info->buildGlobalLookup(level->getmat().N());
920 const_cast<Matrix&
>(level.getRedistributed().getmat()),
921 *info,info.getRedistributed(), *riIter);
922 info->freeGlobalLookup();
925 for(; level!=coarsest; ++amap) {
926 const Matrix& fine = (level.isRedistributed() ? level.getRedistributed() : *level).getmat();
930 productBuilder.
calculate(fine, *(*amap),
const_cast<Matrix&
>(level->getmat()), *info, copyFlags);
931 if(level.isRedistributed()) {
932 info->buildGlobalLookup(level->getmat().N());
934 const_cast<Matrix&
>(level.getRedistributed().getmat()), *info,
935 info.getRedistributed(), *riIter);
936 info->freeGlobalLookup();
941 template<
class M,
class IS,
class A>
944 return matrices_.levels();
947 template<
class M,
class IS,
class A>
953 template<
class M,
class IS,
class A>
956 return levels()==maxlevels() &&
957 (!matrices_.coarsest().isRedistributed() ||matrices_.coarsest()->getmat().N()>0);
960 template<
class M,
class IS,
class A>
970 #endif // end DUNE_AMG_MATRIXHIERARCHY_HH