5#ifndef DUNE_ISTL_NOVLPSCHWARZ_HH
6#define DUNE_ISTL_NOVLPSCHWARZ_HH
15#include <dune/common/timer.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;
81 typedef std::multimap<int,int>
MM;
82 typedef std::multimap<int,std::pair<int,RILIterator> >
RIMap;
93 : _A_(stackobject_to_shared_ptr(A)), communication(com), buildcomm(true)
97 : _A_(A), communication(com), buildcomm(true)
101 virtual void apply (
const X& x, Y& y)
const
105 communication.addOwnerCopyToOwnerCopy(y,y);
116 communication.addOwnerCopyToOwnerCopy(y,y);
129 const PIS& pis=communication.indexSet();
130 const RI& ri = communication.remoteIndices();
135 if (buildcomm ==
true) {
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)
144 mask[i->local().local()] = 0;
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)
163 if (rindex->localIndexPair().local().local() == i.index()) {
165 (std::make_pair(i.index(),
166 std::pair<int,RILIterator>(remote->first, rindex)));
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)
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;
303 : _preconditioner(stackobject_to_shared_ptr(p)), _communication(c)
314 : _preconditioner(p), _communication(c)
324 _preconditioner->pre(x,b);
337 _preconditioner->apply(v,d);
338 _communication.addOwnerCopyToOwnerCopy(v,v);
341 template<
bool forward>
344 _preconditioner->template apply<forward>(v,d);
345 _communication.addOwnerCopyToOwnerCopy(v,v);
355 _preconditioner->post(x);
366 std::shared_ptr<P> _preconditioner;
The incomplete LU factorization kernels.
Classes providing communication interfaces for overlapping Schwarz methods.
Implementations of the inverse operator interface.
Some generic functions for pretty printing vectors and matrices.
Implementation of the BCRSMatrix class.
Define general, extensible interface for operators. The available implementation wraps a matrix.
Define general preconditioner interface.
Define base class for scalar product and norm.
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Definition allocator.hh:11
A nonoverlapping operator with communication object.
Definition novlpschwarz.hh:61
C::PIS PIS
Definition novlpschwarz.hh:74
C communication_type
The type of the communication object.
Definition novlpschwarz.hh:72
std::multimap< int, std::pair< int, RILIterator > > RIMap
Definition novlpschwarz.hh:82
C::RI RI
Definition novlpschwarz.hh:75
void novlp_op_apply(const X &x, Y &y, field_type alpha) const
Definition novlpschwarz.hh:126
NonoverlappingSchwarzOperator(std::shared_ptr< const matrix_type > A, const communication_type &com)
Definition novlpschwarz.hh:96
M::ConstColIterator ColIterator
Definition novlpschwarz.hh:79
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
RIMap::iterator RIMapit
Definition novlpschwarz.hh:83
virtual const matrix_type & getmat() const
get matrix via *
Definition novlpschwarz.hh:121
RIL::const_iterator RILIterator
Definition novlpschwarz.hh:78
std::multimap< int, int > MM
Definition novlpschwarz.hh:81
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
M::ConstRowIterator RowIterator
Definition novlpschwarz.hh:80
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
RI::RemoteIndexList RIL
Definition novlpschwarz.hh:76
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
RI::const_iterator RIIterator
Definition novlpschwarz.hh:77
Nonoverlapping parallel preconditioner.
Definition novlpschwarz.hh:276
NonoverlappingBlockPreconditioner(P &p, const communication_type &c)
Constructor.
Definition novlpschwarz.hh:302
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition novlpschwarz.hh:359
virtual void apply(domain_type &v, const range_type &d)
Apply the preconditioner.
Definition novlpschwarz.hh:332
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:313
virtual void post(domain_type &x)
Clean up.
Definition novlpschwarz.hh:353
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:322
void apply(X &v, const Y &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition novlpschwarz.hh:342
P::domain_type domain_type
The domain type of the preconditioner.
Definition novlpschwarz.hh:282
A linear operator exporting itself in matrix form.
Definition operators.hh:109
@ owner
Definition owneroverlapcopy.hh:61
@ copy
Definition owneroverlapcopy.hh:61
@ overlap
Definition owneroverlapcopy.hh:61
Base class for matrix free definition of preconditioners.
Definition preconditioner.hh:32
Category
Definition solvercategory.hh:23
@ nonoverlapping
Category for non-overlapping solvers.
Definition solvercategory.hh:27