dune-istl 2.9.0
Loading...
Searching...
No Matches
novlpschwarz.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright (C) DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_ISTL_NOVLPSCHWARZ_HH
6#define DUNE_ISTL_NOVLPSCHWARZ_HH
7
8#include <iostream> // for input/output to shell
9#include <fstream> // for input/output to files
10#include <vector> // STL vector class
11#include <sstream>
12
13#include <cmath> // Yes, we do some math here
14
15#include <dune/common/timer.hh>
16
17#include <dune/common/hybridutilities.hh>
18
19#include "io.hh"
20#include "bvector.hh"
21#include "vbvector.hh"
22#include "bcrsmatrix.hh"
23#include "io.hh"
24#include "gsetc.hh"
25#include "ilu.hh"
26#include "operators.hh"
27#include "solvers.hh"
28#include "preconditioners.hh"
29#include "scalarproducts.hh"
30#include "owneroverlapcopy.hh"
31
32namespace Dune {
33
59 template<class M, class X, class Y, class C>
61 {
62 public:
64 typedef M matrix_type;
66 typedef X domain_type;
68 typedef Y range_type;
70 typedef typename X::field_type field_type;
73
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;
84
93 : _A_(stackobject_to_shared_ptr(A)), communication(com), buildcomm(true)
94 {}
95
96 NonoverlappingSchwarzOperator (std::shared_ptr<const matrix_type> A, const communication_type& com)
97 : _A_(A), communication(com), buildcomm(true)
98 {}
99
101 virtual void apply (const X& x, Y& y) const
102 {
103 y = 0;
104 novlp_op_apply(x,y,1);
105 communication.addOwnerCopyToOwnerCopy(y,y);
106 }
107
109 virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
110 {
111 // only apply communication to alpha*A*x to make it consistent,
112 // y already has to be consistent.
113 Y y1(y);
114 y = 0;
115 novlp_op_apply(x,y,alpha);
116 communication.addOwnerCopyToOwnerCopy(y,y);
117 y += y1;
118 }
119
121 virtual const matrix_type& getmat () const
122 {
123 return *_A_;
124 }
125
126 void novlp_op_apply (const X& x, Y& y, field_type alpha) const
127 {
128 //get index sets
129 const PIS& pis=communication.indexSet();
130 const RI& ri = communication.remoteIndices();
131
132 // at the beginning make a multimap "bordercontribution".
133 // process has i and j as border dofs but is not the owner
134 // => only contribute to Ax if i,j is in bordercontribution
135 if (buildcomm == true) {
136
137 // set up mask vector
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++)
141 mask[i] = 1;
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;
147 }
148
149 for (MM::iterator iter = bordercontribution.begin();
150 iter != bordercontribution.end(); ++iter)
151 bordercontribution.erase(iter);
152 std::map<int,int> owner; //key: local index i, value: process, that owns i
153 RIMap rimap;
154
155 // for each local index make multimap rimap:
156 // key: local index i, data: pair of process that knows i and pointer to RI entry
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()) {
164 rimap.insert
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));
169 }
170 }
171
172 int iowner = 0;
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());
176 iowner = it->second;
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) {
180 bool flag = true;
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()) {
188 flag = false;
189 continue;
190 }
191 if (flag == false)
192 continue;
193 }
194 // don“t contribute to Ax if
195 // 1. the owner of j has i as interior/border dof
196 // 2. iowner has j as interior/border dof
197 // 3. there is another process with smaller rank that has i and j
198 // as interor/border dofs
199 // if the owner of j does not have i as interior/border dof,
200 // it will not be taken into account
201 if (flag==true)
202 bordercontribution.insert(std::pair<int,int>(i.index(),j.index()));
203 }
204 }
205 }
206 }
207 buildcomm = false;
208 }
209
210 //compute alpha*A*x nonoverlapping case
211 for (RowIterator i = _A_->begin(); i != _A_->end(); ++i) {
212 if (mask[i.index()] == 0) {
213 //dof doesn't belong to process but is border (not ghost)
214 for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()].end(); ++j) {
215 if (mask[j.index()] == 1) //j is owner => then sum entries
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()]);
223 }
224 }
225 }
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()]);
230 }
231 }
232 }
233
236 {
238 }
239
242 {
243 return communication;
244 }
245 private:
246 std::shared_ptr<const matrix_type> _A_;
247 const communication_type& communication;
248 mutable bool buildcomm;
249 mutable std::vector<double> mask;
250 mutable std::multimap<int,int> bordercontribution;
251 };
252
255 namespace Amg
256 {
257 template<class T> struct ConstructionTraits;
258 }
259
274 template<class C, class P>
276 : public Preconditioner<typename P::domain_type,typename P::range_type> {
277 friend struct Amg::ConstructionTraits<NonoverlappingBlockPreconditioner<C,P> >;
278 using X = typename P::domain_type;
279 using Y = typename P::range_type;
280 public:
282 typedef typename P::domain_type domain_type;
284 typedef typename P::range_type range_type;
287
303 : _preconditioner(stackobject_to_shared_ptr(p)), _communication(c)
304 { }
305
313 NonoverlappingBlockPreconditioner (const std::shared_ptr<P>& p, const communication_type& c)
314 : _preconditioner(p), _communication(c)
315 { }
316
322 virtual void pre (domain_type& x, range_type& b)
323 {
324 _preconditioner->pre(x,b);
325 }
326
332 virtual void apply (domain_type& v, const range_type& d)
333 {
334 // block preconditioner equivalent to WrappedPreconditioner from
335 // pdelab/backend/ovlpistsolverbackend.hh,
336 // but not to BlockPreconditioner from schwarz.hh
337 _preconditioner->apply(v,d);
338 _communication.addOwnerCopyToOwnerCopy(v,v);
339 }
340
341 template<bool forward>
342 void apply (X& v, const Y& d)
343 {
344 _preconditioner->template apply<forward>(v,d);
345 _communication.addOwnerCopyToOwnerCopy(v,v);
346 }
347
353 virtual void post (domain_type& x)
354 {
355 _preconditioner->post(x);
356 }
357
360 {
362 }
363
364 private:
366 std::shared_ptr<P> _preconditioner;
367
369 const communication_type& _communication;
370 };
371
374} // end namespace
375
376#endif
Define base class for scalar product and norm.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
The incomplete LU factorization kernels.
Implementation of the BCRSMatrix class.
Classes providing communication interfaces for overlapping Schwarz methods.
Define general preconditioner interface.
Implementations of the inverse operator interface.
Some generic functions for pretty printing vectors and matrices.
Define general, extensible interface for operators. The available implementation wraps a matrix.
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