Dune Core Modules (2.7.1)

novlpschwarz.hh
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ISTL_NOVLPSCHWARZ_HH
4 #define DUNE_ISTL_NOVLPSCHWARZ_HH
5 
6 #include <iostream> // for input/output to shell
7 #include <fstream> // for input/output to files
8 #include <vector> // STL vector class
9 #include <sstream>
10 
11 #include <cmath> // Yes, we do some math here
12 
13 #include <dune/common/timer.hh>
14 
15 #include <dune/common/hybridutilities.hh>
16 
17 #include "io.hh"
18 #include "bvector.hh"
19 #include "vbvector.hh"
20 #include "bcrsmatrix.hh"
21 #include "io.hh"
22 #include "gsetc.hh"
23 #include "ilu.hh"
24 #include "operators.hh"
25 #include "solvers.hh"
26 #include "preconditioners.hh"
27 #include "scalarproducts.hh"
28 #include "owneroverlapcopy.hh"
29 
30 namespace Dune {
31 
60  template<class M, class X, class Y, class C>
62  {
63  public:
65  typedef M matrix_type;
67  typedef X domain_type;
69  typedef Y range_type;
71  typedef typename X::field_type field_type;
73  typedef C communication_type;
74 
75  typedef typename C::PIS PIS;
76  typedef typename C::RI RI;
77  typedef typename RI::RemoteIndexList RIL;
78  typedef typename RI::const_iterator RIIterator;
79  typedef typename RIL::const_iterator RILIterator;
80  typedef typename M::ConstColIterator ColIterator;
81  typedef typename M::ConstRowIterator RowIterator;
82  typedef std::multimap<int,int> MM;
83  typedef std::multimap<int,std::pair<int,RILIterator> > RIMap;
84  typedef typename RIMap::iterator RIMapit;
85 
94  : _A_(stackobject_to_shared_ptr(A)), communication(com), buildcomm(true)
95  {}
96 
97  NonoverlappingSchwarzOperator (std::shared_ptr<const matrix_type> A, const communication_type& com)
98  : _A_(A), communication(com), buildcomm(true)
99  {}
100 
102  virtual void apply (const X& x, Y& y) const
103  {
104  y = 0;
105  novlp_op_apply(x,y,1);
106  communication.addOwnerCopyToOwnerCopy(y,y);
107  }
108 
110  virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
111  {
112  // only apply communication to alpha*A*x to make it consistent,
113  // y already has to be consistent.
114  Y y1(y);
115  y = 0;
116  novlp_op_apply(x,y,alpha);
117  communication.addOwnerCopyToOwnerCopy(y,y);
118  y += y1;
119  }
120 
122  virtual const matrix_type& getmat () const
123  {
124  return *_A_;
125  }
126 
127  void novlp_op_apply (const X& x, Y& y, field_type alpha) const
128  {
129  //get index sets
130  const PIS& pis=communication.indexSet();
131  const RI& ri = communication.remoteIndices();
132 
133  // at the beginning make a multimap "bordercontribution".
134  // process has i and j as border dofs but is not the owner
135  // => only contribute to Ax if i,j is in bordercontribution
136  if (buildcomm == true) {
137 
138  // set up mask vector
139  if (mask.size()!=static_cast<typename std::vector<double>::size_type>(x.size())) {
140  mask.resize(x.size());
141  for (typename std::vector<double>::size_type i=0; i<mask.size(); i++)
142  mask[i] = 1;
143  for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
144  if (i->local().attribute()==OwnerOverlapCopyAttributeSet::copy)
145  mask[i->local().local()] = 0;
146  else if (i->local().attribute()==OwnerOverlapCopyAttributeSet::overlap)
147  mask[i->local().local()] = 2;
148  }
149 
150  for (MM::iterator iter = bordercontribution.begin();
151  iter != bordercontribution.end(); ++iter)
152  bordercontribution.erase(iter);
153  std::map<int,int> owner; //key: local index i, value: process, that owns i
154  RIMap rimap;
155 
156  // for each local index make multimap rimap:
157  // key: local index i, data: pair of process that knows i and pointer to RI entry
158  for (RowIterator i = _A_->begin(); i != _A_->end(); ++i)
159  if (mask[i.index()] == 0)
160  for (RIIterator remote = ri.begin(); remote != ri.end(); ++remote) {
161  RIL& ril = *(remote->second.first);
162  for (RILIterator rindex = ril.begin(); rindex != ril.end(); ++rindex)
163  if (rindex->attribute() != OwnerOverlapCopyAttributeSet::overlap)
164  if (rindex->localIndexPair().local().local() == i.index()) {
165  rimap.insert
166  (std::make_pair(i.index(),
167  std::pair<int,RILIterator>(remote->first, rindex)));
168  if(rindex->attribute()==OwnerOverlapCopyAttributeSet::owner)
169  owner.insert(std::make_pair(i.index(),remote->first));
170  }
171  }
172 
173  int iowner = 0;
174  for (RowIterator i = _A_->begin(); i != _A_->end(); ++i) {
175  if (mask[i.index()] == 0) {
176  std::map<int,int>::iterator it = owner.find(i.index());
177  iowner = it->second;
178  std::pair<RIMapit, RIMapit> foundiit = rimap.equal_range(i.index());
179  for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()].end(); ++j) {
180  if (mask[j.index()] == 0) {
181  bool flag = true;
182  for (RIMapit foundi = foundiit.first; foundi != foundiit.second; ++foundi) {
183  std::pair<RIMapit, RIMapit> foundjit = rimap.equal_range(j.index());
184  for (RIMapit foundj = foundjit.first; foundj != foundjit.second; ++foundj)
185  if (foundj->second.first == foundi->second.first)
186  if (foundj->second.second->attribute() == OwnerOverlapCopyAttributeSet::owner
187  || foundj->second.first == iowner
188  || foundj->second.first < communication.communicator().rank()) {
189  flag = false;
190  continue;
191  }
192  if (flag == false)
193  continue;
194  }
195  // donĀ“t contribute to Ax if
196  // 1. the owner of j has i as interior/border dof
197  // 2. iowner has j as interior/border dof
198  // 3. there is another process with smaller rank that has i and j
199  // as interor/border dofs
200  // if the owner of j does not have i as interior/border dof,
201  // it will not be taken into account
202  if (flag==true)
203  bordercontribution.insert(std::pair<int,int>(i.index(),j.index()));
204  }
205  }
206  }
207  }
208  buildcomm = false;
209  }
210 
211  //compute alpha*A*x nonoverlapping case
212  for (RowIterator i = _A_->begin(); i != _A_->end(); ++i) {
213  if (mask[i.index()] == 0) {
214  //dof doesn't belong to process but is border (not ghost)
215  for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()].end(); ++j) {
216  if (mask[j.index()] == 1) //j is owner => then sum entries
217  Impl::asMatrix(*j).usmv(alpha,x[j.index()],y[i.index()]);
218  else if (mask[j.index()] == 0) {
219  std::pair<MM::iterator, MM::iterator> itp =
220  bordercontribution.equal_range(i.index());
221  for (MM::iterator it = itp.first; it != itp.second; ++it)
222  if ((*it).second == (int)j.index())
223  Impl::asMatrix(*j).usmv(alpha,x[j.index()],y[i.index()]);
224  }
225  }
226  }
227  else if (mask[i.index()] == 1) {
228  for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()].end(); ++j)
229  if (mask[j.index()] != 2)
230  Impl::asMatrix(*j).usmv(alpha,x[j.index()],y[i.index()]);
231  }
232  }
233  }
234 
237  {
239  }
240 
241  private:
242  std::shared_ptr<const matrix_type> _A_;
243  const communication_type& communication;
244  mutable bool buildcomm;
245  mutable std::vector<double> mask;
246  mutable std::multimap<int,int> bordercontribution;
247  };
248 
251  namespace Amg
252  {
253  template<class T> class ConstructionTraits;
254  }
255 
270  template<class C, class P>
272  : public Preconditioner<typename P::domain_type,typename P::range_type> {
274  using X = typename P::domain_type;
275  using Y = typename P::range_type;
276  public:
278  typedef typename P::domain_type domain_type;
280  typedef typename P::range_type range_type;
283 
299  : _preconditioner(stackobject_to_shared_ptr(p)), _communication(c)
300  { }
301 
309  NonoverlappingBlockPreconditioner (const std::shared_ptr<P>& p, const communication_type& c)
310  : _preconditioner(p), _communication(c)
311  { }
312 
318  virtual void pre (domain_type& x, range_type& b)
319  {
320  _preconditioner->pre(x,b);
321  }
322 
328  virtual void apply (domain_type& v, const range_type& d)
329  {
330  // block preconditioner equivalent to WrappedPreconditioner from
331  // pdelab/backend/ovlpistsolverbackend.hh,
332  // but not to BlockPreconditioner from schwarz.hh
333  _preconditioner->apply(v,d);
334  _communication.addOwnerCopyToOwnerCopy(v,v);
335  }
336 
337  template<bool forward>
338  void apply (X& v, const Y& d)
339  {
340  _preconditioner->template apply<forward>(v,d);
341  _communication.addOwnerCopyToOwnerCopy(v,v);
342  }
343 
349  virtual void post (domain_type& x)
350  {
351  _preconditioner->post(x);
352  }
353 
356  {
358  }
359 
360  private:
362  std::shared_ptr<P> _preconditioner;
363 
365  const communication_type& _communication;
366  };
367 
370 } // end namespace
371 
372 #endif
Implementation of the BCRSMatrix class.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Traits class for generically constructing non default constructable types.
Definition: construction.hh:38
A linear operator exporting itself in matrix form.
Definition: operators.hh:107
Nonoverlapping parallel preconditioner.
Definition: novlpschwarz.hh:272
NonoverlappingBlockPreconditioner(P &p, const communication_type &c)
Constructor.
Definition: novlpschwarz.hh:298
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: novlpschwarz.hh:355
virtual void apply(domain_type &v, const range_type &d)
Apply the preconditioner.
Definition: novlpschwarz.hh:328
P::range_type range_type
The range type of the preconditioner.
Definition: novlpschwarz.hh:280
NonoverlappingBlockPreconditioner(const std::shared_ptr< P > &p, const communication_type &c)
Constructor.
Definition: novlpschwarz.hh:309
virtual void post(domain_type &x)
Clean up.
Definition: novlpschwarz.hh:349
C communication_type
The type of the communication object.
Definition: novlpschwarz.hh:282
virtual void pre(domain_type &x, range_type &b)
Prepare the preconditioner.
Definition: novlpschwarz.hh:318
void apply(X &v, const Y &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition: novlpschwarz.hh:338
P::domain_type domain_type
The domain type of the preconditioner.
Definition: novlpschwarz.hh:278
A nonoverlapping operator with communication object.
Definition: novlpschwarz.hh:62
C communication_type
The type of the communication object.
Definition: novlpschwarz.hh:73
virtual void apply(const X &x, Y &y) const
apply operator to x:
Definition: novlpschwarz.hh:102
X domain_type
The type of the domain.
Definition: novlpschwarz.hh:67
virtual SolverCategory::Category category() const
Category of the linear operator (see SolverCategory::Category)
Definition: novlpschwarz.hh:236
Y range_type
The type of the range.
Definition: novlpschwarz.hh:69
M matrix_type
The type of the matrix we operate on.
Definition: novlpschwarz.hh:65
virtual const matrix_type & getmat() const
get matrix via *
Definition: novlpschwarz.hh:122
virtual void applyscaleadd(field_type alpha, const X &x, Y &y) const
apply operator to x, scale and add:
Definition: novlpschwarz.hh:110
X::field_type field_type
The field type of the range.
Definition: novlpschwarz.hh:71
NonoverlappingSchwarzOperator(const matrix_type &A, const communication_type &com)
constructor: just store a reference to a matrix.
Definition: novlpschwarz.hh:93
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:30
Some generic functions for pretty printing vectors and matrices.
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
Dune namespace.
Definition: alignedallocator.hh:14
shared_ptr< T > stackobject_to_shared_ptr(T &t)
Create a shared_ptr for a stack-allocated object.
Definition: shared_ptr.hh:75
Define general, extensible interface for operators. The available implementation wraps a matrix.
Classes providing communication interfaces for overlapping Schwarz methods.
Define general preconditioner interface.
Define base class for scalar product and norm.
Implementations of the inverse operator interface.
Category
Definition: solvercategory.hh:21
@ nonoverlapping
Category for non-overlapping solvers.
Definition: solvercategory.hh:25
A simple timing class.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.80.0 (May 12, 22:29, 2024)