dune-istl  2.7.0
novlpschwarz.hh
Go to the documentation of this file.
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 "io.hh"
16 #include "bvector.hh"
17 #include "vbvector.hh"
18 #include "bcrsmatrix.hh"
19 #include "io.hh"
20 #include "gsetc.hh"
21 #include "ilu.hh"
22 #include "operators.hh"
23 #include "solvers.hh"
24 #include "preconditioners.hh"
25 #include "scalarproducts.hh"
26 #include "owneroverlapcopy.hh"
27 
28 namespace Dune {
29 
58  template<class M, class X, class Y, class C>
60  {
61  public:
63  typedef M matrix_type;
65  typedef X domain_type;
67  typedef Y range_type;
69  typedef typename X::field_type field_type;
71  typedef C communication_type;
72 
73  typedef typename C::PIS PIS;
74  typedef typename C::RI RI;
75  typedef typename RI::RemoteIndexList RIL;
76  typedef typename RI::const_iterator RIIterator;
77  typedef typename RIL::const_iterator RILIterator;
78  typedef typename M::ConstColIterator ColIterator;
79  typedef typename M::ConstRowIterator RowIterator;
80  typedef std::multimap<int,int> MM;
81  typedef std::multimap<int,std::pair<int,RILIterator> > RIMap;
82  typedef typename RIMap::iterator RIMapit;
83 
92  : _A_(stackobject_to_shared_ptr(A)), communication(com), buildcomm(true)
93  {}
94 
95  NonoverlappingSchwarzOperator (std::shared_ptr<const matrix_type> A, const communication_type& com)
96  : _A_(A), communication(com), buildcomm(true)
97  {}
98 
100  virtual void apply (const X& x, Y& y) const
101  {
102  y = 0;
103  novlp_op_apply(x,y,1);
104  communication.addOwnerCopyToOwnerCopy(y,y);
105  }
106 
108  virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
109  {
110  // only apply communication to alpha*A*x to make it consistent,
111  // y already has to be consistent.
112  Y y1(y);
113  y = 0;
114  novlp_op_apply(x,y,alpha);
115  communication.addOwnerCopyToOwnerCopy(y,y);
116  y += y1;
117  }
118 
120  virtual const matrix_type& getmat () const
121  {
122  return *_A_;
123  }
124 
125  void novlp_op_apply (const X& x, Y& y, field_type alpha) const
126  {
127  //get index sets
128  const PIS& pis=communication.indexSet();
129  const RI& ri = communication.remoteIndices();
130 
131  // at the beginning make a multimap "bordercontribution".
132  // process has i and j as border dofs but is not the owner
133  // => only contribute to Ax if i,j is in bordercontribution
134  if (buildcomm == true) {
135 
136  // set up mask vector
137  if (mask.size()!=static_cast<typename std::vector<double>::size_type>(x.size())) {
138  mask.resize(x.size());
139  for (typename std::vector<double>::size_type i=0; i<mask.size(); i++)
140  mask[i] = 1;
141  for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
142  if (i->local().attribute()==OwnerOverlapCopyAttributeSet::copy)
143  mask[i->local().local()] = 0;
144  else if (i->local().attribute()==OwnerOverlapCopyAttributeSet::overlap)
145  mask[i->local().local()] = 2;
146  }
147 
148  for (MM::iterator iter = bordercontribution.begin();
149  iter != bordercontribution.end(); ++iter)
150  bordercontribution.erase(iter);
151  std::map<int,int> owner; //key: local index i, value: process, that owns i
152  RIMap rimap;
153 
154  // for each local index make multimap rimap:
155  // key: local index i, data: pair of process that knows i and pointer to RI entry
156  for (RowIterator i = _A_->begin(); i != _A_->end(); ++i)
157  if (mask[i.index()] == 0)
158  for (RIIterator remote = ri.begin(); remote != ri.end(); ++remote) {
159  RIL& ril = *(remote->second.first);
160  for (RILIterator rindex = ril.begin(); rindex != ril.end(); ++rindex)
161  if (rindex->attribute() != OwnerOverlapCopyAttributeSet::overlap)
162  if (rindex->localIndexPair().local().local() == i.index()) {
163  rimap.insert
164  (std::make_pair(i.index(),
165  std::pair<int,RILIterator>(remote->first, rindex)));
166  if(rindex->attribute()==OwnerOverlapCopyAttributeSet::owner)
167  owner.insert(std::make_pair(i.index(),remote->first));
168  }
169  }
170 
171  int iowner = 0;
172  for (RowIterator i = _A_->begin(); i != _A_->end(); ++i) {
173  if (mask[i.index()] == 0) {
174  std::map<int,int>::iterator it = owner.find(i.index());
175  iowner = it->second;
176  std::pair<RIMapit, RIMapit> foundiit = rimap.equal_range(i.index());
177  for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()].end(); ++j) {
178  if (mask[j.index()] == 0) {
179  bool flag = true;
180  for (RIMapit foundi = foundiit.first; foundi != foundiit.second; ++foundi) {
181  std::pair<RIMapit, RIMapit> foundjit = rimap.equal_range(j.index());
182  for (RIMapit foundj = foundjit.first; foundj != foundjit.second; ++foundj)
183  if (foundj->second.first == foundi->second.first)
184  if (foundj->second.second->attribute() == OwnerOverlapCopyAttributeSet::owner
185  || foundj->second.first == iowner
186  || foundj->second.first < communication.communicator().rank()) {
187  flag = false;
188  continue;
189  }
190  if (flag == false)
191  continue;
192  }
193  // donĀ“t contribute to Ax if
194  // 1. the owner of j has i as interior/border dof
195  // 2. iowner has j as interior/border dof
196  // 3. there is another process with smaller rank that has i and j
197  // as interor/border dofs
198  // if the owner of j does not have i as interior/border dof,
199  // it will not be taken into account
200  if (flag==true)
201  bordercontribution.insert(std::pair<int,int>(i.index(),j.index()));
202  }
203  }
204  }
205  }
206  buildcomm = false;
207  }
208 
209  //compute alpha*A*x nonoverlapping case
210  for (RowIterator i = _A_->begin(); i != _A_->end(); ++i) {
211  if (mask[i.index()] == 0) {
212  //dof doesn't belong to process but is border (not ghost)
213  for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()].end(); ++j) {
214  if (mask[j.index()] == 1) //j is owner => then sum entries
215  (*j).usmv(alpha,x[j.index()],y[i.index()]);
216  else if (mask[j.index()] == 0) {
217  std::pair<MM::iterator, MM::iterator> itp =
218  bordercontribution.equal_range(i.index());
219  for (MM::iterator it = itp.first; it != itp.second; ++it)
220  if ((*it).second == (int)j.index())
221  (*j).usmv(alpha,x[j.index()],y[i.index()]);
222  }
223  }
224  }
225  else if (mask[i.index()] == 1) {
226  for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()].end(); ++j)
227  if (mask[j.index()] != 2)
228  (*j).usmv(alpha,x[j.index()],y[i.index()]);
229  }
230  }
231  }
232 
235  {
237  }
238 
239  private:
240  std::shared_ptr<const matrix_type> _A_;
241  const communication_type& communication;
242  mutable bool buildcomm;
243  mutable std::vector<double> mask;
244  mutable std::multimap<int,int> bordercontribution;
245  };
246 
249  namespace Amg
250  {
251  template<class T> class ConstructionTraits;
252  }
253 
268  template<class C, class P>
270  : public Preconditioner<typename P::domain_type,typename P::range_type> {
272  using X = typename P::domain_type;
273  using Y = typename P::range_type;
274  public:
276  typedef typename P::domain_type domain_type;
278  typedef typename P::range_type range_type;
281 
297  : _preconditioner(stackobject_to_shared_ptr(p)), _communication(c)
298  { }
299 
307  NonoverlappingBlockPreconditioner (const std::shared_ptr<P>& p, const communication_type& c)
308  : _preconditioner(p), _communication(c)
309  { }
310 
316  virtual void pre (domain_type& x, range_type& b)
317  {
318  _preconditioner->pre(x,b);
319  }
320 
326  virtual void apply (domain_type& v, const range_type& d)
327  {
328  // block preconditioner equivalent to WrappedPreconditioner from
329  // pdelab/backend/ovlpistsolverbackend.hh,
330  // but not to BlockPreconditioner from schwarz.hh
331  _preconditioner->apply(v,d);
332  _communication.addOwnerCopyToOwnerCopy(v,v);
333  }
334 
335  template<bool forward>
336  void apply (X& v, const Y& d)
337  {
338  _preconditioner->template apply<forward>(v,d);
339  _communication.addOwnerCopyToOwnerCopy(v,v);
340  }
341 
347  virtual void post (domain_type& x)
348  {
349  _preconditioner->post(x);
350  }
351 
354  {
356  }
357 
358  private:
360  std::shared_ptr<P> _preconditioner;
361 
363  const communication_type& _communication;
364  };
365 
368 } // end namespace
369 
370 #endif
gsetc.hh
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
Dune::NonoverlappingSchwarzOperator::range_type
Y range_type
The type of the range.
Definition: novlpschwarz.hh:67
Dune::NonoverlappingSchwarzOperator::PIS
C::PIS PIS
Definition: novlpschwarz.hh:73
scalarproducts.hh
Define base class for scalar product and norm.
Dune::NonoverlappingBlockPreconditioner::pre
virtual void pre(domain_type &x, range_type &b)
Prepare the preconditioner.
Definition: novlpschwarz.hh:316
Dune::AssembledLinearOperator
A linear operator exporting itself in matrix form.
Definition: operators.hh:107
Dune::Amg::ConstructionTraits
Traits class for generically constructing non default constructable types.
Definition: novlpschwarz.hh:251
Dune::NonoverlappingBlockPreconditioner::NonoverlappingBlockPreconditioner
NonoverlappingBlockPreconditioner(P &p, const communication_type &c)
Constructor.
Definition: novlpschwarz.hh:296
Dune::NonoverlappingSchwarzOperator::getmat
virtual const matrix_type & getmat() const
get matrix via *
Definition: novlpschwarz.hh:120
Dune::NonoverlappingSchwarzOperator::MM
std::multimap< int, int > MM
Definition: novlpschwarz.hh:80
Dune::OwnerOverlapCopyAttributeSet::copy
@ copy
Definition: owneroverlapcopy.hh:59
Dune::SolverCategory::nonoverlapping
@ nonoverlapping
Category for non-overlapping solvers.
Definition: solvercategory.hh:25
Dune::NonoverlappingSchwarzOperator::domain_type
X domain_type
The type of the domain.
Definition: novlpschwarz.hh:65
Dune::NonoverlappingSchwarzOperator::RI
C::RI RI
Definition: novlpschwarz.hh:74
Dune::NonoverlappingSchwarzOperator::RowIterator
M::ConstRowIterator RowIterator
Definition: novlpschwarz.hh:79
Dune::NonoverlappingSchwarzOperator::NonoverlappingSchwarzOperator
NonoverlappingSchwarzOperator(const matrix_type &A, const communication_type &com)
constructor: just store a reference to a matrix.
Definition: novlpschwarz.hh:91
Dune::NonoverlappingSchwarzOperator::NonoverlappingSchwarzOperator
NonoverlappingSchwarzOperator(std::shared_ptr< const matrix_type > A, const communication_type &com)
Definition: novlpschwarz.hh:95
vbvector.hh
???
Dune::NonoverlappingSchwarzOperator::novlp_op_apply
void novlp_op_apply(const X &x, Y &y, field_type alpha) const
Definition: novlpschwarz.hh:125
Dune::NonoverlappingSchwarzOperator::communication_type
C communication_type
The type of the communication object.
Definition: novlpschwarz.hh:71
bcrsmatrix.hh
Implementation of the BCRSMatrix class.
io.hh
Some generic functions for pretty printing vectors and matrices.
Dune::NonoverlappingSchwarzOperator::RIMap
std::multimap< int, std::pair< int, RILIterator > > RIMap
Definition: novlpschwarz.hh:81
solvers.hh
Implementations of the inverse operator interface.
Dune::NonoverlappingBlockPreconditioner
Nonoverlapping parallel preconditioner.
Definition: novlpschwarz.hh:269
Dune
Definition: allocator.hh:7
Dune::NonoverlappingSchwarzOperator::field_type
X::field_type field_type
The field type of the range.
Definition: novlpschwarz.hh:69
Dune::NonoverlappingBlockPreconditioner::category
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: novlpschwarz.hh:353
Dune::NonoverlappingBlockPreconditioner::apply
void apply(X &v, const Y &d)
Definition: novlpschwarz.hh:336
Dune::NonoverlappingSchwarzOperator
A nonoverlapping operator with communication object.
Definition: novlpschwarz.hh:59
Dune::NonoverlappingSchwarzOperator::applyscaleadd
virtual void applyscaleadd(field_type alpha, const X &x, Y &y) const
apply operator to x, scale and add:
Definition: novlpschwarz.hh:108
Dune::NonoverlappingSchwarzOperator::category
virtual SolverCategory::Category category() const
Category of the linear operator (see SolverCategory::Category)
Definition: novlpschwarz.hh:234
Dune::NonoverlappingSchwarzOperator::RIMapit
RIMap::iterator RIMapit
Definition: novlpschwarz.hh:82
Dune::NonoverlappingSchwarzOperator::ColIterator
M::ConstColIterator ColIterator
Definition: novlpschwarz.hh:78
ilu.hh
???
Dune::NonoverlappingBlockPreconditioner::communication_type
C communication_type
The type of the communication object.
Definition: novlpschwarz.hh:280
Dune::OwnerOverlapCopyAttributeSet::owner
@ owner
Definition: owneroverlapcopy.hh:59
Dune::NonoverlappingBlockPreconditioner::NonoverlappingBlockPreconditioner
NonoverlappingBlockPreconditioner(const std::shared_ptr< P > &p, const communication_type &c)
Constructor.
Definition: novlpschwarz.hh:307
Dune::Preconditioner
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:30
Dune::NonoverlappingSchwarzOperator::RIIterator
RI::const_iterator RIIterator
Definition: novlpschwarz.hh:76
owneroverlapcopy.hh
Classes providing communication interfaces for overlapping Schwarz methods.
Dune::NonoverlappingBlockPreconditioner::domain_type
P::domain_type domain_type
The domain type of the preconditioner.
Definition: novlpschwarz.hh:276
Dune::OwnerOverlapCopyAttributeSet::overlap
@ overlap
Definition: owneroverlapcopy.hh:59
Dune::NonoverlappingSchwarzOperator::RIL
RI::RemoteIndexList RIL
Definition: novlpschwarz.hh:75
preconditioners.hh
Define general preconditioner interface.
Dune::NonoverlappingSchwarzOperator::RILIterator
RIL::const_iterator RILIterator
Definition: novlpschwarz.hh:77
Dune::NonoverlappingBlockPreconditioner::post
virtual void post(domain_type &x)
Clean up.
Definition: novlpschwarz.hh:347
Dune::SolverCategory::Category
Category
Definition: solvercategory.hh:21
bvector.hh
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Dune::NonoverlappingBlockPreconditioner::range_type
P::range_type range_type
The range type of the preconditioner.
Definition: novlpschwarz.hh:278
Dune::NonoverlappingSchwarzOperator::apply
virtual void apply(const X &x, Y &y) const
apply operator to x:
Definition: novlpschwarz.hh:100
Dune::NonoverlappingBlockPreconditioner::apply
virtual void apply(domain_type &v, const range_type &d)
Apply the preconditioner.
Definition: novlpschwarz.hh:326
Dune::NonoverlappingSchwarzOperator::matrix_type
M matrix_type
The type of the matrix we operate on.
Definition: novlpschwarz.hh:63
operators.hh
Define general, extensible interface for operators. The available implementation wraps a matrix.