3 #ifndef DUNE_ISTL_OWNEROVERLAPCOPY_HH
4 #define DUNE_ISTL_OWNEROVERLAPCOPY_HH
21 #include <dune/common/enumset.hh>
24 #include <dune/common/parallel/indexset.hh>
25 #include <dune/common/parallel/communicator.hh>
26 #include <dune/common/parallel/remoteindices.hh>
27 #include <dune/common/parallel/mpicommunication.hh>
32 #include <dune/common/parallel/communication.hh>
35 template<
int dim,
template<
class,
class>
class Comm>
74 template <
class G,
class L>
109 DUNE_THROW(
ISTLError,
"OwnerOverlapCopyCommunication: global index not in index set");
110 localindices.insert(x);
123 DUNE_THROW(
ISTLError,
"OwnerOverlapCopyCommunication: global index not in index set");
124 remoteindices.insert(x);
142 return remoteindices;
150 localindices.clear();
151 remoteindices.clear();
156 std::set<IndexTripel> localindices;
158 std::set<RemoteIndexTripel> remoteindices;
170 template <
class GlobalIdType,
class LocalIdType=
int>
171 class OwnerOverlapCopyCommunication
173 template<
typename M,
typename G,
typename L>
176 OwnerOverlapCopyCommunication<G,L>&,
181 typedef typename std::set<IndexTripel>::const_iterator localindex_iterator;
182 typedef typename std::set<RemoteIndexTripel>::const_iterator remoteindex_iterator;
184 typedef Dune::ParallelLocalIndex<AttributeSet> LI;
186 typedef Dune::ParallelIndexSet<GlobalIdType,LI,512> PIS;
187 typedef Dune::RemoteIndices<PIS> RI;
188 typedef Dune::RemoteIndexListModifier<PIS,typename RI::Allocator,false> RILM;
189 typedef typename RI::RemoteIndex RX;
190 typedef Dune::BufferedCommunicator BC;
191 typedef Dune::Interface IF;
192 typedef EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::owner> OwnerSet;
193 typedef EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::copy> CopySet;
194 typedef Combine<EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::owner>,EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::overlap>,AttributeSet> OwnerOverlapSet;
195 typedef Dune::AllSet<AttributeSet> AllSet;
201 struct CopyGatherScatter
203 typedef typename CommPolicy<T>::IndexedType V;
205 static V gather(
const T& a, std::size_t i)
210 static void scatter(T& a, V v, std::size_t i)
216 struct AddGatherScatter
218 typedef typename CommPolicy<T>::IndexedType V;
220 static V gather(
const T& a, std::size_t i)
225 static void scatter(T& a, V v, std::size_t i)
231 void buildOwnerOverlapToAllInterface ()
const
233 if (OwnerOverlapToAllInterfaceBuilt)
234 OwnerOverlapToAllInterface.free();
235 typedef Combine<EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::owner>,EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::overlap>,AttributeSet> OwnerOverlapSet;
236 typedef Combine<OwnerOverlapSet,EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::copy>,AttributeSet> AllSet;
237 OwnerOverlapSet sourceFlags;
239 OwnerOverlapToAllInterface.build(ri,sourceFlags,destFlags);
240 OwnerOverlapToAllInterfaceBuilt =
true;
243 void buildOwnerToAllInterface ()
const
245 if (OwnerToAllInterfaceBuilt)
246 OwnerToAllInterface.free();
247 OwnerSet sourceFlags;
249 OwnerToAllInterface.build(ri,sourceFlags,destFlags);
250 OwnerToAllInterfaceBuilt =
true;
253 void buildOwnerCopyToAllInterface ()
const
255 if (OwnerCopyToAllInterfaceBuilt)
256 OwnerCopyToAllInterface.free();
257 typedef Combine<EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::owner>,EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::copy>,AttributeSet> OwnerCopySet;
258 typedef Combine<OwnerCopySet,EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::overlap>,AttributeSet> AllSet;
259 OwnerCopySet sourceFlags;
261 OwnerCopyToAllInterface.build(ri,sourceFlags,destFlags);
262 OwnerCopyToAllInterfaceBuilt =
true;
265 void buildOwnerCopyToOwnerCopyInterface ()
const
267 if (OwnerCopyToOwnerCopyInterfaceBuilt)
268 OwnerCopyToOwnerCopyInterface.free();
269 typedef Combine<EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::owner>,EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::copy>,AttributeSet> OwnerCopySet;
270 OwnerCopySet sourceFlags;
271 OwnerCopySet destFlags;
272 OwnerCopyToOwnerCopyInterface.build(ri,sourceFlags,destFlags);
273 OwnerCopyToOwnerCopyInterfaceBuilt =
true;
276 void buildCopyToAllInterface ()
const
278 if (CopyToAllInterfaceBuilt)
279 CopyToAllInterface.free();
282 CopyToAllInterface.build(ri,sourceFlags,destFlags);
283 CopyToAllInterfaceBuilt =
true;
292 DUNE_DEPRECATED_MSG(
"The solver category can only be set in the constructor. This method is deprecated and will be removed after Dune 2.7")
293 setSolverCategory (SolverCategory::Category set) {
298 DUNE_DEPRECATED_MSG(
"This method is deprecated and will be removed after Dune 2.7, use category() instead.")
299 getSolverCategory ()
const {
311 const CollectiveCommunication<MPI_Comm>& communicator()
const
323 void copyOwnerToAll (
const T& source, T& dest)
const
325 if (!OwnerToAllInterfaceBuilt)
326 buildOwnerToAllInterface ();
328 communicator.template build<T>(OwnerToAllInterface);
329 communicator.template forward<CopyGatherScatter<T> >(source,dest);
340 void copyCopyToAll (
const T& source, T& dest)
const
342 if (!CopyToAllInterfaceBuilt)
343 buildCopyToAllInterface ();
345 communicator.template build<T>(CopyToAllInterface);
346 communicator.template forward<CopyGatherScatter<T> >(source,dest);
357 void addOwnerOverlapToAll (
const T& source, T& dest)
const
359 if (!OwnerOverlapToAllInterfaceBuilt)
360 buildOwnerOverlapToAllInterface ();
362 communicator.template build<T>(OwnerOverlapToAllInterface);
363 communicator.template forward<AddGatherScatter<T> >(source,dest);
374 void addOwnerCopyToAll (
const T& source, T& dest)
const
376 if (!OwnerCopyToAllInterfaceBuilt)
377 buildOwnerCopyToAllInterface ();
379 communicator.template build<T>(OwnerCopyToAllInterface);
380 communicator.template forward<AddGatherScatter<T> >(source,dest);
391 void addOwnerCopyToOwnerCopy (
const T& source, T& dest)
const
393 if (!OwnerCopyToOwnerCopyInterfaceBuilt)
394 buildOwnerCopyToOwnerCopyInterface ();
396 communicator.template build<T>(OwnerCopyToOwnerCopyInterface);
397 communicator.template forward<AddGatherScatter<T> >(source,dest);
409 template<
class T1,
class T2>
410 void dot (
const T1& x,
const T1& y, T2& result)
const
413 if (mask.size()!=
static_cast<typename std::vector<double>::size_type
>(x.size()))
415 mask.resize(x.size());
416 for (
typename std::vector<double>::size_type i=0; i<mask.size(); i++)
418 for (
typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
420 mask[i->local().local()] = 0;
424 for (
typename T1::size_type i=0; i<x.size(); i++)
425 result += x[i]*(y[i])*mask[i];
426 result = cc.sum(result);
436 typename FieldTraits<typename T1::field_type>::real_type norm (
const T1& x)
const
438 using real_type =
typename FieldTraits<typename T1::field_type>::real_type;
441 if (mask.size()!=
static_cast<typename std::vector<double>::size_type
>(x.size()))
443 mask.resize(x.size());
444 for (
typename std::vector<double>::size_type i=0; i<mask.size(); i++)
446 for (
typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
448 mask[i->local().local()] = 0;
450 auto result = real_type(0.0);
451 for (
typename T1::size_type i=0; i<x.size(); i++)
452 result += Impl::asVector(x[i]).two_norm2()*mask[i];
453 return sqrt(cc.sum(result));
456 typedef Dune::EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::copy> CopyFlags;
459 typedef Dune::ParallelIndexSet<GlobalIdType,LI,512> ParallelIndexSet;
462 typedef Dune::RemoteIndices<PIS> RemoteIndices;
466 typedef Dune::GlobalLookupIndexSet<ParallelIndexSet> GlobalLookupIndexSet;
472 const ParallelIndexSet& indexSet()
const
481 const RemoteIndices& remoteIndices()
const
490 ParallelIndexSet& indexSet()
500 RemoteIndices& remoteIndices()
505 void buildGlobalLookup()
508 if(pis.seqNo()==oldseqNo)
511 delete globalLookup_;
514 globalLookup_ =
new GlobalLookupIndexSet(pis);
515 oldseqNo = pis.seqNo();
518 void buildGlobalLookup(std::size_t size)
521 if(pis.seqNo()==oldseqNo)
524 delete globalLookup_;
526 globalLookup_ =
new GlobalLookupIndexSet(pis, size);
527 oldseqNo = pis.seqNo();
530 void freeGlobalLookup()
532 delete globalLookup_;
536 const GlobalLookupIndexSet& globalLookup()
const
538 assert(globalLookup_ != 0);
539 return *globalLookup_;
548 void project (T1& x)
const
550 for (
typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
552 x[i->local().local()] = 0;
564 OwnerOverlapCopyCommunication (MPI_Comm comm_,
566 bool freecomm_ =
false)
567 : comm(comm_), cc(comm_), pis(), ri(pis,pis,comm_),
568 OwnerToAllInterfaceBuilt(false), OwnerOverlapToAllInterfaceBuilt(false),
569 OwnerCopyToAllInterfaceBuilt(false), OwnerCopyToOwnerCopyInterfaceBuilt(false),
570 CopyToAllInterfaceBuilt(false), globalLookup_(0), category_(cat_),
583 : comm(MPI_COMM_WORLD), cc(MPI_COMM_WORLD), pis(), ri(pis,pis,MPI_COMM_WORLD),
584 OwnerToAllInterfaceBuilt(false), OwnerOverlapToAllInterfaceBuilt(false),
585 OwnerCopyToAllInterfaceBuilt(false), OwnerCopyToOwnerCopyInterfaceBuilt(false),
586 CopyToAllInterfaceBuilt(false), globalLookup_(0), category_(cat_), freecomm(false)
596 OwnerOverlapCopyCommunication (
const IndexInfoFromGrid<GlobalIdType, LocalIdType>& indexinfo,
599 bool freecomm_ =
false)
600 : comm(comm_), cc(comm_), OwnerToAllInterfaceBuilt(false),
601 OwnerOverlapToAllInterfaceBuilt(false), OwnerCopyToAllInterfaceBuilt(false),
602 OwnerCopyToOwnerCopyInterfaceBuilt(false), CopyToAllInterfaceBuilt(false),
603 globalLookup_(0), category_(cat_), freecomm(freecomm_)
607 for (localindex_iterator i=indexinfo.localIndices().begin(); i!=indexinfo.localIndices().end(); ++i)
621 ri.setIndexSets(pis,pis,cc);
622 if (indexinfo.remoteIndices().size()>0)
624 remoteindex_iterator i=indexinfo.remoteIndices().begin();
625 int p = std::get<0>(*i);
626 RILM modifier = ri.template getModifier<false,true>(p);
627 typename PIS::const_iterator pi=pis.begin();
628 for ( ; i!=indexinfo.remoteIndices().end(); ++i)
631 if (p!=std::get<0>(*i))
634 modifier = ri.template getModifier<false,true>(p);
639 while (pi->global()!=std::get<1>(*i) && pi!=pis.end())
642 DUNE_THROW(ISTLError,
"OwnerOverlapCopyCommunication: global index not in index set");
655 ri.template getModifier<false,true>(0);
660 ~OwnerOverlapCopyCommunication ()
663 if (OwnerToAllInterfaceBuilt) OwnerToAllInterface.free();
664 if (OwnerOverlapToAllInterfaceBuilt) OwnerOverlapToAllInterface.free();
665 if (OwnerCopyToAllInterfaceBuilt) OwnerCopyToAllInterface.free();
666 if (OwnerCopyToOwnerCopyInterfaceBuilt) OwnerCopyToOwnerCopyInterface.free();
667 if (CopyToAllInterfaceBuilt) CopyToAllInterface.free();
668 if (globalLookup_)
delete globalLookup_;
670 if(comm!=MPI_COMM_NULL)
676 int wasFinalized = 0;
677 MPI_Finalized( &wasFinalized );
680 MPI_Comm_free(&comm);
685 OwnerOverlapCopyCommunication (
const OwnerOverlapCopyCommunication&)
688 CollectiveCommunication<MPI_Comm> cc;
691 mutable IF OwnerToAllInterface;
692 mutable bool OwnerToAllInterfaceBuilt;
693 mutable IF OwnerOverlapToAllInterface;
694 mutable bool OwnerOverlapToAllInterfaceBuilt;
695 mutable IF OwnerCopyToAllInterface;
696 mutable bool OwnerCopyToAllInterfaceBuilt;
697 mutable IF OwnerCopyToOwnerCopyInterface;
698 mutable bool OwnerCopyToOwnerCopyInterfaceBuilt;
699 mutable IF CopyToAllInterface;
700 mutable bool CopyToAllInterfaceBuilt;
701 mutable std::vector<double> mask;
703 GlobalLookupIndexSet* globalLookup_;