3#ifndef DUNE_ISTL_NOVLPSCHWARZ_HH
4#define DUNE_ISTL_NOVLPSCHWARZ_HH
15#include <dune/common/hybridutilities.hh>
60 template<
class M,
class X,
class Y,
class C>
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;
98 : _A_(A), communication(com), buildcomm(true)
102 virtual void apply (
const X& x, Y& y)
const
105 novlp_op_apply(x,y,1);
106 communication.addOwnerCopyToOwnerCopy(y,y);
116 novlp_op_apply(x,y,alpha);
117 communication.addOwnerCopyToOwnerCopy(y,y);
127 void novlp_op_apply (
const X& x, Y& y,
field_type alpha)
const
130 const PIS& pis=communication.indexSet();
131 const RI& ri = communication.remoteIndices();
136 if (buildcomm ==
true) {
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++)
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;
150 for (MM::iterator iter = bordercontribution.begin();
151 iter != bordercontribution.end(); ++iter)
152 bordercontribution.erase(iter);
153 std::map<int,int> owner;
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()) {
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));
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());
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) {
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()) {
203 bordercontribution.insert(std::pair<int,int>(i.index(),j.index()));
212 for (RowIterator i = _A_->begin(); i != _A_->end(); ++i) {
213 if (mask[i.index()] == 0) {
215 for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()].end(); ++j) {
216 if (mask[j.index()] == 1)
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()]);
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()]);
244 return communication;
247 std::shared_ptr<const matrix_type> _A_;
249 mutable bool buildcomm;
250 mutable std::vector<double> mask;
251 mutable std::multimap<int,int> bordercontribution;
258 template<
class T>
struct ConstructionTraits;
275 template<
class C,
class P>
277 :
public Preconditioner<typename P::domain_type,typename P::range_type> {
279 using X =
typename P::domain_type;
280 using Y =
typename P::range_type;
315 : _preconditioner(p), _communication(c)
325 _preconditioner->pre(x,b);
338 _preconditioner->apply(v,d);
339 _communication.addOwnerCopyToOwnerCopy(v,v);
342 template<
bool forward>
345 _preconditioner->template apply<forward>(v,d);
346 _communication.addOwnerCopyToOwnerCopy(v,v);
356 _preconditioner->post(x);
367 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:107
Nonoverlapping parallel preconditioner.
Definition: novlpschwarz.hh:277
NonoverlappingBlockPreconditioner(P &p, const communication_type &c)
Constructor.
Definition: novlpschwarz.hh:303
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: novlpschwarz.hh:360
virtual void apply(domain_type &v, const range_type &d)
Apply the preconditioner.
Definition: novlpschwarz.hh:333
P::range_type range_type
The range type of the preconditioner.
Definition: novlpschwarz.hh:285
NonoverlappingBlockPreconditioner(const std::shared_ptr< P > &p, const communication_type &c)
Constructor.
Definition: novlpschwarz.hh:314
virtual void post(domain_type &x)
Clean up.
Definition: novlpschwarz.hh:354
C communication_type
The type of the communication object.
Definition: novlpschwarz.hh:287
virtual void pre(domain_type &x, range_type &b)
Prepare the preconditioner.
Definition: novlpschwarz.hh:323
void apply(X &v, const Y &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition: novlpschwarz.hh:343
P::domain_type domain_type
The domain type of the preconditioner.
Definition: novlpschwarz.hh:283
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
virtual const matrix_type & getmat() const
get matrix via *
Definition: novlpschwarz.hh:122
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
const communication_type & getCommunication() const
Get the object responsible for communication.
Definition: novlpschwarz.hh:242
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.
The incomplete LU factorization kernels.
Dune namespace.
Definition: alignedallocator.hh:11
std::shared_ptr< T > stackobject_to_shared_ptr(T &t)
Create a shared_ptr for a stack-allocated object.
Definition: shared_ptr.hh:70
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:37
Category
Definition: solvercategory.hh:21
@ nonoverlapping
Category for non-overlapping solvers.
Definition: solvercategory.hh:25