dune-istl
2.2.0
|
00001 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*- 00002 // vi: set ts=8 sw=2 et sts=2: 00003 // $Id: galerkin.hh 1553 2012-04-26 12:18:19Z mblatt $ 00004 #ifndef DUNE_GALERKIN_HH 00005 #define DUNE_GALERKIN_HH 00006 00007 #include"aggregates.hh" 00008 #include"pinfo.hh" 00009 #include<dune/common/poolallocator.hh> 00010 #include<dune/common/enumset.hh> 00011 #include<set> 00012 #include<limits> 00013 #include<algorithm> 00014 00015 namespace Dune 00016 { 00017 namespace Amg 00018 { 00030 template<class T> 00031 struct OverlapVertex 00032 { 00036 typedef T Aggregate; 00037 00041 typedef T Vertex; 00042 00046 Aggregate* aggregate; 00047 00051 Vertex vertex; 00052 }; 00053 00054 00055 00060 template<class M> 00061 class SparsityBuilder 00062 { 00063 public: 00069 SparsityBuilder(M& matrix); 00070 00071 void insert(const typename M::size_type& index); 00072 00073 void operator++(); 00074 00075 std::size_t minRowSize(); 00076 00077 std::size_t maxRowSize(); 00078 00079 std::size_t sumRowSize(); 00080 std::size_t index() 00081 { 00082 return row_.index(); 00083 } 00084 private: 00086 typename M::CreateIterator row_; 00088 std::size_t minRowSize_; 00090 std::size_t maxRowSize_; 00091 std::size_t sumRowSize_; 00092 #ifdef DUNE_ISTL_WITH_CHECKING 00093 bool diagonalInserted; 00094 #endif 00095 }; 00096 00097 class BaseGalerkinProduct 00098 { 00099 public: 00108 template<class M, class V, class I, class O> 00109 void calculate(const M& fine, const AggregatesMap<V>& aggregates, M& coarse, 00110 const I& pinfo, const O& copy); 00111 00112 }; 00113 00114 template<class T> 00115 class GalerkinProduct 00116 : public BaseGalerkinProduct 00117 { 00118 public: 00119 typedef T ParallelInformation; 00120 00131 template<class M, class G, class V, class Set> 00132 M* build(const M& fine, G& fineGraph, V& visitedMap, 00133 const ParallelInformation& pinfo, 00134 AggregatesMap<typename G::VertexDescriptor>& aggregates, 00135 const typename M::size_type& size, 00136 const Set& copy); 00137 private: 00138 00145 template<class G, class I, class Set> 00146 const OverlapVertex<typename G::VertexDescriptor>* 00147 buildOverlapVertices(const G& graph, const I& pinfo, 00148 AggregatesMap<typename G::VertexDescriptor>& aggregates, 00149 const Set& overlap, 00150 std::size_t& overlapCount); 00151 00152 template<class A> 00153 struct OVLess 00154 { 00155 bool operator()(const OverlapVertex<A>& o1, const OverlapVertex<A>& o2) 00156 { 00157 return *o1.aggregate < *o2.aggregate; 00158 } 00159 }; 00160 }; 00161 00162 template<> 00163 class GalerkinProduct<SequentialInformation> 00164 : public BaseGalerkinProduct 00165 { 00166 public: 00177 template<class M, class G, class V, class Set> 00178 M* build(const M& fine, G& fineGraph, V& visitedMap, 00179 const SequentialInformation& pinfo, 00180 const AggregatesMap<typename G::VertexDescriptor>& aggregates, 00181 const typename M::size_type& size, 00182 const Set& copy); 00183 }; 00184 00185 struct BaseConnectivityConstructor 00186 { 00187 template<class R, class G, class V> 00188 static void constructOverlapConnectivity(R& row, G& graph, V& visitedMap, 00189 const AggregatesMap<typename G::VertexDescriptor>& aggregates, 00190 const OverlapVertex<typename G::VertexDescriptor>*& seed, 00191 const OverlapVertex<typename G::VertexDescriptor>* overlapEnd); 00192 00196 template<class R, class G, class V> 00197 static void constructNonOverlapConnectivity(R& row, G& graph, V& visitedMap, 00198 const AggregatesMap<typename G::VertexDescriptor>& aggregates, 00199 const typename G::VertexDescriptor& seed); 00200 00201 00205 template<class G, class S, class V> 00206 class ConnectedBuilder 00207 { 00208 public: 00212 typedef G Graph; 00216 typedef typename Graph::ConstEdgeIterator ConstEdgeIterator; 00217 00221 typedef S Set; 00222 00226 typedef V VisitedMap; 00227 00231 typedef typename Graph::VertexDescriptor Vertex; 00232 00240 ConnectedBuilder(const AggregatesMap<Vertex>& aggregates, Graph& graph, 00241 VisitedMap& visitedMap, Set& connected); 00242 00247 void operator()(const ConstEdgeIterator& edge); 00248 00249 private: 00253 const AggregatesMap<Vertex>& aggregates_; 00254 00255 Graph& graph_; 00256 00260 VisitedMap& visitedMap_; 00261 00265 Set& connected_; 00266 }; 00267 00268 }; 00269 00270 template<class G, class T> 00271 struct ConnectivityConstructor: public BaseConnectivityConstructor 00272 { 00273 typedef typename G::VertexDescriptor Vertex; 00274 00275 template<class V, class O, class R> 00276 static void examine(G& graph, 00277 V& visitedMap, 00278 const T& pinfo, 00279 const AggregatesMap<Vertex>& aggregates, 00280 const O& overlap, 00281 const OverlapVertex<Vertex>* overlapVertices, 00282 const OverlapVertex<Vertex>* overlapEnd, 00283 R& row); 00284 }; 00285 00286 template<class G> 00287 struct ConnectivityConstructor<G,SequentialInformation> : public BaseConnectivityConstructor 00288 { 00289 typedef typename G::VertexDescriptor Vertex; 00290 00291 template<class V, class R> 00292 static void examine(G& graph, 00293 V& visitedMap, 00294 const SequentialInformation& pinfo, 00295 const AggregatesMap<Vertex>& aggregates, 00296 R& row); 00297 }; 00298 00299 template<class T> 00300 struct DirichletBoundarySetter 00301 { 00302 template<class M, class O> 00303 static void set(M& coarse, const T& pinfo, const O& copy); 00304 }; 00305 00306 template<> 00307 struct DirichletBoundarySetter<SequentialInformation> 00308 { 00309 template<class M, class O> 00310 static void set(M& coarse, const SequentialInformation& pinfo, const O& copy); 00311 }; 00312 00313 template<class R, class G, class V> 00314 void BaseConnectivityConstructor::constructNonOverlapConnectivity(R& row, G& graph, V& visitedMap, 00315 const AggregatesMap<typename G::VertexDescriptor>& aggregates, 00316 const typename G::VertexDescriptor& seed) 00317 { 00318 assert(row.index()==aggregates[seed]); 00319 row.insert(aggregates[seed]); 00320 ConnectedBuilder<G,R,V> conBuilder(aggregates, graph, visitedMap, row); 00321 typedef typename G::VertexDescriptor Vertex; 00322 typedef std::allocator<Vertex> Allocator; 00323 typedef SLList<Vertex,Allocator> VertexList; 00324 typedef typename AggregatesMap<Vertex>::DummyEdgeVisitor DummyVisitor; 00325 VertexList vlist; 00326 DummyVisitor dummy; 00327 aggregates.template breadthFirstSearch<true,false>(seed,aggregates[seed], graph, vlist, dummy, 00328 conBuilder, visitedMap); 00329 } 00330 00331 template<class R, class G, class V> 00332 void BaseConnectivityConstructor::constructOverlapConnectivity(R& row, G& graph, V& visitedMap, 00333 const AggregatesMap<typename G::VertexDescriptor>& aggregates, 00334 const OverlapVertex<typename G::VertexDescriptor>*& seed, 00335 const OverlapVertex<typename G::VertexDescriptor>* overlapEnd) 00336 { 00337 ConnectedBuilder<G,R,V> conBuilder(aggregates, graph, visitedMap, row); 00338 const typename G::VertexDescriptor aggregate=*seed->aggregate; 00339 00340 if (row.index()==*seed->aggregate){ 00341 while(seed != overlapEnd && aggregate == *seed->aggregate){ 00342 row.insert(*seed->aggregate); 00343 // Walk over all neighbours and add them to the connected array. 00344 visitNeighbours(graph, seed->vertex, conBuilder); 00345 // Mark vertex as visited 00346 put(visitedMap, seed->vertex, true); 00347 ++seed; 00348 } 00349 } 00350 } 00351 00352 template<class G, class S, class V> 00353 BaseConnectivityConstructor::ConnectedBuilder<G,S,V>::ConnectedBuilder(const AggregatesMap<Vertex>& aggregates, 00354 Graph& graph, VisitedMap& visitedMap, 00355 Set& connected) 00356 : aggregates_(aggregates), graph_(graph), visitedMap_(visitedMap), connected_(connected) 00357 {} 00358 00359 template<class G, class S, class V> 00360 void BaseConnectivityConstructor::ConnectedBuilder<G,S,V>::operator()(const ConstEdgeIterator& edge) 00361 { 00362 typedef typename G::VertexDescriptor Vertex; 00363 const Vertex& vertex = aggregates_[edge.target()]; 00364 assert(vertex!= AggregatesMap<Vertex>::UNAGGREGATED); 00365 if(vertex!= AggregatesMap<Vertex>::ISOLATED) 00366 connected_.insert(vertex); 00367 } 00368 00369 template<class T> 00370 template<class G, class I, class Set> 00371 const OverlapVertex<typename G::VertexDescriptor>* 00372 GalerkinProduct<T>::buildOverlapVertices(const G& graph, const I& pinfo, 00373 AggregatesMap<typename G::VertexDescriptor>& aggregates, 00374 const Set& overlap, 00375 std::size_t& overlapCount) 00376 { 00377 // count the overlap vertices. 00378 typedef typename G::ConstVertexIterator ConstIterator; 00379 typedef typename I::GlobalLookupIndexSet GlobalLookup; 00380 typedef typename GlobalLookup::IndexPair IndexPair; 00381 00382 const ConstIterator end = graph.end(); 00383 overlapCount = 0; 00384 00385 const GlobalLookup& lookup=pinfo.globalLookup(); 00386 00387 for(ConstIterator vertex=graph.begin(); vertex != end; ++vertex){ 00388 const IndexPair* pair = lookup.pair(*vertex); 00389 00390 if(pair!=0 && overlap.contains(pair->local().attribute())) 00391 ++overlapCount; 00392 } 00393 // Allocate space 00394 typedef typename G::VertexDescriptor Vertex; 00395 00396 OverlapVertex<Vertex>* overlapVertices = new OverlapVertex<Vertex>[overlapCount=0?1:overlapCount]; 00397 if(overlapCount==0) 00398 return overlapVertices; 00399 00400 // Initialize them 00401 overlapCount=0; 00402 for(ConstIterator vertex=graph.begin(); vertex != end; ++vertex){ 00403 const IndexPair* pair = lookup.pair(*vertex); 00404 00405 if(pair!=0 && overlap.contains(pair->local().attribute())){ 00406 overlapVertices[overlapCount].aggregate = &aggregates[pair->local()]; 00407 overlapVertices[overlapCount].vertex = pair->local(); 00408 ++overlapCount; 00409 } 00410 } 00411 00412 dverb << overlapCount<<" overlap vertices"<<std::endl; 00413 00414 std::sort(overlapVertices, overlapVertices+overlapCount, OVLess<Vertex>()); 00415 // due to the sorting the isolated aggregates (to be skipped) are at the end. 00416 00417 return overlapVertices; 00418 } 00419 00420 template<class G, class T> 00421 template<class V, class O, class R> 00422 void ConnectivityConstructor<G,T>::examine(G& graph, 00423 V& visitedMap, 00424 const T& pinfo, 00425 const AggregatesMap<Vertex>& aggregates, 00426 const O& overlap, 00427 const OverlapVertex<Vertex>* overlapVertices, 00428 const OverlapVertex<Vertex>* overlapEnd, 00429 R& row) 00430 { 00431 typedef typename T::GlobalLookupIndexSet GlobalLookup; 00432 const GlobalLookup& lookup = pinfo.globalLookup(); 00433 00434 typedef typename G::VertexIterator VertexIterator; 00435 00436 VertexIterator vend=graph.end(); 00437 00438 #ifdef DUNE_ISTL_WITH_CHECKING 00439 std::set<Vertex> examined; 00440 #endif 00441 00442 // The aggregates owned by the process have lower local indices 00443 // then those not owned. We process them in the first pass. 00444 // They represent the rows 0, 1, ..., n of the coarse matrix 00445 for(VertexIterator vertex = graph.begin(); vertex != vend; ++vertex) 00446 if(!get(visitedMap, *vertex)){ 00447 // In the first pass we only process owner nodes 00448 typedef typename GlobalLookup::IndexPair IndexPair; 00449 const IndexPair* pair = lookup.pair(*vertex); 00450 if(pair==0 || !overlap.contains(pair->local().attribute())){ 00451 #ifdef DUNE_ISTL_WITH_CHECKING 00452 assert(examined.find(aggregates[*vertex])==examined.end()); 00453 examined.insert(aggregates[*vertex]); 00454 #endif 00455 constructNonOverlapConnectivity(row, graph, visitedMap, aggregates, *vertex); 00456 00457 // only needed for ALU 00458 // (ghosts with same global id as owners on the same process) 00459 if (pinfo.getSolverCategory() == static_cast<int>(SolverCategory::nonoverlapping)){ 00460 if(overlapVertices != overlapEnd){ 00461 if(*overlapVertices->aggregate!=AggregatesMap<Vertex>::ISOLATED){ 00462 constructOverlapConnectivity(row, graph, visitedMap, aggregates, overlapVertices, overlapEnd);} 00463 else{ 00464 ++overlapVertices; 00465 } 00466 } 00467 } 00468 ++row; 00469 } 00470 } 00471 00472 dvverb<<"constructed "<<row.index()<<" non-overlapping rows"<<std::endl; 00473 00474 // Now come the aggregates not owned by use. 00475 // They represent the rows n+1, ..., N 00476 while(overlapVertices != overlapEnd) 00477 if(*overlapVertices->aggregate!=AggregatesMap<Vertex>::ISOLATED){ 00478 00479 #ifdef DUNE_ISTL_WITH_CHECKING 00480 typedef typename GlobalLookup::IndexPair IndexPair; 00481 const IndexPair* pair = lookup.pair(overlapVertices->vertex); 00482 assert(pair!=0 && overlap.contains(pair->local().attribute())); 00483 assert(examined.find(aggregates[overlapVertices->vertex])==examined.end()); 00484 examined.insert(aggregates[overlapVertices->vertex]); 00485 #endif 00486 constructOverlapConnectivity(row, graph, visitedMap, aggregates, overlapVertices, overlapEnd); 00487 ++row; 00488 }else{ 00489 ++overlapVertices; 00490 } 00491 } 00492 00493 template<class G> 00494 template<class V, class R> 00495 void ConnectivityConstructor<G,SequentialInformation>::examine(G& graph, 00496 V& visitedMap, 00497 const SequentialInformation& pinfo, 00498 const AggregatesMap<Vertex>& aggregates, 00499 R& row) 00500 { 00501 typedef typename G::VertexIterator VertexIterator; 00502 00503 VertexIterator vend=graph.end(); 00504 for(VertexIterator vertex = graph.begin(); vertex != vend; ++vertex){ 00505 if(!get(visitedMap, *vertex)){ 00506 constructNonOverlapConnectivity(row, graph, visitedMap, aggregates, *vertex); 00507 ++row; 00508 } 00509 } 00510 00511 } 00512 00513 template<class M> 00514 SparsityBuilder<M>::SparsityBuilder(M& matrix) 00515 : row_(matrix.createbegin()), 00516 minRowSize_(std::numeric_limits<std::size_t>::max()), 00517 maxRowSize_(0), sumRowSize_(0) 00518 { 00519 #ifdef DUNE_ISTL_WITH_CHECKING 00520 diagonalInserted = false; 00521 #endif 00522 } 00523 template<class M> 00524 std::size_t SparsityBuilder<M>::maxRowSize() 00525 { 00526 return maxRowSize_; 00527 } 00528 template<class M> 00529 std::size_t SparsityBuilder<M>::minRowSize() 00530 { 00531 return minRowSize_; 00532 } 00533 00534 template<class M> 00535 std::size_t SparsityBuilder<M>::sumRowSize() 00536 { 00537 return sumRowSize_; 00538 } 00539 template<class M> 00540 void SparsityBuilder<M>::operator++() 00541 { 00542 sumRowSize_ += row_.size(); 00543 minRowSize_=std::min(minRowSize_, row_.size()); 00544 maxRowSize_=std::max(maxRowSize_, row_.size()); 00545 ++row_; 00546 #ifdef DUNE_ISTL_WITH_CHECKING 00547 assert(diagonalInserted); 00548 diagonalInserted = false; 00549 #endif 00550 } 00551 00552 template<class M> 00553 void SparsityBuilder<M>::insert(const typename M::size_type& index) 00554 { 00555 row_.insert(index); 00556 #ifdef DUNE_ISTL_WITH_CHECKING 00557 diagonalInserted = diagonalInserted || row_.index()==index; 00558 #endif 00559 } 00560 00561 template<class T> 00562 template<class M, class G, class V, class Set> 00563 M* GalerkinProduct<T>::build(const M& fine, G& fineGraph, V& visitedMap, 00564 const ParallelInformation& pinfo, 00565 AggregatesMap<typename G::VertexDescriptor>& aggregates, 00566 const typename M::size_type& size, 00567 const Set& overlap) 00568 { 00569 00570 typedef OverlapVertex<typename G::VertexDescriptor> OverlapVertex; 00571 00572 std::size_t count; 00573 00574 const OverlapVertex* overlapVertices = buildOverlapVertices(fineGraph, 00575 pinfo, 00576 aggregates, 00577 overlap, 00578 count); 00579 M* coarseMatrix = new M(size, size, M::row_wise); 00580 00581 // Reset the visited flags of all vertices. 00582 // As the isolated nodes will be skipped we simply mark them as visited 00583 00584 typedef typename G::VertexIterator Vertex; 00585 Vertex vend = fineGraph.end(); 00586 for(Vertex vertex = fineGraph.begin(); vertex != vend; ++vertex){ 00587 assert(aggregates[*vertex] != AggregatesMap<typename G::VertexDescriptor>::UNAGGREGATED); 00588 put(visitedMap, *vertex, aggregates[*vertex]==AggregatesMap<typename G::VertexDescriptor>::ISOLATED); 00589 } 00590 00591 SparsityBuilder<M> sparsityBuilder(*coarseMatrix); 00592 00593 ConnectivityConstructor<G,T>::examine(fineGraph, visitedMap, pinfo, 00594 aggregates, overlap, 00595 overlapVertices, 00596 overlapVertices+count, 00597 sparsityBuilder); 00598 00599 dinfo<<pinfo.communicator().rank()<<": Matrix ("<<coarseMatrix->N()<<"x"<<coarseMatrix->M()<<" row: min="<<sparsityBuilder.minRowSize()<<" max=" 00600 <<sparsityBuilder.maxRowSize()<<" avg=" 00601 <<static_cast<double>(sparsityBuilder.sumRowSize())/coarseMatrix->N() 00602 <<std::endl; 00603 00604 delete[] overlapVertices; 00605 00606 //calculate(fine, aggregates, *coarse, overlap); 00607 00608 return coarseMatrix; 00609 } 00610 00611 template<class M, class G, class V, class Set> 00612 M* GalerkinProduct<SequentialInformation>::build(const M& fine, G& fineGraph, V& visitedMap, 00613 const SequentialInformation& pinfo, 00614 const AggregatesMap<typename G::VertexDescriptor>& aggregates, 00615 const typename M::size_type& size, 00616 const Set& overlap) 00617 { 00618 M* coarseMatrix = new M(size, size, M::row_wise); 00619 00620 // Reset the visited flags of all vertices. 00621 // As the isolated nodes will be skipped we simply mark them as visited 00622 00623 typedef typename G::VertexIterator Vertex; 00624 Vertex vend = fineGraph.end(); 00625 for(Vertex vertex = fineGraph.begin(); vertex != vend; ++vertex){ 00626 assert(aggregates[*vertex] != AggregatesMap<typename G::VertexDescriptor>::UNAGGREGATED); 00627 put(visitedMap, *vertex, aggregates[*vertex]==AggregatesMap<typename G::VertexDescriptor>::ISOLATED); 00628 } 00629 00630 SparsityBuilder<M> sparsityBuilder(*coarseMatrix); 00631 00632 ConnectivityConstructor<G,SequentialInformation>::examine(fineGraph, visitedMap, pinfo, 00633 aggregates, sparsityBuilder); 00634 dinfo<<"Matrix row: min="<<sparsityBuilder.minRowSize()<<" max=" 00635 <<sparsityBuilder.maxRowSize()<<" average=" 00636 <<static_cast<double>(sparsityBuilder.sumRowSize())/coarseMatrix->N()<<std::endl; 00637 return coarseMatrix; 00638 } 00639 00640 template<class M, class V, class P, class O> 00641 void BaseGalerkinProduct::calculate(const M& fine, const AggregatesMap<V>& aggregates, M& coarse, 00642 const P& pinfo, const O& copy) 00643 { 00644 coarse = static_cast<typename M::field_type>(0); 00645 00646 typedef typename M::ConstIterator RowIterator; 00647 RowIterator endRow = fine.end(); 00648 00649 for(RowIterator row = fine.begin(); row != endRow; ++row) 00650 if(aggregates[row.index()] != AggregatesMap<V>::ISOLATED){ 00651 assert(aggregates[row.index()]!=AggregatesMap<V>::UNAGGREGATED); 00652 typedef typename M::ConstColIterator ColIterator; 00653 ColIterator endCol = row->end(); 00654 00655 for(ColIterator col = row->begin(); col != endCol; ++col) 00656 if(aggregates[col.index()] != AggregatesMap<V>::ISOLATED){ 00657 assert(aggregates[row.index()]!=AggregatesMap<V>::UNAGGREGATED); 00658 coarse[aggregates[row.index()]][aggregates[col.index()]]+=*col; 00659 } 00660 } 00661 00662 // get the right diagonal matrix values on copy lines from owner processes 00663 typedef typename M::block_type BlockType; 00664 std::vector<BlockType> rowsize(coarse.N(),BlockType(0)); 00665 for (RowIterator row = coarse.begin(); row != coarse.end(); ++row) 00666 rowsize[row.index()]=coarse[row.index()][row.index()]; 00667 pinfo.copyOwnerToAll(rowsize,rowsize); 00668 for (RowIterator row = coarse.begin(); row != coarse.end(); ++row) 00669 coarse[row.index()][row.index()] = rowsize[row.index()]; 00670 00671 // don't set dirichlet boundaries for copy lines to make novlp case work, 00672 // the preconditioner yields slightly different results now. 00673 00674 // Set the dirichlet border 00675 //DirichletBoundarySetter<P>::template set<M>(coarse, pinfo, copy); 00676 00677 } 00678 00679 template<class T> 00680 template<class M, class O> 00681 void DirichletBoundarySetter<T>::set(M& coarse, const T& pinfo, const O& copy) 00682 { 00683 typedef typename T::ParallelIndexSet::const_iterator ConstIterator; 00684 ConstIterator end = pinfo.indexSet().end(); 00685 typedef typename M::block_type Block; 00686 Block identity=Block(0.0); 00687 for(typename Block::RowIterator b=identity.begin(); b != identity.end(); ++b) 00688 b->operator[](b.index())=1.0; 00689 00690 for(ConstIterator index = pinfo.indexSet().begin(); 00691 index != end; ++index){ 00692 if(copy.contains(index->local().attribute())){ 00693 typedef typename M::ColIterator ColIterator; 00694 typedef typename M::row_type Row; 00695 Row row = coarse[index->local()]; 00696 ColIterator cend = row.find(index->local()); 00697 ColIterator col = row.begin(); 00698 for(; col != cend; ++col) 00699 *col = 0; 00700 00701 cend = row.end(); 00702 00703 assert(col != cend); // There should be a diagonal entry 00704 *col = identity; 00705 00706 for(++col; col != cend; ++col) 00707 *col = 0; 00708 } 00709 } 00710 } 00711 00712 template<class M, class O> 00713 void DirichletBoundarySetter<SequentialInformation>::set(M& coarse, 00714 const SequentialInformation& pinfo, 00715 const O& overlap) 00716 { 00717 } 00718 00719 }// namespace Amg 00720 }// namespace Dune 00721 #endif