Dune Core Modules (2.6.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 "io.hh"
16#include "bvector.hh"
17#include "vbvector.hh"
18#include "bcrsmatrix.hh"
19#include "io.hh"
20#include "gsetc.hh"
21#include "ilu.hh"
22#include "operators.hh"
23#include "solvers.hh"
24#include "preconditioners.hh"
25#include "scalarproducts.hh"
26#include "owneroverlapcopy.hh"
27
28namespace Dune {
29
58 template<class M, class X, class Y, class C>
60 {
61 public:
63 typedef M matrix_type;
65 typedef X domain_type;
67 typedef Y range_type;
69 typedef typename X::field_type field_type;
72
73 typedef typename C::PIS PIS;
74 typedef typename C::RI RI;
75 typedef typename RI::RemoteIndexList RIL;
76 typedef typename RI::const_iterator RIIterator;
77 typedef typename RIL::const_iterator RILIterator;
78 typedef typename M::ConstColIterator ColIterator;
79 typedef typename M::ConstRowIterator RowIterator;
80 typedef std::multimap<int,int> MM;
81 typedef std::multimap<int,std::pair<int,RILIterator> > RIMap;
82 typedef typename RIMap::iterator RIMapit;
83
92 : _A_(A), communication(com), buildcomm(true)
93 {}
94
96 virtual void apply (const X& x, Y& y) const
97 {
98 y = 0;
99 novlp_op_apply(x,y,1);
100 communication.addOwnerCopyToOwnerCopy(y,y);
101 }
102
104 virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
105 {
106 // only apply communication to alpha*A*x to make it consistent,
107 // y already has to be consistent.
108 Y y1(y);
109 y = 0;
110 novlp_op_apply(x,y,alpha);
111 communication.addOwnerCopyToOwnerCopy(y,y);
112 y += y1;
113 }
114
116 virtual const matrix_type& getmat () const
117 {
118 return _A_;
119 }
120
121 void novlp_op_apply (const X& x, Y& y, field_type alpha) const
122 {
123 //get index sets
124 const PIS& pis=communication.indexSet();
125 const RI& ri = communication.remoteIndices();
126
127 // at the beginning make a multimap "bordercontribution".
128 // process has i and j as border dofs but is not the owner
129 // => only contribute to Ax if i,j is in bordercontribution
130 if (buildcomm == true) {
131
132 // set up mask vector
133 if (mask.size()!=static_cast<typename std::vector<double>::size_type>(x.size())) {
134 mask.resize(x.size());
135 for (typename std::vector<double>::size_type i=0; i<mask.size(); i++)
136 mask[i] = 1;
137 for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
138 if (i->local().attribute()==OwnerOverlapCopyAttributeSet::copy)
139 mask[i->local().local()] = 0;
140 else if (i->local().attribute()==OwnerOverlapCopyAttributeSet::overlap)
141 mask[i->local().local()] = 2;
142 }
143
144 for (MM::iterator iter = bordercontribution.begin();
145 iter != bordercontribution.end(); ++iter)
146 bordercontribution.erase(iter);
147 std::map<int,int> owner; //key: local index i, value: process, that owns i
148 RIMap rimap;
149
150 // for each local index make multimap rimap:
151 // key: local index i, data: pair of process that knows i and pointer to RI entry
152 for (RowIterator i = _A_.begin(); i != _A_.end(); ++i)
153 if (mask[i.index()] == 0)
154 for (RIIterator remote = ri.begin(); remote != ri.end(); ++remote) {
155 RIL& ril = *(remote->second.first);
156 for (RILIterator rindex = ril.begin(); rindex != ril.end(); ++rindex)
157 if (rindex->attribute() != OwnerOverlapCopyAttributeSet::overlap)
158 if (rindex->localIndexPair().local().local() == i.index()) {
159 rimap.insert
160 (std::make_pair(i.index(),
161 std::pair<int,RILIterator>(remote->first, rindex)));
162 if(rindex->attribute()==OwnerOverlapCopyAttributeSet::owner)
163 owner.insert(std::make_pair(i.index(),remote->first));
164 }
165 }
166
167 int iowner = 0;
168 for (RowIterator i = _A_.begin(); i != _A_.end(); ++i) {
169 if (mask[i.index()] == 0) {
170 std::map<int,int>::iterator it = owner.find(i.index());
171 iowner = it->second;
172 std::pair<RIMapit, RIMapit> foundiit = rimap.equal_range(i.index());
173 for (ColIterator j = _A_[i.index()].begin(); j != _A_[i.index()].end(); ++j) {
174 if (mask[j.index()] == 0) {
175 bool flag = true;
176 for (RIMapit foundi = foundiit.first; foundi != foundiit.second; ++foundi) {
177 std::pair<RIMapit, RIMapit> foundjit = rimap.equal_range(j.index());
178 for (RIMapit foundj = foundjit.first; foundj != foundjit.second; ++foundj)
179 if (foundj->second.first == foundi->second.first)
180 if (foundj->second.second->attribute() == OwnerOverlapCopyAttributeSet::owner
181 || foundj->second.first == iowner
182 || foundj->second.first < communication.communicator().rank()) {
183 flag = false;
184 continue;
185 }
186 if (flag == false)
187 continue;
188 }
189 // donĀ“t contribute to Ax if
190 // 1. the owner of j has i as interior/border dof
191 // 2. iowner has j as interior/border dof
192 // 3. there is another process with smaller rank that has i and j
193 // as interor/border dofs
194 // if the owner of j does not have i as interior/border dof,
195 // it will not be taken into account
196 if (flag==true)
197 bordercontribution.insert(std::pair<int,int>(i.index(),j.index()));
198 }
199 }
200 }
201 }
202 buildcomm = false;
203 }
204
205 //compute alpha*A*x nonoverlapping case
206 for (RowIterator i = _A_.begin(); i != _A_.end(); ++i) {
207 if (mask[i.index()] == 0) {
208 //dof doesn't belong to process but is border (not ghost)
209 for (ColIterator j = _A_[i.index()].begin(); j != _A_[i.index()].end(); ++j) {
210 if (mask[j.index()] == 1) //j is owner => then sum entries
211 (*j).usmv(alpha,x[j.index()],y[i.index()]);
212 else if (mask[j.index()] == 0) {
213 std::pair<MM::iterator, MM::iterator> itp =
214 bordercontribution.equal_range(i.index());
215 for (MM::iterator it = itp.first; it != itp.second; ++it)
216 if ((*it).second == (int)j.index())
217 (*j).usmv(alpha,x[j.index()],y[i.index()]);
218 }
219 }
220 }
221 else if (mask[i.index()] == 1) {
222 for (ColIterator j = _A_[i.index()].begin(); j != _A_[i.index()].end(); ++j)
223 if (mask[j.index()] != 2)
224 (*j).usmv(alpha,x[j.index()],y[i.index()]);
225 }
226 }
227 }
228
231 {
233 }
234
235 private:
236 const matrix_type& _A_;
237 const communication_type& communication;
238 mutable bool buildcomm;
239 mutable std::vector<double> mask;
240 mutable std::multimap<int,int> bordercontribution;
241 };
242
245 namespace Amg
246 {
247 template<class T> class ConstructionTraits;
248 }
249
264 template<class C, class P>
266 : public Dune::Preconditioner<typename P::domain_type,typename P::range_type> {
268 public:
270 typedef typename P::domain_type domain_type;
272 typedef typename P::range_type range_type;
275
284 : preconditioner(prec), communication(c)
285 {}
286
292 virtual void pre (domain_type& x, range_type& b)
293 {
294 preconditioner.pre(x,b);
295 }
296
302 virtual void apply (domain_type& v, const range_type& d)
303 {
304 // block preconditioner equivalent to WrappedPreconditioner from
305 // pdelab/backend/ovlpistsolverbackend.hh,
306 // but not to BlockPreconditioner from schwarz.hh
307 preconditioner.apply(v,d);
308 communication.addOwnerCopyToOwnerCopy(v,v);
309 }
310
316 virtual void post (domain_type& x)
317 {
318 preconditioner.post(x);
319 }
320
323 {
325 }
326
327 private:
329 P& preconditioner;
330
332 const communication_type& communication;
333 };
334
337} // end namespace
338
339#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...
Traits class for generically constructing non default constructable types.
Definition: construction.hh:38
A linear operator exporting itself in matrix form.
Definition: operators.hh:106
Nonoverlapping parallel preconditioner.
Definition: novlpschwarz.hh:266
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: novlpschwarz.hh:322
virtual void apply(domain_type &v, const range_type &d)
Apply the preconditioner.
Definition: novlpschwarz.hh:302
P::range_type range_type
The range type of the preconditioner.
Definition: novlpschwarz.hh:272
virtual void post(domain_type &x)
Clean up.
Definition: novlpschwarz.hh:316
C communication_type
The type of the communication object.
Definition: novlpschwarz.hh:274
NonoverlappingBlockPreconditioner(P &prec, const communication_type &c)
Constructor.
Definition: novlpschwarz.hh:283
virtual void pre(domain_type &x, range_type &b)
Prepare the preconditioner.
Definition: novlpschwarz.hh:292
P::domain_type domain_type
The domain type of the preconditioner.
Definition: novlpschwarz.hh:270
A nonoverlapping operator with communication object.
Definition: novlpschwarz.hh:60
C communication_type
The type of the communication object.
Definition: novlpschwarz.hh:71
virtual void apply(const X &x, Y &y) const
apply operator to x:
Definition: novlpschwarz.hh:96
X domain_type
The type of the domain.
Definition: novlpschwarz.hh:65
virtual SolverCategory::Category category() const
Category of the linear operator (see SolverCategory::Category)
Definition: novlpschwarz.hh:230
virtual const matrix_type & getmat() const
get matrix via *
Definition: novlpschwarz.hh:116
Y range_type
The type of the range.
Definition: novlpschwarz.hh:67
M matrix_type
The type of the matrix we operate on.
Definition: novlpschwarz.hh:63
virtual void applyscaleadd(field_type alpha, const X &x, Y &y) const
apply operator to x, scale and add:
Definition: novlpschwarz.hh:104
X::field_type field_type
The field type of the range.
Definition: novlpschwarz.hh:69
NonoverlappingSchwarzOperator(const matrix_type &A, const communication_type &com)
constructor: just store a reference to a matrix.
Definition: novlpschwarz.hh:91
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:30
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
Some generic functions for pretty printing vectors and matrices.
Dune namespace.
Definition: alignedallocator.hh:10
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.
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)