Dune Core Modules (2.8.0)

novlpschwarz.hh
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 <dune/common/hybridutilities.hh>
16
17#include "io.hh"
18#include "bvector.hh"
19#include "vbvector.hh"
20#include "bcrsmatrix.hh"
21#include "io.hh"
22#include "gsetc.hh"
23#include "ilu.hh"
24#include "operators.hh"
25#include "solvers.hh"
26#include "preconditioners.hh"
27#include "scalarproducts.hh"
28#include "owneroverlapcopy.hh"
29
30namespace Dune {
31
60 template<class M, class X, class Y, class C>
62 {
63 public:
65 typedef M matrix_type;
67 typedef X domain_type;
69 typedef Y range_type;
71 typedef typename X::field_type field_type;
74
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;
85
94 : _A_(stackobject_to_shared_ptr(A)), communication(com), buildcomm(true)
95 {}
96
97 NonoverlappingSchwarzOperator (std::shared_ptr<const matrix_type> A, const communication_type& com)
98 : _A_(A), communication(com), buildcomm(true)
99 {}
100
102 virtual void apply (const X& x, Y& y) const
103 {
104 y = 0;
105 novlp_op_apply(x,y,1);
106 communication.addOwnerCopyToOwnerCopy(y,y);
107 }
108
110 virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
111 {
112 // only apply communication to alpha*A*x to make it consistent,
113 // y already has to be consistent.
114 Y y1(y);
115 y = 0;
116 novlp_op_apply(x,y,alpha);
117 communication.addOwnerCopyToOwnerCopy(y,y);
118 y += y1;
119 }
120
122 virtual const matrix_type& getmat () const
123 {
124 return *_A_;
125 }
126
127 void novlp_op_apply (const X& x, Y& y, field_type alpha) const
128 {
129 //get index sets
130 const PIS& pis=communication.indexSet();
131 const RI& ri = communication.remoteIndices();
132
133 // at the beginning make a multimap "bordercontribution".
134 // process has i and j as border dofs but is not the owner
135 // => only contribute to Ax if i,j is in bordercontribution
136 if (buildcomm == true) {
137
138 // set up mask vector
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++)
142 mask[i] = 1;
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;
148 }
149
150 for (MM::iterator iter = bordercontribution.begin();
151 iter != bordercontribution.end(); ++iter)
152 bordercontribution.erase(iter);
153 std::map<int,int> owner; //key: local index i, value: process, that owns i
154 RIMap rimap;
155
156 // for each local index make multimap rimap:
157 // key: local index i, data: pair of process that knows i and pointer to RI entry
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()) {
165 rimap.insert
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));
170 }
171 }
172
173 int iowner = 0;
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());
177 iowner = it->second;
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) {
181 bool flag = true;
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()) {
189 flag = false;
190 continue;
191 }
192 if (flag == false)
193 continue;
194 }
195 // donĀ“t contribute to Ax if
196 // 1. the owner of j has i as interior/border dof
197 // 2. iowner has j as interior/border dof
198 // 3. there is another process with smaller rank that has i and j
199 // as interor/border dofs
200 // if the owner of j does not have i as interior/border dof,
201 // it will not be taken into account
202 if (flag==true)
203 bordercontribution.insert(std::pair<int,int>(i.index(),j.index()));
204 }
205 }
206 }
207 }
208 buildcomm = false;
209 }
210
211 //compute alpha*A*x nonoverlapping case
212 for (RowIterator i = _A_->begin(); i != _A_->end(); ++i) {
213 if (mask[i.index()] == 0) {
214 //dof doesn't belong to process but is border (not ghost)
215 for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()].end(); ++j) {
216 if (mask[j.index()] == 1) //j is owner => then sum entries
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()]);
224 }
225 }
226 }
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()]);
231 }
232 }
233 }
234
237 {
239 }
240
243 {
244 return communication;
245 }
246 private:
247 std::shared_ptr<const matrix_type> _A_;
248 const communication_type& communication;
249 mutable bool buildcomm;
250 mutable std::vector<double> mask;
251 mutable std::multimap<int,int> bordercontribution;
252 };
253
256 namespace Amg
257 {
258 template<class T> struct ConstructionTraits;
259 }
260
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;
281 public:
283 typedef typename P::domain_type domain_type;
285 typedef typename P::range_type range_type;
288
304 : _preconditioner(stackobject_to_shared_ptr(p)), _communication(c)
305 { }
306
314 NonoverlappingBlockPreconditioner (const std::shared_ptr<P>& p, const communication_type& c)
315 : _preconditioner(p), _communication(c)
316 { }
317
323 virtual void pre (domain_type& x, range_type& b)
324 {
325 _preconditioner->pre(x,b);
326 }
327
333 virtual void apply (domain_type& v, const range_type& d)
334 {
335 // block preconditioner equivalent to WrappedPreconditioner from
336 // pdelab/backend/ovlpistsolverbackend.hh,
337 // but not to BlockPreconditioner from schwarz.hh
338 _preconditioner->apply(v,d);
339 _communication.addOwnerCopyToOwnerCopy(v,v);
340 }
341
342 template<bool forward>
343 void apply (X& v, const Y& d)
344 {
345 _preconditioner->template apply<forward>(v,d);
346 _communication.addOwnerCopyToOwnerCopy(v,v);
347 }
348
354 virtual void post (domain_type& x)
355 {
356 _preconditioner->post(x);
357 }
358
361 {
363 }
364
365 private:
367 std::shared_ptr<P> _preconditioner;
368
370 const communication_type& _communication;
371 };
372
375} // end namespace
376
377#endif
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
A simple timing class.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)