dune-istl  2.2.1
owneroverlapcopy.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set ts=4 sw=2 et sts=2:
3 // $Id: owneroverlapcopy.hh 1504 2011-10-31 09:26:48Z rebecca $
4 #ifndef DUNE_OWNEROVERLAPCOPY_HH
5 #define DUNE_OWNEROVERLAPCOPY_HH
6 
7 #include<new>
8 #include<iostream>
9 #include<vector>
10 #include<list>
11 #include<map>
12 #include<set>
13 
14 #include"cmath"
15 
16 // MPI header
17 #if HAVE_MPI
18 #include<mpi.h>
19 #endif
20 
21 
22 #include<dune/common/tuples.hh>
23 #include<dune/common/enumset.hh>
24 
25 #if HAVE_MPI
26 #include <dune/common/parallel/indexset.hh>
27 #include <dune/common/parallel/communicator.hh>
28 #include <dune/common/parallel/remoteindices.hh>
29 #include<dune/common/mpicollectivecommunication.hh>
30 #endif
31 
32 #include"solvercategory.hh"
33 #include"istlexception.hh"
34 #include<dune/common/collectivecommunication.hh>
35 
36 template<int dim, template<class,class> class Comm>
37 void testRedistributed(int s);
38 
39 
40 namespace Dune {
41 
58  {
59  enum AttributeSet {
60  owner=1, overlap=2, copy=3 };
61  };
62 
74  template <class G, class L>
76  {
77  public:
79  typedef G GlobalIdType;
80 
82  typedef L LocalIdType;
83 
90  typedef tuple<GlobalIdType,LocalIdType,int> IndexTripel;
97  typedef tuple<int,GlobalIdType,int> RemoteIndexTripel;
98 
104  void addLocalIndex (const IndexTripel& x)
105  {
106  if (get<2>(x)!=OwnerOverlapCopyAttributeSet::owner &&
109  DUNE_THROW(ISTLError,"OwnerOverlapCopyCommunication: global index not in index set");
110  localindices.insert(x);
111  }
112 
119  {
120  if (get<2>(x)!=OwnerOverlapCopyAttributeSet::owner &&
123  DUNE_THROW(ISTLError,"OwnerOverlapCopyCommunication: global index not in index set");
124  remoteindices.insert(x);
125  }
126 
131  const std::set<IndexTripel>& localIndices () const
132  {
133  return localindices;
134  }
135 
140  const std::set<RemoteIndexTripel>& remoteIndices () const
141  {
142  return remoteindices;
143  }
144 
148  void clear ()
149  {
150  localindices.clear();
151  remoteindices.clear();
152  }
153 
154  private:
156  std::set<IndexTripel> localindices;
158  std::set<RemoteIndexTripel> remoteindices;
159  };
160 
161 
162 #if HAVE_MPI
163 
170  template <class GlobalIdType, class LocalIdType=int>
172  {
173  template<typename M, typename G, typename L>
174  friend void loadMatrixMarket(M&,
175  const std::string&,
177  bool);
178  // used types
180  typedef typename IndexInfoFromGrid<GlobalIdType,LocalIdType>::RemoteIndexTripel RemoteIndexTripel;
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;
185  public:
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;
196  protected:
197 
198 
200  template<typename T>
202  {
203  typedef typename CommPolicy<T>::IndexedType V;
204 
205  static V gather(const T& a, std::size_t i)
206  {
207  return a[i];
208  }
209 
210  static void scatter(T& a, V v, std::size_t i)
211  {
212  a[i] = v;
213  }
214  };
215  template<typename T>
217  {
218  typedef typename CommPolicy<T>::IndexedType V;
219 
220  static V gather(const T& a, std::size_t i)
221  {
222  return a[i];
223  }
224 
225  static void scatter(T& a, V v, std::size_t i)
226  {
227  a[i] += v;
228  }
229  };
230 
232  {
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;
238  AllSet destFlags;
239  OwnerOverlapToAllInterface.build(ri,sourceFlags,destFlags);
240  OwnerOverlapToAllInterfaceBuilt = true;
241  }
242 
244  {
245  if (OwnerToAllInterfaceBuilt)
246  OwnerToAllInterface.free();
247  OwnerSet sourceFlags;
248  AllSet destFlags;
249  OwnerToAllInterface.build(ri,sourceFlags,destFlags);
250  OwnerToAllInterfaceBuilt = true;
251  }
252 
254  {
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;
260  AllSet destFlags;
261  OwnerCopyToAllInterface.build(ri,sourceFlags,destFlags);
262  OwnerCopyToAllInterfaceBuilt = true;
263  }
264 
266  {
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;
274  }
275 
277  {
278  if (CopyToAllInterfaceBuilt)
279  CopyToAllInterface.free();
280  CopySet sourceFlags;
281  AllSet destFlags;
282  CopyToAllInterface.build(ri,sourceFlags,destFlags);
283  CopyToAllInterfaceBuilt = true;
284  }
285 
286  public:
287 
293  category = set;
294  }
295 
302  return category;
303  }
304 
305  const CollectiveCommunication<MPI_Comm>& communicator() const
306  {
307  return cc;
308  }
309 
316  template<class T>
317  void copyOwnerToAll (const T& source, T& dest) const
318  {
319  if (!OwnerToAllInterfaceBuilt)
322  communicator.template build<T>(OwnerToAllInterface);
323  communicator.template forward<CopyGatherScatter<T> >(source,dest);
324  communicator.free();
325  }
326 
333  template<class T>
334  void copyCopyToAll (const T& source, T& dest) const
335  {
336  if (!CopyToAllInterfaceBuilt)
339  communicator.template build<T>(CopyToAllInterface);
340  communicator.template forward<CopyGatherScatter<T> >(source,dest);
341  communicator.free();
342  }
343 
350  template<class T>
351  void addOwnerOverlapToAll (const T& source, T& dest) const
352  {
353  if (!OwnerOverlapToAllInterfaceBuilt)
356  communicator.template build<T>(OwnerOverlapToAllInterface);
357  communicator.template forward<AddGatherScatter<T> >(source,dest);
358  communicator.free();
359  }
360 
367  template<class T>
368  void addOwnerCopyToAll (const T& source, T& dest) const
369  {
370  if (!OwnerCopyToAllInterfaceBuilt)
373  communicator.template build<T>(OwnerCopyToAllInterface);
374  communicator.template forward<AddGatherScatter<T> >(source,dest);
375  communicator.free();
376  }
377 
384  template<class T>
385  void addOwnerCopyToOwnerCopy (const T& source, T& dest) const
386  {
387  if (!OwnerCopyToOwnerCopyInterfaceBuilt)
390  communicator.template build<T>(OwnerCopyToOwnerCopyInterface);
391  communicator.template forward<AddGatherScatter<T> >(source,dest);
392  communicator.free();
393  }
394 
395 
403  template<class T1, class T2>
404  void dot (const T1& x, const T1& y, T2& result) const
405  {
406  // set up mask vector
407  if (mask.size()!=static_cast<typename std::vector<double>::size_type>(x.size()))
408  {
409  mask.resize(x.size());
410  for (typename std::vector<double>::size_type i=0; i<mask.size(); i++)
411  mask[i] = 1;
412  for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
413  if (i->local().attribute()!=OwnerOverlapCopyAttributeSet::owner)
414  mask[i->local().local()] = 0;
415  }
416  result = 0;
417 
418  for (typename T1::size_type i=0; i<x.size(); i++)
419  result += x[i]*(y[i])*mask[i];
420  result = cc.sum(result);
421  return;
422  }
423 
430  template<class T1>
431  double norm (const T1& x) const
432  {
433  // set up mask vector
434  if (mask.size()!=static_cast<typename std::vector<double>::size_type>(x.size()))
435  {
436  mask.resize(x.size());
437  for (typename std::vector<double>::size_type i=0; i<mask.size(); i++)
438  mask[i] = 1;
439  for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
440  if (i->local().attribute()!=OwnerOverlapCopyAttributeSet::owner)
441  mask[i->local().local()] = 0;
442  }
443  double result = 0;
444  for (typename T1::size_type i=0; i<x.size(); i++)
445  result += x[i].two_norm2()*mask[i];
446  return sqrt(cc.sum(result));
447  }
448 
449  typedef Dune::EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::copy> CopyFlags;
450 
452  typedef Dune::ParallelIndexSet<GlobalIdType,LI,512> ParallelIndexSet;
453 
455  typedef Dune::RemoteIndices<PIS> RemoteIndices;
456 
459  typedef Dune::GlobalLookupIndexSet<ParallelIndexSet> GlobalLookupIndexSet;
460 
465  const ParallelIndexSet& indexSet() const
466  {
467  return pis;
468  }
469 
475  {
476  return ri;
477  }
478 
484  {
485  return pis;
486  }
487 
488 
494  {
495  return ri;
496  }
497 
499  {
500  if(globalLookup_){
501  if(pis.seqNo()==oldseqNo)
502  // Nothing changed!
503  return;
504  delete globalLookup_;
505  }
506 
507  globalLookup_ = new GlobalLookupIndexSet(pis);
508  oldseqNo = pis.seqNo();
509  }
510 
511  void buildGlobalLookup(std::size_t size)
512  {
513  if(globalLookup_){
514  if(pis.seqNo()==oldseqNo)
515  // Nothing changed!
516  return;
517  delete globalLookup_;
518  }
519  globalLookup_ = new GlobalLookupIndexSet(pis, size);
520  oldseqNo = pis.seqNo();
521  }
522 
524  {
525  delete globalLookup_;
526  globalLookup_=0;
527  }
528 
530  {
531  assert(globalLookup_ != 0);
532  return *globalLookup_;
533  }
534 
540  template<class T1>
541  void project (T1& x) const
542  {
543  for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
544  if (i->local().attribute()==OwnerOverlapCopyAttributeSet::copy)
545  x[i->local().local()] = 0;
546  }
547 
559  bool freecomm_ = false)
560  : comm(comm_), cc(comm_), pis(), ri(pis,pis,comm_),
561  OwnerToAllInterfaceBuilt(false), OwnerOverlapToAllInterfaceBuilt(false),
562  OwnerCopyToAllInterfaceBuilt(false), OwnerCopyToOwnerCopyInterfaceBuilt(false),
563  CopyToAllInterfaceBuilt(false), globalLookup_(0), category(cat_),
564  freecomm(freecomm_)
565  {}
566 
576  : comm(MPI_COMM_WORLD), cc(MPI_COMM_WORLD), pis(), ri(pis,pis,MPI_COMM_WORLD),
577  OwnerToAllInterfaceBuilt(false), OwnerOverlapToAllInterfaceBuilt(false),
578  OwnerCopyToAllInterfaceBuilt(false), OwnerCopyToOwnerCopyInterfaceBuilt(false),
579  CopyToAllInterfaceBuilt(false), globalLookup_(0), category(cat_), freecomm(false)
580  {}
581 
590  MPI_Comm comm_,
592  bool freecomm_ = false)
593  : comm(comm_), cc(comm_), OwnerToAllInterfaceBuilt(false),
594  OwnerOverlapToAllInterfaceBuilt(false), OwnerCopyToAllInterfaceBuilt(false),
595  OwnerCopyToOwnerCopyInterfaceBuilt(false), CopyToAllInterfaceBuilt(false),
596  globalLookup_(0), category(cat_), freecomm(freecomm_)
597  {
598  // set up an ISTL index set
599  pis.beginResize();
600  for (localindex_iterator i=indexinfo.localIndices().begin(); i!=indexinfo.localIndices().end(); ++i)
601  {
602  if (get<2>(*i)==OwnerOverlapCopyAttributeSet::owner)
603  pis.add(get<0>(*i),LI(get<1>(*i),OwnerOverlapCopyAttributeSet::owner,true));
605  pis.add(get<0>(*i),LI(get<1>(*i),OwnerOverlapCopyAttributeSet::overlap,true));
606  if (get<2>(*i)==OwnerOverlapCopyAttributeSet::copy)
607  pis.add(get<0>(*i),LI(get<1>(*i),OwnerOverlapCopyAttributeSet::copy,true));
608 // std::cout << cc.rank() << ": adding index " << get<0>(*i) << " " << get<1>(*i) << " " << get<2>(*i) << std::endl;
609  }
610  pis.endResize();
611 
612  // build remote indices WITHOUT communication
613 // std::cout << cc.rank() << ": build remote indices" << std::endl;
614  ri.setIndexSets(pis,pis,cc);
615  if (indexinfo.remoteIndices().size()>0)
616  {
617  remoteindex_iterator i=indexinfo.remoteIndices().begin();
618  int p = get<0>(*i);
619  RILM modifier = ri.template getModifier<false,true>(p);
620  typename PIS::const_iterator pi=pis.begin();
621  for ( ; i!=indexinfo.remoteIndices().end(); ++i)
622  {
623  // handle processor change
624  if (p!=get<0>(*i))
625  {
626  p = get<0>(*i);
627  modifier = ri.template getModifier<false,true>(p);
628  pi=pis.begin();
629  }
630 
631  // position to correct entry in parallel index set
632  while (pi->global()!=get<1>(*i) && pi!=pis.end())
633  ++pi;
634  if (pi==pis.end())
635  DUNE_THROW(ISTLError,"OwnerOverlapCopyCommunication: global index not in index set");
636 
637  // insert entry
638 // std::cout << cc.rank() << ": adding remote index " << get<0>(*i) << " " << get<1>(*i) << " " << get<2>(*i) << std::endl;
639  if (get<2>(*i)==OwnerOverlapCopyAttributeSet::owner)
640  modifier.insert(RX(OwnerOverlapCopyAttributeSet::owner,&(*pi)));
642  modifier.insert(RX(OwnerOverlapCopyAttributeSet::overlap,&(*pi)));
643  if (get<2>(*i)==OwnerOverlapCopyAttributeSet::copy)
644  modifier.insert(RX(OwnerOverlapCopyAttributeSet::copy,&(*pi)));
645  }
646  }else{
647  // Force remote indices to be synced!
648  ri.template getModifier<false,true>(0);
649  }
650  }
651 
652  // destructor: free memory in some objects
654  {
655  ri.free();
656  if (OwnerToAllInterfaceBuilt) OwnerToAllInterface.free();
657  if (OwnerOverlapToAllInterfaceBuilt) OwnerOverlapToAllInterface.free();
658  if (OwnerCopyToAllInterfaceBuilt) OwnerCopyToAllInterface.free();
659  if (OwnerCopyToOwnerCopyInterfaceBuilt) OwnerCopyToOwnerCopyInterface.free();
660  if (CopyToAllInterfaceBuilt) CopyToAllInterface.free();
661  if (globalLookup_) delete globalLookup_;
662  if (freecomm==true)
663  if(comm!=MPI_COMM_NULL)
664  MPI_Comm_free(&comm);
665  }
666 
667  private:
669  {}
670  MPI_Comm comm;
671  CollectiveCommunication<MPI_Comm> cc;
672  PIS pis;
673  RI ri;
674  mutable IF OwnerToAllInterface;
675  mutable bool OwnerToAllInterfaceBuilt;
676  mutable IF OwnerOverlapToAllInterface;
677  mutable bool OwnerOverlapToAllInterfaceBuilt;
678  mutable IF OwnerCopyToAllInterface;
679  mutable bool OwnerCopyToAllInterfaceBuilt;
680  mutable IF OwnerCopyToOwnerCopyInterface;
681  mutable bool OwnerCopyToOwnerCopyInterfaceBuilt;
682  mutable IF CopyToAllInterface;
683  mutable bool CopyToAllInterfaceBuilt;
684  mutable std::vector<double> mask;
685  int oldseqNo;
686  GlobalLookupIndexSet* globalLookup_;
687  SolverCategory::Category category;
688  bool freecomm;
689  };
690 
691 #endif
692 
693 
696 } // end namespace
697 
698 #endif