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 #ifndef DUNE_MATRIXREDIST_HH 00004 #define DUNE_MATRIXREDIST_HH 00005 #include"repartition.hh" 00006 #include<dune/common/exceptions.hh> 00007 #include<dune/common/parallel/indexset.hh> 00008 #include<dune/istl/owneroverlapcopy.hh> 00014 namespace Dune 00015 { 00016 template<typename T> 00017 struct RedistributeInformation 00018 { 00019 bool isSetup() const 00020 { 00021 return false; 00022 } 00023 template<class D> 00024 void redistribute(const D& from, D& to) const 00025 {} 00026 00027 template<class D> 00028 void redistributeBackward(D& from, const D& to) const 00029 {} 00030 00031 void resetSetup() 00032 {} 00033 00034 void setNoRows(std::size_t size) 00035 {} 00036 00037 void setNoCopyRows(std::size_t size) 00038 {} 00039 00040 void setNoBackwardsCopyRows(std::size_t size) 00041 {} 00042 00043 std::size_t getRowSize(std::size_t index) const 00044 { 00045 return -1; 00046 } 00047 00048 std::size_t getCopyRowSize(std::size_t index) const 00049 { 00050 return -1; 00051 } 00052 00053 std::size_t getBackwardsCopyRowSize(std::size_t index) const 00054 { 00055 return -1; 00056 } 00057 00058 }; 00059 00060 #if HAVE_MPI 00061 template<typename T, typename T1> 00062 class RedistributeInformation<OwnerOverlapCopyCommunication<T,T1> > 00063 { 00064 public: 00065 typedef OwnerOverlapCopyCommunication<T,T1> Comm; 00066 00067 RedistributeInformation() 00068 : interface(), setup_(false) 00069 {} 00070 00071 RedistributeInterface& getInterface() 00072 { 00073 return interface; 00074 } 00075 template<typename IS> 00076 void checkInterface(const IS& source, 00077 const IS& target, MPI_Comm comm) 00078 { 00079 RemoteIndices<IS> *ri=new RemoteIndices<IS>(source, target, comm); 00080 ri->template rebuild<true>(); 00081 Interface inf; 00082 typename OwnerOverlapCopyCommunication<int>::OwnerSet flags; 00083 int rank; 00084 MPI_Comm_rank(MPI_COMM_WORLD, &rank); 00085 inf.free(); 00086 inf.build(*ri, flags, flags); 00087 00088 00089 #ifdef DEBUG_REPART 00090 if(inf!=interface){ 00091 00092 int rank; 00093 MPI_Comm_rank(MPI_COMM_WORLD, &rank); 00094 if(rank==0) 00095 std::cout<<"Interfaces do not match!"<<std::endl; 00096 std::cout<<rank<<": redist interface new :"<<inf<<std::endl; 00097 std::cout<<rank<<": redist interface :"<<interface<<std::endl; 00098 00099 throw "autsch!"; 00100 delete ri; 00101 }else 00102 00103 #endif 00104 delete ri; 00105 } 00106 void setSetup() 00107 { 00108 setup_=true; 00109 interface.strip(); 00110 } 00111 00112 void resetSetup() 00113 { 00114 setup_=false; 00115 } 00116 00117 template<class GatherScatter, class D> 00118 void redistribute(const D& from, D& to) const 00119 { 00120 BufferedCommunicator communicator; 00121 communicator.template build<D>(from,to, interface); 00122 communicator.template forward<GatherScatter>(from, to); 00123 communicator.free(); 00124 } 00125 template<class GatherScatter, class D> 00126 void redistributeBackward(D& from, const D& to) const 00127 { 00128 00129 BufferedCommunicator communicator; 00130 communicator.template build<D>(from,to, interface); 00131 communicator.template backward<GatherScatter>(from, to); 00132 communicator.free(); 00133 } 00134 00135 template<class D> 00136 void redistribute(const D& from, D& to) const 00137 { 00138 redistribute<CopyGatherScatter<D> >(from,to); 00139 } 00140 template<class D> 00141 void redistributeBackward(D& from, const D& to) const 00142 { 00143 redistributeBackward<CopyGatherScatter<D> >(from,to); 00144 } 00145 bool isSetup() const 00146 { 00147 return setup_; 00148 } 00149 00150 void reserve(std::size_t size) 00151 {} 00152 00153 std::size_t& getRowSize(std::size_t index) 00154 { 00155 return rowSize[index]; 00156 } 00157 00158 std::size_t getRowSize(std::size_t index) const 00159 { 00160 return rowSize[index]; 00161 } 00162 00163 std::size_t& getCopyRowSize(std::size_t index) 00164 { 00165 return copyrowSize[index]; 00166 } 00167 00168 std::size_t getCopyRowSize(std::size_t index) const 00169 { 00170 return copyrowSize[index]; 00171 } 00172 00173 std::size_t& getBackwardsCopyRowSize(std::size_t index) 00174 { 00175 return backwardscopyrowSize[index]; 00176 } 00177 00178 std::size_t getBackwardsCopyRowSize(std::size_t index) const 00179 { 00180 return backwardscopyrowSize[index]; 00181 } 00182 00183 void setNoRows(std::size_t rows) 00184 { 00185 rowSize.resize(rows, 0); 00186 } 00187 00188 void setNoCopyRows(std::size_t rows) 00189 { 00190 copyrowSize.resize(rows, 0); 00191 } 00192 00193 void setNoBackwardsCopyRows(std::size_t rows) 00194 { 00195 backwardscopyrowSize.resize(rows, 0); 00196 } 00197 00198 private: 00199 std::vector<std::size_t> rowSize; 00200 std::vector<std::size_t> copyrowSize; 00201 std::vector<std::size_t> backwardscopyrowSize; 00202 RedistributeInterface interface; 00203 bool setup_; 00204 }; 00205 00214 template<class M, class RI> 00215 struct CommMatrixRowSize 00216 { 00217 // Make the default communication policy work. 00218 typedef typename M::size_type value_type; 00219 typedef typename M::size_type size_type; 00220 00226 CommMatrixRowSize(const M& m_, RI& rowsize_) 00227 : matrix(m_), rowsize(rowsize_) 00228 {} 00229 const M& matrix; 00230 RI& rowsize; 00231 00232 }; 00233 00234 00243 template<class M, class I> 00244 struct CommMatrixSparsityPattern 00245 { 00246 typedef typename M::size_type size_type; 00247 00254 CommMatrixSparsityPattern(const M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_) 00255 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize() 00256 {} 00257 00265 CommMatrixSparsityPattern(const M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_, 00266 const std::vector<typename M::size_type>& rowsize_) 00267 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), sparsity(aggidxset_.size()), rowsize(&rowsize_) 00268 {} 00269 00276 void storeSparsityPattern(M& m) 00277 { 00278 // insert diagonal to overlap rows 00279 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator IIter; 00280 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet; 00281 std::size_t nnz=0; 00282 #ifdef DEBUG_REPART 00283 int rank; 00284 00285 MPI_Comm_rank(MPI_COMM_WORLD, &rank); 00286 #endif 00287 for(IIter i= aggidxset.begin(), end=aggidxset.end(); i!=end; ++i){ 00288 if(!OwnerSet::contains(i->local().attribute())){ 00289 #ifdef DEBUG_REPART 00290 std::cout<<rank<<" Inserting diagonal for"<<i->local()<<std::endl; 00291 #endif 00292 sparsity[i->local()].insert(i->local()); 00293 } 00294 00295 nnz+=sparsity[i->local()].size(); 00296 } 00297 assert( aggidxset.size()==sparsity.size()); 00298 00299 if(nnz>0){ 00300 m.setSize(aggidxset.size(), aggidxset.size(), nnz); 00301 m.setBuildMode(M::row_wise); 00302 typename M::CreateIterator citer=m.createbegin(); 00303 #ifdef DEBUG_REPART 00304 std::size_t idx=0; 00305 bool correct=true; 00306 Dune::GlobalLookupIndexSet<I> global(aggidxset); 00307 #endif 00308 typedef typename std::vector<std::set<size_type> >::const_iterator Iter; 00309 for(Iter i=sparsity.begin(), end=sparsity.end(); i!=end; ++i, ++citer) 00310 { 00311 typedef typename std::set<size_type>::const_iterator SIter; 00312 for(SIter si=i->begin(), send=i->end(); si!=send; ++si) 00313 citer.insert(*si); 00314 #ifdef DEBUG_REPART 00315 if(i->find(idx)==i->end()){ 00316 const typename I::IndexPair* gi=global.pair(idx); 00317 assert(gi); 00318 std::cout<<rank<<": row "<<idx<<" is missing a diagonal entry! global="<<gi->global()<<" attr="<<gi->local().attribute()<<" "<< 00319 OwnerSet::contains(gi->local().attribute())<< 00320 " row size="<<i->size()<<std::endl; 00321 correct=false; 00322 } 00323 ++idx; 00324 #endif 00325 } 00326 #ifdef DEBUG_REPART 00327 if(!correct) 00328 throw "bla"; 00329 #endif 00330 } 00331 } 00332 00340 void completeSparsityPattern(std::vector<std::set<size_type> > add_sparsity) 00341 { 00342 typedef typename std::vector<std::set<size_type> >::const_iterator Iter; 00343 for (unsigned int i = 0; i != sparsity.size(); ++i) { 00344 if (add_sparsity[i].size() != 0){ 00345 typedef std::set<size_type> Set; 00346 Set tmp_set; 00347 std::insert_iterator<Set> tmp_insert (tmp_set, tmp_set.begin()); 00348 std::set_union(add_sparsity[i].begin(), add_sparsity[i].end(), 00349 sparsity[i].begin(), sparsity[i].end(), tmp_insert); 00350 sparsity[i].swap(tmp_set); 00351 } 00352 } 00353 } 00354 00355 const M& matrix; 00356 typedef Dune::GlobalLookupIndexSet<I> LookupIndexSet; 00357 const Dune::GlobalLookupIndexSet<I>& idxset; 00358 const I& aggidxset; 00359 std::vector<std::set<size_type> > sparsity; 00360 const std::vector<size_type>* rowsize; 00361 }; 00362 00363 template<class M, class I> 00364 struct CommPolicy<CommMatrixSparsityPattern<M,I> > 00365 { 00366 typedef CommMatrixSparsityPattern<M,I> Type; 00367 00372 typedef typename I::GlobalIndex IndexedType; 00373 00375 typedef VariableSize IndexedTypeFlag; 00376 00377 static typename M::size_type getSize(const Type& t, std::size_t i) 00378 { 00379 if(!t.rowsize) 00380 return t.matrix[i].size(); 00381 else 00382 { 00383 assert((*t.rowsize)[i]>0); 00384 return (*t.rowsize)[i]; 00385 } 00386 } 00387 }; 00388 00395 template<class M, class I> 00396 struct CommMatrixRow 00397 { 00406 CommMatrixRow(M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_) 00407 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize() 00408 {} 00409 00413 CommMatrixRow(M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_, 00414 std::vector<size_t>& rowsize_) 00415 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize(&rowsize_) 00416 {} 00422 void setOverlapRowsToDirichlet() 00423 { 00424 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator Iter; 00425 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet; 00426 00427 for(Iter i= aggidxset.begin(), end=aggidxset.end(); i!=end; ++i) 00428 if(!OwnerSet::contains(i->local().attribute())){ 00429 // Set to Dirchlet 00430 typedef typename M::ColIterator CIter; 00431 for(CIter c=matrix[i->local()].begin(), cend= matrix[i->local()].end(); 00432 c!= cend; ++c) 00433 { 00434 *c=0; 00435 if(c.index()==i->local()){ 00436 typedef typename M::block_type::RowIterator RIter; 00437 for(RIter r=c->begin(), rend=c->end(); 00438 r != rend; ++r) 00439 (*r)[r.index()]=1; 00440 } 00441 } 00442 } 00443 } 00445 M& matrix; 00447 const Dune::GlobalLookupIndexSet<I>& idxset; 00449 const I& aggidxset; 00451 std::vector<size_t>* rowsize; // row sizes differ from sender side in overlap! 00452 }; 00453 00454 template<class M, class I> 00455 struct CommPolicy<CommMatrixRow<M,I> > 00456 { 00457 typedef CommMatrixRow<M,I> Type; 00458 00463 typedef std::pair<typename I::GlobalIndex,typename M::block_type> IndexedType; 00464 00466 typedef VariableSize IndexedTypeFlag; 00467 00468 static std::size_t getSize(const Type& t, std::size_t i) 00469 { 00470 if(!t.rowsize) 00471 return t.matrix[i].size(); 00472 else 00473 { 00474 assert((*t.rowsize)[i]>0); 00475 return (*t.rowsize)[i]; 00476 } 00477 } 00478 }; 00479 00480 template<class M, class I, class RI> 00481 struct MatrixRowSizeGatherScatter 00482 { 00483 typedef CommMatrixRowSize<M,RI> Container; 00484 00485 static const typename M::size_type gather(const Container& cont, std::size_t i) 00486 { 00487 return cont.matrix[i].size(); 00488 } 00489 static void scatter(Container& cont, const typename M::size_type& rowsize, std::size_t i) 00490 { 00491 assert(rowsize); 00492 cont.rowsize.getRowSize(i)=rowsize; 00493 } 00494 00495 }; 00496 00497 template<class M, class I, class RI> 00498 struct MatrixCopyRowSizeGatherScatter 00499 { 00500 typedef CommMatrixRowSize<M,RI> Container; 00501 00502 static const typename M::size_type gather(const Container& cont, std::size_t i) 00503 { 00504 return cont.matrix[i].size(); 00505 } 00506 static void scatter(Container& cont, const typename M::size_type& rowsize, std::size_t i) 00507 { 00508 assert(rowsize); 00509 if (rowsize > cont.rowsize.getCopyRowSize(i)) 00510 cont.rowsize.getCopyRowSize(i)=rowsize; 00511 } 00512 00513 }; 00514 00515 template<class M, class I> 00516 struct MatrixSparsityPatternGatherScatter 00517 { 00518 typedef typename I::GlobalIndex GlobalIndex; 00519 typedef CommMatrixSparsityPattern<M,I> Container; 00520 typedef typename M::ConstColIterator ColIter; 00521 00522 static ColIter col; 00523 static GlobalIndex numlimits; 00524 00525 static const GlobalIndex& gather(const Container& cont, std::size_t i, std::size_t j) 00526 { 00527 if(j==0) 00528 col=cont.matrix[i].begin(); 00529 else if (col!=cont.matrix[i].end()) 00530 ++col; 00531 00532 //copy communication: different row sizes for copy rows with the same global index 00533 //are possible. If all values of current matrix row are sent, send 00534 //std::numeric_limits<GlobalIndex>::max() 00535 //and receiver will ignore it. 00536 if (col==cont.matrix[i].end()){ 00537 numlimits = std::numeric_limits<GlobalIndex>::max(); 00538 return numlimits; 00539 } 00540 else { 00541 const typename I::IndexPair* index=cont.idxset.pair(col.index()); 00542 assert(index); 00543 // Only send index if col is no ghost 00544 if ( index->local().attribute() != 2) 00545 return index->global(); 00546 else { 00547 numlimits = std::numeric_limits<GlobalIndex>::max(); 00548 return numlimits; 00549 } 00550 } 00551 } 00552 static void scatter(Container& cont, const GlobalIndex& gi, std::size_t i, std::size_t j) 00553 { 00554 try{ 00555 if (gi != std::numeric_limits<GlobalIndex>::max()){ 00556 const typename I::IndexPair& ip=cont.aggidxset.at(gi); 00557 assert(ip.global()==gi); 00558 std::size_t col = ip.local(); 00559 cont.sparsity[i].insert(col); 00560 00561 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet; 00562 if(!OwnerSet::contains(ip.local().attribute())) 00563 // preserve symmetry for overlap 00564 cont.sparsity[col].insert(i); 00565 } 00566 } 00567 catch(Dune::RangeError er){ 00568 // Entry not present in the new index set. Ignore! 00569 #ifdef DEBUG_REPART 00570 typedef typename Container::LookupIndexSet GlobalLookup; 00571 typedef typename GlobalLookup::IndexPair IndexPair; 00572 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet; 00573 00574 GlobalLookup lookup(cont.aggidxset); 00575 const IndexPair* pi=lookup.pair(i); 00576 assert(pi); 00577 if(OwnerSet::contains(pi->local().attribute())){ 00578 int rank; 00579 MPI_Comm_rank(MPI_COMM_WORLD,&rank); 00580 std::cout<<rank<<cont.aggidxset<<std::endl; 00581 std::cout<<rank<<": row "<<i<<" (global="<<gi <<") not in index set for owner index "<<pi->global()<<std::endl; 00582 throw er; 00583 } 00584 #endif 00585 } 00586 } 00587 00588 }; 00589 template<class M, class I> 00590 typename MatrixSparsityPatternGatherScatter<M,I>::ColIter MatrixSparsityPatternGatherScatter<M,I>::col; 00591 00592 template<class M, class I> 00593 typename MatrixSparsityPatternGatherScatter<M,I>::GlobalIndex MatrixSparsityPatternGatherScatter<M,I>::numlimits; 00594 00595 00596 template<class M, class I> 00597 struct MatrixRowGatherScatter 00598 { 00599 typedef typename I::GlobalIndex GlobalIndex; 00600 typedef CommMatrixRow<M,I> Container; 00601 typedef typename M::ConstColIterator ColIter; 00602 typedef typename std::pair<GlobalIndex,typename M::block_type> Data; 00603 static ColIter col; 00604 static Data datastore; 00605 static GlobalIndex numlimits; 00606 00607 static const Data& gather(const Container& cont, std::size_t i, std::size_t j) 00608 { 00609 if(j==0) 00610 col=cont.matrix[i].begin(); 00611 else if (col!=cont.matrix[i].end()) 00612 ++col; 00613 // copy communication: different row sizes for copy rows with the same global index 00614 // are possible. If all values of current matrix row are sent, send 00615 // std::numeric_limits<GlobalIndex>::max() 00616 // and receiver will ignore it. 00617 if (col==cont.matrix[i].end()){ 00618 numlimits = std::numeric_limits<GlobalIndex>::max(); 00619 datastore = std::make_pair(numlimits,*col); 00620 return datastore; 00621 } 00622 else { 00623 // convert local column index to global index 00624 const typename I::IndexPair* index=cont.idxset.pair(col.index()); 00625 assert(index); 00626 // Store the data to prevent reference to temporary 00627 // Only send index if col is no ghost 00628 if ( index->local().attribute() != 2) 00629 datastore = std::make_pair(index->global(),*col); 00630 else { 00631 numlimits = std::numeric_limits<GlobalIndex>::max(); 00632 datastore = std::make_pair(numlimits,*col); 00633 } 00634 return datastore; 00635 } 00636 } 00637 static void scatter(Container& cont, const Data& data, std::size_t i, std::size_t j) 00638 { 00639 try{ 00640 if (data.first != std::numeric_limits<GlobalIndex>::max()){ 00641 typename M::size_type column=cont.aggidxset.at(data.first).local(); 00642 cont.matrix[i][column]=data.second; 00643 } 00644 } 00645 catch(Dune::RangeError er){ 00646 // This an overlap row and might therefore lack some entries! 00647 } 00648 00649 } 00650 }; 00651 00652 template<class M, class I> 00653 typename MatrixRowGatherScatter<M,I>::ColIter MatrixRowGatherScatter<M,I>::col; 00654 00655 template<class M, class I> 00656 typename MatrixRowGatherScatter<M,I>::Data MatrixRowGatherScatter<M,I>::datastore; 00657 00658 template<class M, class I> 00659 typename MatrixRowGatherScatter<M,I>::GlobalIndex MatrixRowGatherScatter<M,I>::numlimits; 00660 00661 template<typename M, typename C> 00662 void redistributeSparsityPattern(M& origMatrix, M& newMatrix, C& origComm, C& newComm, 00663 RedistributeInformation<C>& ri) 00664 { 00665 typename C::CopySet copyflags; 00666 typename C::OwnerSet ownerflags; 00667 typedef typename C::ParallelIndexSet IndexSet; 00668 typedef RedistributeInformation<C> RI; 00669 std::vector<typename M::size_type> rowsize(newComm.indexSet().size(), 0); 00670 std::vector<typename M::size_type> copyrowsize(newComm.indexSet().size(), 0); 00671 std::vector<typename M::size_type> backwardscopyrowsize(origComm.indexSet().size(), 0); 00672 00673 // get owner rowsizes 00674 CommMatrixRowSize<M,RI> commRowSize(origMatrix, ri); 00675 ri.template redistribute<MatrixRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize,commRowSize); 00676 00677 origComm.buildGlobalLookup(); 00678 00679 for (std::size_t i=0; i < newComm.indexSet().size(); i++){ 00680 rowsize[i] = ri.getRowSize(i); 00681 } 00682 // get sparsity pattern from owner rows 00683 CommMatrixSparsityPattern<M,IndexSet> 00684 origsp(origMatrix, origComm.globalLookup(), newComm.indexSet()); 00685 CommMatrixSparsityPattern<M,IndexSet> 00686 newsp(origMatrix, origComm.globalLookup(), newComm.indexSet(), rowsize); 00687 00688 ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp,newsp); 00689 00690 // build copy to owner interface to get missing matrix values for novlp case 00691 if (origComm.getSolverCategory() == SolverCategory::nonoverlapping){ 00692 RemoteIndices<IndexSet> *ris = new RemoteIndices<IndexSet>(origComm.indexSet(), 00693 newComm.indexSet(), 00694 origComm.communicator()); 00695 ris->template rebuild<true>(); 00696 00697 ri.getInterface().free(); 00698 ri.getInterface().build(*ris,copyflags,ownerflags); 00699 00700 // get copy rowsizes 00701 CommMatrixRowSize<M,RI> commRowSize_copy(origMatrix, ri); 00702 ri.template redistribute<MatrixCopyRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize_copy, 00703 commRowSize_copy); 00704 00705 for (std::size_t i=0; i < newComm.indexSet().size(); i++){ 00706 copyrowsize[i] = ri.getCopyRowSize(i); 00707 } 00708 //get copy rowsizes for sender 00709 ri.redistributeBackward(backwardscopyrowsize,copyrowsize); 00710 for (std::size_t i=0; i < origComm.indexSet().size(); i++){ 00711 ri.getBackwardsCopyRowSize(i) = backwardscopyrowsize[i]; 00712 } 00713 00714 // get sparsity pattern from copy rows 00715 CommMatrixSparsityPattern<M,IndexSet> origsp_copy(origMatrix, 00716 origComm.globalLookup(), 00717 newComm.indexSet(), 00718 backwardscopyrowsize); 00719 CommMatrixSparsityPattern<M,IndexSet> newsp_copy(origMatrix, origComm.globalLookup(), 00720 newComm.indexSet(), copyrowsize); 00721 ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp_copy, 00722 newsp_copy); 00723 00724 newsp.completeSparsityPattern(newsp_copy.sparsity); 00725 newsp.storeSparsityPattern(newMatrix); 00726 } 00727 else 00728 newsp.storeSparsityPattern(newMatrix); 00729 00730 #ifdef DUNE_ISTL_WITH_CHECKING 00731 // Check for symmetry 00732 int ret=0; 00733 typedef typename M::ConstRowIterator RIter; 00734 for(RIter row=newMatrix.begin(), rend=newMatrix.end(); row != rend; ++row){ 00735 typedef typename M::ConstColIterator CIter; 00736 for(CIter col=row->begin(), cend=row->end(); col!=cend; ++col) 00737 { 00738 try{ 00739 newMatrix[col.index()][row.index()]; 00740 }catch(Dune::ISTLError e){ 00741 std::cerr<<newComm.communicator().rank()<<": entry (" 00742 <<col.index()<<","<<row.index()<<") missing! for symmetry!"<<std::endl; 00743 ret=1; 00744 00745 } 00746 00747 } 00748 } 00749 00750 if(ret) 00751 DUNE_THROW(ISTLError, "Matrix not symmetric!"); 00752 #endif 00753 } 00754 00755 template<typename M, typename C> 00756 void redistributeMatrixEntries(M& origMatrix, M& newMatrix, C& origComm, C& newComm, 00757 RedistributeInformation<C>& ri) 00758 { 00759 typedef typename C::ParallelIndexSet IndexSet; 00760 typedef RedistributeInformation<C> RI; 00761 typename C::OwnerSet ownerflags; 00762 std::vector<typename M::size_type> rowsize(newComm.indexSet().size(), 0); 00763 std::vector<typename M::size_type> copyrowsize(newComm.indexSet().size(), 0); 00764 std::vector<typename M::size_type> backwardscopyrowsize(origComm.indexSet().size(), 0); 00765 00766 for (std::size_t i=0; i < newComm.indexSet().size(); i++){ 00767 rowsize[i] = ri.getRowSize(i); 00768 if (origComm.getSolverCategory() == SolverCategory::nonoverlapping){ 00769 copyrowsize[i] = ri.getCopyRowSize(i); 00770 } 00771 } 00772 00773 for (std::size_t i=0; i < origComm.indexSet().size(); i++) 00774 if (origComm.getSolverCategory() == SolverCategory::nonoverlapping) 00775 backwardscopyrowsize[i] = ri.getBackwardsCopyRowSize(i); 00776 00777 00778 if (origComm.getSolverCategory() == SolverCategory::nonoverlapping){ 00779 // fill sparsity pattern from copy rows 00780 CommMatrixRow<M,IndexSet> origrow_copy(origMatrix, origComm.globalLookup(), 00781 newComm.indexSet(), backwardscopyrowsize); 00782 CommMatrixRow<M,IndexSet> newrow_copy(newMatrix, origComm.globalLookup(), 00783 newComm.indexSet(),copyrowsize); 00784 ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow_copy, 00785 newrow_copy); 00786 ri.getInterface().free(); 00787 RemoteIndices<IndexSet> *ris = new RemoteIndices<IndexSet>(origComm.indexSet(), 00788 newComm.indexSet(), 00789 origComm.communicator()); 00790 ris->template rebuild<true>(); 00791 ri.getInterface().build(*ris,ownerflags,ownerflags); 00792 } 00793 00794 CommMatrixRow<M,IndexSet> 00795 origrow(origMatrix, origComm.globalLookup(), newComm.indexSet()); 00796 CommMatrixRow<M,IndexSet> 00797 newrow(newMatrix, origComm.globalLookup(), newComm.indexSet(),rowsize); 00798 ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow,newrow); 00799 if (origComm.getSolverCategory() != static_cast<int>(SolverCategory::nonoverlapping)) 00800 newrow.setOverlapRowsToDirichlet(); 00801 if(newMatrix.N()>0&&newMatrix.N()<20) 00802 printmatrix(std::cout, newMatrix, "redist", "row"); 00803 } 00804 00821 template<typename M, typename C> 00822 void redistributeMatrix(M& origMatrix, M& newMatrix, C& origComm, C& newComm, 00823 RedistributeInformation<C>& ri) 00824 { 00825 ri.setNoRows(newComm.indexSet().size()); 00826 ri.setNoCopyRows(newComm.indexSet().size()); 00827 ri.setNoBackwardsCopyRows(origComm.indexSet().size()); 00828 redistributeSparsityPattern(origMatrix, newMatrix, origComm, newComm, ri); 00829 redistributeMatrixEntries(origMatrix, newMatrix, origComm, newComm, ri); 00830 } 00831 #endif 00832 } 00833 #endif