5#ifndef DUNE_ISTL_NOVLPSCHWARZ_HH
6#define DUNE_ISTL_NOVLPSCHWARZ_HH
17#include <dune/common/hybridutilities.hh>
59 template<
class M,
class X,
class Y,
class C>
74 typedef typename C::PIS PIS;
75 typedef typename C::RI RI;
76 typedef typename RI::RemoteIndexList RIL;
77 typedef typename RI::const_iterator RIIterator;
78 typedef typename RIL::const_iterator RILIterator;
79 typedef typename M::ConstColIterator ColIterator;
80 typedef typename M::ConstRowIterator RowIterator;
81 typedef std::multimap<int,int> MM;
82 typedef std::multimap<int,std::pair<int,RILIterator> > RIMap;
83 typedef typename RIMap::iterator RIMapit;
97 : _A_(A), communication(com), buildcomm(true)
101 virtual void apply (
const X& x, Y& y)
const
104 novlp_op_apply(x,y,1);
105 communication.addOwnerCopyToOwnerCopy(y,y);
115 novlp_op_apply(x,y,alpha);
116 communication.addOwnerCopyToOwnerCopy(y,y);
126 void novlp_op_apply (
const X& x, Y& y,
field_type alpha)
const
129 const PIS& pis=communication.indexSet();
130 const RI& ri = communication.remoteIndices();
138 if (mask.size()!=
static_cast<typename std::vector<double>::size_type
>(x.size())) {
139 mask.resize(x.size());
140 for (
typename std::vector<double>::size_type i=0; i<mask.size(); i++)
142 for (
typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
143 if (i->local().attribute()==OwnerOverlapCopyAttributeSet::copy)
144 mask[i->local().local()] = 0;
145 else if (i->local().attribute()==OwnerOverlapCopyAttributeSet::overlap)
146 mask[i->local().local()] = 2;
149 for (MM::iterator iter = bordercontribution.begin();
150 iter != bordercontribution.end(); ++iter)
151 bordercontribution.erase(iter);
152 std::map<int,int> owner;
157 for (RowIterator i = _A_->begin(); i != _A_->end(); ++i)
158 if (mask[i.index()] == 0)
159 for (RIIterator remote = ri.begin(); remote != ri.end(); ++remote) {
160 RIL& ril = *(remote->second.first);
161 for (RILIterator rindex = ril.begin(); rindex != ril.end(); ++rindex)
162 if (rindex->attribute() != OwnerOverlapCopyAttributeSet::overlap)
163 if (rindex->localIndexPair().local().local() == i.index()) {
165 (std::make_pair(i.index(),
166 std::pair<int,RILIterator>(remote->first, rindex)));
167 if(rindex->attribute()==OwnerOverlapCopyAttributeSet::owner)
168 owner.insert(std::make_pair(i.index(),remote->first));
173 for (RowIterator i = _A_->begin(); i != _A_->end(); ++i) {
174 if (mask[i.index()] == 0) {
175 std::map<int,int>::iterator it = owner.find(i.index());
177 std::pair<RIMapit, RIMapit> foundiit = rimap.equal_range(i.index());
178 for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()].end(); ++j) {
179 if (mask[j.index()] == 0) {
181 for (RIMapit foundi = foundiit.first; foundi != foundiit.second; ++foundi) {
182 std::pair<RIMapit, RIMapit> foundjit = rimap.equal_range(j.index());
183 for (RIMapit foundj = foundjit.first; foundj != foundjit.second; ++foundj)
184 if (foundj->second.first == foundi->second.first)
185 if (foundj->second.second->attribute() == OwnerOverlapCopyAttributeSet::owner
186 || foundj->second.first == iowner
187 || foundj->second.first < communication.communicator().rank()) {
202 bordercontribution.insert(std::pair<int,int>(i.index(),j.index()));
211 for (RowIterator i = _A_->begin(); i != _A_->end(); ++i) {
212 if (mask[i.index()] == 0) {
214 for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()].end(); ++j) {
215 if (mask[j.index()] == 1)
216 Impl::asMatrix(*j).usmv(alpha,x[j.index()],y[i.index()]);
217 else if (mask[j.index()] == 0) {
218 std::pair<MM::iterator, MM::iterator> itp =
219 bordercontribution.equal_range(i.index());
220 for (MM::iterator it = itp.first; it != itp.second; ++it)
221 if ((*it).second == (
int)j.index())
222 Impl::asMatrix(*j).usmv(alpha,x[j.index()],y[i.index()]);
226 else if (mask[i.index()] == 1) {
227 for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()].end(); ++j)
228 if (mask[j.index()] != 2)
229 Impl::asMatrix(*j).usmv(alpha,x[j.index()],y[i.index()]);
243 return communication;
246 std::shared_ptr<const matrix_type> _A_;
248 mutable bool buildcomm;
249 mutable std::vector<double> mask;
250 mutable std::multimap<int,int> bordercontribution;
257 template<
class T>
struct ConstructionTraits;
274 template<
class C,
class P>
276 :
public Preconditioner<typename P::domain_type,typename P::range_type> {
278 using X =
typename P::domain_type;
279 using Y =
typename P::range_type;
307 : _preconditioner(p), _communication(c)
317 _preconditioner->pre(x,b);
330 _preconditioner->apply(v,d);
331 _communication.addOwnerCopyToOwnerCopy(v,v);
334 template<
bool forward>
337 _preconditioner->template apply<forward>(v,d);
338 _communication.addOwnerCopyToOwnerCopy(v,v);
348 _preconditioner->post(x);
359 std::shared_ptr<P> _preconditioner;
Implementation of the BCRSMatrix class.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
A linear operator exporting itself in matrix form.
Definition: operators.hh:111
Nonoverlapping parallel preconditioner.
Definition: novlpschwarz.hh:276
NonoverlappingBlockPreconditioner(P &p, const communication_type &c)
Constructor.
Definition: novlpschwarz.hh:295
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: novlpschwarz.hh:352
virtual void apply(domain_type &v, const range_type &d)
Apply the preconditioner.
Definition: novlpschwarz.hh:325
P::range_type range_type
The range type of the preconditioner.
Definition: novlpschwarz.hh:284
NonoverlappingBlockPreconditioner(const std::shared_ptr< P > &p, const communication_type &c)
Constructor.
Definition: novlpschwarz.hh:306
virtual void post(domain_type &x)
Clean up.
Definition: novlpschwarz.hh:346
C communication_type
The type of the communication object.
Definition: novlpschwarz.hh:286
virtual void pre(domain_type &x, range_type &b)
Prepare the preconditioner.
Definition: novlpschwarz.hh:315
void apply(X &v, const Y &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition: novlpschwarz.hh:335
P::domain_type domain_type
The domain type of the preconditioner.
Definition: novlpschwarz.hh:282
A nonoverlapping operator with communication object.
Definition: novlpschwarz.hh:61
C communication_type
The type of the communication object.
Definition: novlpschwarz.hh:72
virtual void apply(const X &x, Y &y) const
apply operator to x:
Definition: novlpschwarz.hh:101
X domain_type
The type of the domain.
Definition: novlpschwarz.hh:66
virtual SolverCategory::Category category() const
Category of the linear operator (see SolverCategory::Category)
Definition: novlpschwarz.hh:235
virtual const matrix_type & getmat() const
get matrix via *
Definition: novlpschwarz.hh:121
Y range_type
The type of the range.
Definition: novlpschwarz.hh:68
M matrix_type
The type of the matrix we operate on.
Definition: novlpschwarz.hh:64
const communication_type & getCommunication() const
Get the object responsible for communication.
Definition: novlpschwarz.hh:241
virtual void applyscaleadd(field_type alpha, const X &x, Y &y) const
apply operator to x, scale and add:
Definition: novlpschwarz.hh:109
X::field_type field_type
The field type of the range.
Definition: novlpschwarz.hh:70
NonoverlappingSchwarzOperator(const matrix_type &A, const communication_type &com)
constructor: just store a reference to a matrix.
Definition: novlpschwarz.hh:92
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:33
Some generic functions for pretty printing vectors and matrices.
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
The incomplete LU factorization kernels.
Dune namespace.
Definition: alignedallocator.hh:13
std::shared_ptr< T > stackobject_to_shared_ptr(T &t)
Create a shared_ptr for a stack-allocated object.
Definition: shared_ptr.hh:72
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.
Traits class for generically constructing non default constructable types.
Definition: construction.hh:39
Category
Definition: solvercategory.hh:23
@ nonoverlapping
Category for non-overlapping solvers.
Definition: solvercategory.hh:27