dune-istl  2.2.0
hierarchy.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=8 sw=2 sts=2:
3 // $Id: hierarchy.hh 1610 2012-05-29 16:01:01Z sander $
4 #ifndef DUNE_AMGHIERARCHY_HH
5 #define DUNE_AMGHIERARCHY_HH
6 
7 #include<list>
8 #include<memory>
9 #include<limits>
10 #include<algorithm>
11 #include"aggregates.hh"
12 #include"graph.hh"
13 #include"galerkin.hh"
14 #include"renumberer.hh"
15 #include"graphcreator.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>
20 #include<dune/istl/bvector.hh>
21 #include<dune/common/parallel/indexset.hh>
31 
32 namespace Dune
33 {
34  namespace Amg
35  {
47  enum{
55  MAX_PROCESSES = 72000};
56 
64  template<typename T, typename A=std::allocator<T> >
65  class Hierarchy
66  {
67  public:
71  typedef T MemberType;
72 
73  template<typename T1, typename T2>
74  class LevelIterator;
75 
76  private:
80  struct Element
81  {
82  friend class LevelIterator<Hierarchy<T,A>, T>;
83  friend class LevelIterator<const Hierarchy<T,A>, const T>;
84 
86  Element* coarser_;
87 
89  Element* finer_;
90 
92  MemberType* element_;
93 
95  MemberType* redistributed_;
96  };
97  public:
98 // enum{
99 // /**
100 // * @brief If true only the method addCoarser will be usable
101 // * otherwise only the method addFiner will be usable.
102 // */
103 // coarsen = b
104 // };
105 
109  typedef typename A::template rebind<Element>::other Allocator;
110 
112 
117  Hierarchy(MemberType& first);
118 
122  Hierarchy();
123 
128  void addCoarser(Arguments& args);
129 
131 
136  void addFiner(Arguments& args);
137 
144  template<class C, class T1>
146  : public BidirectionalIteratorFacade<LevelIterator<C,T1>,T1,T1&>
147  {
148  friend class LevelIterator<typename remove_const<C>::type,
149  typename remove_const<T1>::type >;
150  friend class LevelIterator<const typename remove_const<C>::type,
151  const typename remove_const<T1>::type >;
152 
153  public:
156  : element_(0)
157  {}
158 
159  LevelIterator(Element* element)
160  : element_(element)
161  {}
162 
165  typename remove_const<T1>::type>& other)
166  : element_(other.element_)
167  {}
168 
171  const typename remove_const<T1>::type>& other)
172  : element_(other.element_)
173  {}
174 
179  typename remove_const<T1>::type>& other) const
180  {
181  return element_ == other.element_;
182  }
183 
187  bool equals(const LevelIterator<const typename remove_const<C>::type,
188  const typename remove_const<T1>::type>& other) const
189  {
190  return element_ == other.element_;
191  }
192 
194  T1& dereference() const
195  {
196  return *(element_->element_);
197  }
198 
200  void increment()
201  {
202  element_ = element_->coarser_;
203  }
204 
206  void decrement()
207  {
208  element_ = element_->finer_;
209  }
210 
215  bool isRedistributed() const
216  {
217  return element_->redistributed_;
218  }
219 
224  T1& getRedistributed() const
225  {
226  assert(element_->redistributed_);
227  return *element_->redistributed_;
228  }
229  void addRedistributed(T1* t)
230  {
231  element_->redistributed_ = t;
232  }
233 
235  {
236  element_->redistributed_ = 0;
237  }
238 
239  private:
240  Element* element_;
241  };
242 
244  typedef LevelIterator<Hierarchy<T,A>,T> Iterator;
245 
247  typedef LevelIterator<const Hierarchy<T,A>, const T> ConstIterator;
248 
253  Iterator finest();
254 
259  Iterator coarsest();
260 
261 
266  ConstIterator finest() const;
267 
272  ConstIterator coarsest() const;
273 
278  std::size_t levels() const;
279 
281  ~Hierarchy();
282 
283  private:
285  Element* finest_;
287  Element* coarsest_;
289  Element* nonAllocated_;
291  Allocator allocator_;
293  int levels_;
294  };
295 
302  template<class M, class PI, class A=std::allocator<M> >
304  {
305  public:
307  typedef M MatrixOperator;
308 
310  typedef typename MatrixOperator::matrix_type Matrix;
311 
314 
316  typedef A Allocator;
317 
320 
323 
326 
328  typedef typename Allocator::template rebind<AggregatesMap*>::other AAllocator;
329 
331  typedef std::list<AggregatesMap*,AAllocator> AggregatesMapList;
332 
335 
337  typedef typename Allocator::template rebind<RedistributeInfoType>::other RILAllocator;
338 
340  typedef std::list<RedistributeInfoType,RILAllocator> RedistributeInfoList;
341 
347  MatrixHierarchy(const MatrixOperator& fineMatrix,
349 
350 
352 
358  template<typename O, typename T>
359  void build(const T& criterion);
360 
368  template<class F>
369  void recalculateGalerkin(const F& copyFlags);
370 
375  template<class V, class TA>
376  void coarsenVector(Hierarchy<BlockVector<V,TA> >& hierarchy) const;
377 
383  template<class S, class TA>
384  void coarsenSmoother(Hierarchy<S,TA>& smoothers,
385  const typename SmootherTraits<S>::Arguments& args) const;
386 
391  std::size_t levels() const;
392 
397  std::size_t maxlevels() const;
398 
399  bool hasCoarsest() const;
400 
405  bool isBuilt() const;
406 
411  const ParallelMatrixHierarchy& matrices() const;
412 
418 
423  const AggregatesMapList& aggregatesMaps() const;
424 
431 
432 
433  typename MatrixOperator::field_type getProlongationDampingFactor() const
434  {
435  return prolongDamp_;
436  }
437 
448  void getCoarsestAggregatesOnFinest(std::vector<std::size_t>& data) const;
449 
450  private:
451  typedef typename ConstructionTraits<MatrixOperator>::Arguments MatrixArgs;
452  typedef typename ConstructionTraits<ParallelInformation>::Arguments CommunicationArgs;
454  AggregatesMapList aggregatesMaps_;
456  RedistributeInfoList redistributes_;
458  ParallelMatrixHierarchy matrices_;
460  ParallelInformationHierarchy parallelInformation_;
461 
463  bool built_;
464 
466  int maxlevels_;
467 
468  typename MatrixOperator::field_type prolongDamp_;
469 
473  template<class Matrix, bool print>
474  struct MatrixStats
475  {
476 
480  static void stats(const Matrix& matrix)
481  {}
482  };
483 
484  template<class Matrix>
485  struct MatrixStats<Matrix,true>
486  {
487  struct calc
488  {
489  typedef typename Matrix::size_type size_type;
490  typedef typename Matrix::row_type matrix_row;
491 
492  calc()
493  {
494  min=std::numeric_limits<size_type>::max();
495  max=0;
496  sum=0;
497  }
498 
499  void operator()(const matrix_row& row)
500  {
501  min=std::min(min, row.size());
502  max=std::max(max, row.size());
503  sum += row.size();
504  }
505 
509  };
513  static void stats(const Matrix& matrix)
514  {
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()
518  <<std::endl;
519  }
520  };
521  };
522 
526  template<class T>
527  class CoarsenCriterion : public T
528  {
529  public:
535 
546  CoarsenCriterion(int maxLevel=100, int coarsenTarget=1000, double minCoarsenRate=1.2,
547  double prolongDamp=1.6, AccumulationMode accumulate=successiveAccu)
548  : AggregationCriterion(Dune::Amg::Parameters(maxLevel, coarsenTarget, minCoarsenRate, prolongDamp, accumulate))
549  {}
550 
552  : AggregationCriterion(parms)
553  {}
554 
555  };
556 
557  template<typename M, typename C1>
558  bool repartitionAndDistributeMatrix(const M& origMatrix, M& newMatrix,
559  SequentialInformation& origSequentialInformationomm,
560  SequentialInformation*& newComm,
562  int nparts, C1& criterion)
563  {
564  DUNE_THROW(NotImplemented, "Redistribution does not make sense in sequential code!");
565  }
566 
567 
568  template<typename M, typename C, typename C1>
569  bool repartitionAndDistributeMatrix(const M& origMatrix, M& newMatrix, C& origComm, C*& newComm,
571  int nparts, C1& criterion)
572  {
573  Timer time;
574 #ifdef AMG_REPART_ON_COMM_GRAPH
575  // Done not repartition the matrix graph, but a graph of the communication scheme.
576  bool existentOnRedist=Dune::commGraphRepartition(origMatrix, origComm, nparts, newComm,
577  ri.getInterface(),
578  criterion.debugLevel()>1);
579 
580 #else
582  typedef Dune::Amg::PropertiesGraph<MatrixGraph,
585  IdentityMap,
586  IdentityMap> PropertiesGraph;
587  MatrixGraph graph(origMatrix);
588  PropertiesGraph pgraph(graph);
589  buildDependency(pgraph, origMatrix, criterion, false);
590 
591 #ifdef DEBUG_REPART
592  if(origComm.communicator().rank()==0)
593  std::cout<<"Original matrix"<<std::endl;
594  origComm.communicator().barrier();
595  printGlobalSparseMatrix(origMatrix, origComm, std::cout);
596 #endif
597  bool existentOnRedist=Dune::graphRepartition(pgraph, origComm, nparts,
598  newComm, ri.getInterface(),
599  criterion.debugLevel()>1);
600 #endif // if else AMG_REPART
601 
602  if(origComm.communicator().rank()==0 && criterion.debugLevel()>1)
603  std::cout<<"Repartitioning took "<<time.elapsed()<<" seconds."<<std::endl;
604 
605  ri.setSetup();
606 
607 #ifdef DEBUG_REPART
608  ri.checkInterface(origComm.indexSet(), newComm->indexSet(), origComm.communicator());
609 #endif
610 
611  redistributeMatrix(const_cast<M&>(origMatrix), newMatrix, origComm, *newComm, ri);
612 
613 #ifdef DEBUG_REPART
614  if(origComm.communicator().rank()==0)
615  std::cout<<"Original matrix"<<std::endl;
616  origComm.communicator().barrier();
617  if(newComm->communicator().size()>0)
618  printGlobalSparseMatrix(newMatrix, *newComm, std::cout);
619  origComm.communicator().barrier();
620 #endif
621 
622  if(origComm.communicator().rank()==0 && criterion.debugLevel()>1)
623  std::cout<<"Redistributing matrix took "<<time.elapsed()<<" seconds."<<std::endl;
624  return existentOnRedist;
625 
626  }
627 
628  template<typename M>
629  bool repartitionAndDistributeMatrix(M& origMatrix, M& newMatrix,
630  SequentialInformation& origComm,
631  SequentialInformation& newComm,
633  {
634  return true;
635  }
636 
637  template<class M, class IS, class A>
639  const ParallelInformation& pinfo)
640  : matrices_(const_cast<MatrixOperator&>(fineOperator)),
641  parallelInformation_(const_cast<ParallelInformation&>(pinfo))
642  {
643  dune_static_assert((static_cast<int>(MatrixOperator::category) ==
644  static_cast<int>(SolverCategory::sequential) ||
645  static_cast<int>(MatrixOperator::category) ==
646  static_cast<int>(SolverCategory::overlapping) ||
647  static_cast<int>(MatrixOperator::category) ==
648  static_cast<int>(SolverCategory::nonoverlapping)),
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!");
652 
653  }
654 
655  template<class M, class IS, class A>
656  template<typename O, typename T>
657  void MatrixHierarchy<M,IS,A>::build(const T& criterion)
658  {
659  prolongDamp_ = criterion.getProlongationDampingFactor();
660  typedef O OverlapFlags;
661  typedef typename ParallelMatrixHierarchy::Iterator MatIterator;
662  typedef typename ParallelInformationHierarchy::Iterator PInfoIterator;
663 
664  static const int noints=(Dune::Amg::MAX_PROCESSES/4096>0)?(Dune::Amg::MAX_PROCESSES/4096):1;
665 
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());
670 
671  PInfoIterator infoLevel = parallelInformation_.finest();
672  BIGINT finenonzeros=countNonZeros(mlevel->getmat());
673  finenonzeros = infoLevel->communicator().sum(finenonzeros);
674  BIGINT allnonzeros = finenonzeros;
675 
676 
677  int level = 0;
678  int rank = 0;
679 
680  BIGINT unknowns = mlevel->getmat().N();
681 
682  unknowns = infoLevel->communicator().sum(unknowns);
683  double dunknowns=unknowns.todouble();
684  infoLevel->buildGlobalLookup(mlevel->getmat().N());
685  redistributes_.push_back(RedistributeInfoType());
686 
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;
693 
694  MatrixOperator* matrix=&(*mlevel);
695  ParallelInformation* info =&(*infoLevel);
696 
697  if((
698 #if HAVE_PARMETIS
699  criterion.accumulate()==successiveAccu
700 #else
701  false
702 #endif
703  || (criterion.accumulate()==atOnceAccu
704  && dunknowns < 30*infoLevel->communicator().size()))
705  && infoLevel->communicator().size()>1 &&
706  dunknowns/infoLevel->communicator().size() <= criterion.coarsenTarget())
707  {
708  // accumulate to fewer processors
709  Matrix* redistMat= new Matrix();
710  ParallelInformation* redistComm=0;
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() )
715  nodomains=1;
716 
717  bool existentOnNextLevel =
718  repartitionAndDistributeMatrix(mlevel->getmat(), *redistMat, *infoLevel,
719  redistComm, redistributes_.back(), nodomains,
720  criterion);
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);
728  mlevel.addRedistributed(ConstructionTraits<MatrixOperator>::construct(args));
729  assert(mlevel.isRedistributed());
730  infoLevel.addRedistributed(redistComm);
731  infoLevel->freeGlobalLookup();
732 
733  if(!existentOnNextLevel)
734  // We do not hold any data on the redistributed partitioning
735  break;
736 
737  // Work on the redistributed Matrix from now on
738  matrix = &(mlevel.getRedistributed());
739  info = &(infoLevel.getRedistributed());
740  info->buildGlobalLookup(matrix->getmat().N());
741  }
742 
743  rank = info->communicator().rank();
744  if(dunknowns <= criterion.coarsenTarget())
745  // No further coarsening needed
746  break;
747 
748  typedef PropertiesGraphCreator<MatrixOperator> GraphCreator;
749  typedef typename GraphCreator::PropertiesGraph PropertiesGraph;
750  typedef typename GraphCreator::MatrixGraph MatrixGraph;
751  typedef typename GraphCreator::GraphTuple GraphTuple;
752 
753  typedef typename PropertiesGraph::VertexDescriptor Vertex;
754 
755  std::vector<bool> excluded(matrix->getmat().N(), false);
756 
757  GraphTuple graphs = GraphCreator::create(*matrix, excluded, *info, OverlapFlags());
758 
759  AggregatesMap* aggregatesMap=new AggregatesMap(get<1>(graphs)->maxVertex()+1);
760 
761  aggregatesMaps_.push_back(aggregatesMap);
762 
763  Timer watch;
764  watch.reset();
765  int noAggregates, isoAggregates, oneAggregates, skippedAggregates;
766 
767  tie(noAggregates, isoAggregates, oneAggregates, skippedAggregates) =
768  aggregatesMap->buildAggregates(matrix->getmat(), *(get<1>(graphs)), criterion, level==0);
769 
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;
774 #ifdef TEST_AGGLO
775  {
776  // calculate size of local matrix in the distributed direction
777  int start, end, overlapStart, overlapEnd;
778  int procs=info->communicator().rank();
779  int n = UNKNOWNS/procs; // number of unknowns per process
780  int bigger = UNKNOWNS%procs; // number of process with n+1 unknows
781 
782  // Compute owner region
783  if(rank<bigger){
784  start = rank*(n+1);
785  end = (rank+1)*(n+1);
786  }else{
787  start = bigger + rank * n;
788  end = bigger + (rank + 1) * n;
789  }
790 
791  // Compute overlap region
792  if(start>0)
793  overlapStart = start - 1;
794  else
795  overlapStart = start;
796 
797  if(end<UNKNOWNS)
798  overlapEnd = end + 1;
799  else
800  overlapEnd = end;
801 
802  assert((UNKNOWNS)*(overlapEnd-overlapStart)==aggregatesMap->noVertices());
803  for(int j=0; j< UNKNOWNS; ++j)
804  for(int i=0; i < UNKNOWNS; ++i)
805  {
806  if(i>=overlapStart && i<overlapEnd)
807  {
808  int no = (j/2)*((UNKNOWNS)/2)+i/2;
809  (*aggregatesMap)[j*(overlapEnd-overlapStart)+i-overlapStart]=no;
810  }
811  }
812  }
813 #endif
814  if(criterion.debugLevel()>1 && info->communicator().rank()==0)
815  std::cout<<"aggregating finished."<<std::endl;
816 
817  BIGINT gnoAggregates=noAggregates;
818  gnoAggregates = info->communicator().sum(gnoAggregates);
819  double dgnoAggregates = gnoAggregates.todouble();
820 #ifdef TEST_AGGLO
821  BIGINT gnoAggregates=((UNKNOWNS)/2)*((UNKNOWNS)/2);
822 #endif
823 
824  if(criterion.debugLevel()>2 && rank==0)
825  std::cout << "Building "<<dgnoAggregates<<" aggregates took "<<watch.elapsed()<<" seconds."<<std::endl;
826 
827  if(dgnoAggregates==0 || dunknowns/dgnoAggregates<criterion.minCoarsenRate())
828  {
829  if(rank==0)
830  {
831  if(dgnoAggregates>0)
832  std::cerr << "Stopped coarsening because of rate breakdown "<<dunknowns<<"/"<<dgnoAggregates
833  <<"="<<dunknowns/dgnoAggregates<<"<"
834  <<criterion.minCoarsenRate()<<std::endl;
835  else
836  std::cerr<< "Could not build any aggregates. Probably no connected nodes."<<std::endl;
837  }
838  aggregatesMap->free();
839  delete aggregatesMap;
840  aggregatesMaps_.pop_back();
841 
842  if(criterion.accumulate() && mlevel.isRedistributed() && info->communicator().size()>1){
843  // coarse level matrix was already redistributed, but to more than 1 process
844  // Therefore need to delete the redistribution. Further down it will
845  // then be redistributed to 1 process
846  delete &(mlevel.getRedistributed().getmat());
847  mlevel.deleteRedistributed();
848  delete &(infoLevel.getRedistributed());
849  infoLevel.deleteRedistributed();
850  redistributes_.back().resetSetup();
851  }
852 
853  break;
854  }
855  unknowns = noAggregates;
856  dunknowns = dgnoAggregates;
857 
858  CommunicationArgs commargs(info->communicator(),info->getSolverCategory());
859  parallelInformation_.addCoarser(commargs);
860 
861  ++infoLevel; // parallel information on coarse level
862 
863  typename PropertyMapTypeSelector<VertexVisitedTag,PropertiesGraph>::Type visitedMap =
864  get(VertexVisitedTag(), *(get<1>(graphs)));
865 
866  watch.reset();
868  ::coarsen(*info,
869  *(get<1>(graphs)),
870  visitedMap,
871  *aggregatesMap,
872  *infoLevel);
873 
874  GraphCreator::free(graphs);
875 
876  if(criterion.debugLevel()>2){
877  if(rank==0)
878  std::cout<<"Coarsening of index sets took "<<watch.elapsed()<<" seconds."<<std::endl;
879  }
880 
881  watch.reset();
882 
883  infoLevel->buildGlobalLookup(aggregates);
885  *info,
886  infoLevel->globalLookup());
887 
888 
889  if(criterion.debugLevel()>2){
890  if(rank==0)
891  std::cout<<"Communicating global aggregate numbers took "<<watch.elapsed()<<" seconds."<<std::endl;
892  }
893 
894  watch.reset();
895  std::vector<bool>& visited=excluded;
896 
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)
901  *iter=false;
902 
903  VisitedMap2 visitedMap2(visited.begin(), Dune::IdentityMap());
904 
905  typename MatrixOperator::matrix_type* coarseMatrix;
906 
907  coarseMatrix = productBuilder.build(matrix->getmat(), *(get<0>(graphs)), visitedMap2,
908  *info,
909  *aggregatesMap,
910  aggregates,
911  OverlapFlags());
912 
913  info->freeGlobalLookup();
914 
915  delete get<0>(graphs);
916  productBuilder.calculate(matrix->getmat(), *aggregatesMap, *coarseMatrix, *infoLevel, OverlapFlags());
917 
918  if(criterion.debugLevel()>2){
919  if(rank==0)
920  std::cout<<"Calculation of Galerkin product took "<<watch.elapsed()<<" seconds."<<std::endl;
921  }
922 
923  BIGINT nonzeros = countNonZeros(*coarseMatrix);
924  allnonzeros = allnonzeros + infoLevel->communicator().sum(nonzeros);
925  MatrixArgs args(*coarseMatrix, *infoLevel);
926 
927  matrices_.addCoarser(args);
928  redistributes_.push_back(RedistributeInfoType());
929  } // end level loop
930 
931 
932  infoLevel->freeGlobalLookup();
933 
934  built_=true;
935  AggregatesMap* aggregatesMap=new AggregatesMap(0);
936  aggregatesMaps_.push_back(aggregatesMap);
937 
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;
946  }
947  }
948  }
949 
950  if(criterion.accumulate() && !redistributes_.back().isSetup() &&
951  infoLevel->communicator().size()>1){
952 #if HAVE_MPI && !HAVE_PARMETIS
953  if(criterion.accumulate()==successiveAccu &&
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;
957 #endif
958 
959  // accumulate to fewer processors
960  Matrix* redistMat= new Matrix();
961  ParallelInformation* redistComm=0;
962  int nodomains = 1;
963 
964  repartitionAndDistributeMatrix(mlevel->getmat(), *redistMat, *infoLevel,
965  redistComm, redistributes_.back(), nodomains,criterion);
966  MatrixArgs args(*redistMat, *redistComm);
967  BIGINT unknowns = redistMat->N();
968  unknowns = infoLevel->communicator().sum(unknowns);
969 
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;
974  }
975  mlevel.addRedistributed(ConstructionTraits<MatrixOperator>::construct(args));
976  infoLevel.addRedistributed(redistComm);
977  infoLevel->freeGlobalLookup();
978  }
979 
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;
985 
986  }
987 
988  template<class M, class IS, class A>
991  {
992  return matrices_;
993  }
994 
995  template<class M, class IS, class A>
998  {
999  return parallelInformation_;
1000  }
1001 
1002  template<class M, class IS, class A>
1003  void MatrixHierarchy<M,IS,A>::getCoarsestAggregatesOnFinest(std::vector<std::size_t>& data) const
1004  {
1005  int levels=aggregatesMaps().size();
1006  int maxlevels=parallelInformation_.finest()->communicator().max(levels);
1007  std::size_t size=(*(aggregatesMaps().begin()))->noVertices();
1008  // We need an auxiliary vector for the consecutive prolongation.
1009  std::vector<std::size_t> tmp;
1010  std::vector<std::size_t> *coarse, *fine;
1011 
1012  // make sure the allocated space suffices.
1013  tmp.reserve(size);
1014  data.reserve(size);
1015 
1016  // Correctly assign coarse and fine for the first prolongation such that
1017  // we end up in data in the end.
1018  if(levels%2==0){
1019  coarse=&tmp;
1020  fine=&data;
1021  }else{
1022  coarse=&data;
1023  fine=&tmp;
1024  }
1025 
1026  // Number the unknowns on the coarsest level consecutively for each process.
1027  if(levels==maxlevels){
1028  const AggregatesMap& map = *(*(++aggregatesMaps().rbegin()));
1029  std::size_t m=0;
1030 
1031  for(typename AggregatesMap::const_iterator iter = map.begin(); iter != map.end(); ++iter)
1032  if(*iter< AggregatesMap::ISOLATED)
1033  m=std::max(*iter,m);
1034 
1035  coarse->resize(m+1);
1036  std::size_t i=0;
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();
1040  ++iter, ++i)
1041  {
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());
1044 
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);
1048  }
1049  }
1050 
1051  typename ParallelInformationHierarchy::Iterator pinfo = parallelInformation().coarsest();
1052  --pinfo;
1053 
1054  // Now consecutively project the numbers to the finest level.
1055  for(typename AggregatesMapList::const_reverse_iterator aggregates=++aggregatesMaps().rbegin();
1056  aggregates != aggregatesMaps().rend(); ++aggregates,--levels){
1057 
1058  fine->resize((*aggregates)->noVertices());
1059  fine->assign(fine->size(), 0);
1061  ::prolongate(*(*aggregates), *coarse, *fine, static_cast<std::size_t>(1), *pinfo);
1062  --pinfo;
1063  std::swap(coarse, fine);
1064  }
1065 
1066  // Assertion to check that we really projected to data on the last step.
1067  assert(coarse==&data);
1068  }
1069 
1070  template<class M, class IS, class A>
1073  {
1074  return aggregatesMaps_;
1075  }
1076  template<class M, class IS, class A>
1079  {
1080  return redistributes_;
1081  }
1082 
1083  template<class M, class IS, class A>
1085  {
1086  typedef typename AggregatesMapList::reverse_iterator AggregatesMapIterator;
1087  typedef typename ParallelMatrixHierarchy::Iterator Iterator;
1088  typedef typename ParallelInformationHierarchy::Iterator InfoIterator;
1089 
1090  AggregatesMapIterator amap = aggregatesMaps_.rbegin();
1091  InfoIterator info = parallelInformation_.coarsest();
1092  for(Iterator level=matrices_.coarsest(), finest=matrices_.finest(); level != finest; --level, --info, ++amap){
1093  (*amap)->free();
1094  delete *amap;
1095  delete &level->getmat();
1096  if(level.isRedistributed())
1097  delete &(level.getRedistributed().getmat());
1098  }
1099  delete *amap;
1100  }
1101 
1102  template<class M, class IS, class A>
1103  template<class V, class TA>
1105  {
1106  assert(hierarchy.levels()==1);
1107  typedef typename ParallelMatrixHierarchy::ConstIterator Iterator;
1108  typedef typename RedistributeInfoList::const_iterator RIter;
1109  RIter redist = redistributes_.begin();
1110 
1111  Iterator matrix = matrices_.finest(), coarsest = matrices_.coarsest();
1112  int level=0;
1113  if(redist->isSetup())
1114  hierarchy.addRedistributedOnCoarsest(matrix.getRedistributed().getmat().N());
1115  Dune::dvverb<<"Level "<<level<<" has "<<matrices_.finest()->getmat().N()<<" unknowns!"<<std::endl;
1116 
1117  while(matrix != coarsest){
1118  ++matrix; ++level; ++redist;
1119  Dune::dvverb<<"Level "<<level<<" has "<<matrix->getmat().N()<<" unknowns!"<<std::endl;
1120 
1121  hierarchy.addCoarser(matrix->getmat().N());
1122  if(redist->isSetup())
1123  hierarchy.addRedistributedOnCoarsest(matrix.getRedistributed().getmat().N());
1124 
1125  }
1126 
1127  }
1128 
1129  template<class M, class IS, class A>
1130  template<class S, class TA>
1132  const typename SmootherTraits<S>::Arguments& sargs) const
1133  {
1134  assert(smoothers.levels()==0);
1135  typedef typename ParallelMatrixHierarchy::ConstIterator MatrixIterator;
1136  typedef typename ParallelInformationHierarchy::ConstIterator PinfoIterator;
1137  typedef typename AggregatesMapList::const_iterator AggregatesIterator;
1138 
1139  typename ConstructionTraits<S>::Arguments cargs;
1140  cargs.setArgs(sargs);
1141  PinfoIterator pinfo = parallelInformation_.finest();
1142  AggregatesIterator aggregates = aggregatesMaps_.begin();
1143  int level=0;
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);
1148  smoothers.addCoarser(cargs);
1149  }
1150  if(maxlevels()>levels()){
1151  // This is not the globally coarsest level and therefore smoothing is needed
1152  cargs.setMatrix(matrices_.coarsest()->getmat(), **aggregates);
1153  cargs.setComm(*pinfo);
1154  smoothers.addCoarser(cargs);
1155  ++level;
1156  }
1157  }
1158 
1159  template<class M, class IS, class A>
1160  template<class F>
1162  {
1163  typedef typename AggregatesMapList::iterator AggregatesMapIterator;
1164  typedef typename ParallelMatrixHierarchy::Iterator Iterator;
1165  typedef typename ParallelInformationHierarchy::Iterator InfoIterator;
1166 
1167  AggregatesMapIterator amap = aggregatesMaps_.begin();
1168  BaseGalerkinProduct productBuilder;
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());
1174  redistributeMatrixEntries(const_cast<Matrix&>(level->getmat()),
1175  const_cast<Matrix&>(level.getRedistributed().getmat()),
1176  *info,info.getRedistributed(), *riIter);
1177  info->freeGlobalLookup();
1178  }
1179 
1180  for(; level!=coarsest; ++amap){
1181  const Matrix& fine = (level.isRedistributed()?level.getRedistributed():*level).getmat();
1182  ++level;
1183  ++info;
1184  ++riIter;
1185  productBuilder.calculate(fine, *(*amap), const_cast<Matrix&>(level->getmat()), *info, copyFlags);
1186  if(level.isRedistributed()){
1187  info->buildGlobalLookup(info->indexSet().size());
1188  redistributeMatrixEntries(const_cast<Matrix&>(level->getmat()),
1189  const_cast<Matrix&>(level.getRedistributed().getmat()), *info,
1190  info.getRedistributed(), *riIter);
1191  info->freeGlobalLookup();
1192  }
1193  }
1194  }
1195 
1196  template<class M, class IS, class A>
1198  {
1199  return matrices_.levels();
1200  }
1201 
1202  template<class M, class IS, class A>
1204  {
1205  return maxlevels_;
1206  }
1207 
1208  template<class M, class IS, class A>
1210  {
1211  return levels()==maxlevels() &&
1212  (!matrices_.coarsest().isRedistributed() ||matrices_.coarsest()->getmat().N()>0);
1213  }
1214 
1215  template<class M, class IS, class A>
1217  {
1218  return built_;
1219  }
1220 
1221  template<class T, class A>
1223  : finest_(0), coarsest_(0), nonAllocated_(0), allocator_(), levels_(0)
1224  {}
1225 
1226  template<class T, class A>
1228  : allocator_()
1229  {
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;
1236  levels_ = 1;
1237  }
1238 
1239  template<class T, class A>
1241  {
1242  while(coarsest_){
1243  Element* current = coarsest_;
1244  coarsest_ = coarsest_->finer_;
1245  if(current != nonAllocated_){
1246  if(current->redistributed_)
1247  ConstructionTraits<T>::deconstruct(current->redistributed_);
1248  ConstructionTraits<T>::deconstruct(current->element_);
1249  }
1250  allocator_.deallocate(current, 1);
1251  //coarsest_->coarser_ = 0;
1252  }
1253  }
1254 
1255  template<class T, class A>
1256  std::size_t Hierarchy<T,A>::levels() const
1257  {
1258  return levels_;
1259  }
1260 
1261  template<class T, class A>
1263  {
1264  coarsest_->redistributed_ = ConstructionTraits<MemberType>::construct(args);
1265  }
1266 
1267  template<class T, class A>
1269  {
1270  if(!coarsest_){
1271  assert(!finest_);
1272  coarsest_ = allocator_.allocate(1,0);
1273  coarsest_->element_ = ConstructionTraits<MemberType>::construct(args);
1274  finest_ = coarsest_;
1275  coarsest_->finer_ = 0;
1276  }else{
1277  coarsest_->coarser_ = allocator_.allocate(1,0);
1278  coarsest_->coarser_->finer_ = coarsest_;
1279  coarsest_ = coarsest_->coarser_;
1280  coarsest_->element_ = ConstructionTraits<MemberType>::construct(args);
1281  }
1282  coarsest_->redistributed_ = 0;
1283  coarsest_->coarser_=0;
1284  ++levels_;
1285  }
1286 
1287 
1288  template<class T, class A>
1290  {
1291  if(!finest_){
1292  assert(!coarsest_);
1293  finest_ = allocator_.allocate(1,0);
1294  finest_->element = ConstructionTraits<T>::construct(args);
1295  coarsest_ = finest_;
1296  coarsest_->coarser_ = coarsest_->finer_ = 0;
1297  }else{
1298  finest_->finer_ = allocator_.allocate(1,0);
1299  finest_->finer_->coarser_ = finest_;
1300  finest_ = finest_->finer_;
1301  finest_->finer = 0;
1302  finest_->element = ConstructionTraits<T>::construct(args);
1303  }
1304  ++levels_;
1305  }
1306 
1307  template<class T, class A>
1309  {
1310  return Iterator(finest_);
1311  }
1312 
1313  template<class T, class A>
1315  {
1316  return Iterator(coarsest_);
1317  }
1318 
1319  template<class T, class A>
1321  {
1322  return ConstIterator(finest_);
1323  }
1324 
1325  template<class T, class A>
1327  {
1328  return ConstIterator(coarsest_);
1329  }
1331  }// namespace Amg
1332 }// namespace Dune
1333 
1334 #endif