Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Martin Rupp
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__GAUSS_SEIDEL__
34 : #define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__GAUSS_SEIDEL__
35 :
36 : #include "lib_algebra/operator/interface/preconditioner.h"
37 : #include "lib_algebra/algebra_common/core_smoothers.h"
38 : #include "lib_algebra/algebra_common/sparsematrix_util.h"
39 : #ifdef UG_PARALLEL
40 : #include "lib_algebra/parallelization/parallelization.h"
41 : #include "lib_algebra/parallelization/matrix_overlap.h"
42 : #include "lib_algebra/parallelization/parallel_matrix_overlap_impl.h"
43 : #endif
44 :
45 : #include "lib_algebra/ordering_strategies/algorithms/IOrderingAlgorithm.h"
46 : #include "lib_algebra/algebra_common/permutation_util.h"
47 :
48 : namespace ug{
49 :
50 : template<typename TAlgebra>
51 : class GaussSeidelBase : public IPreconditioner<TAlgebra>
52 : {
53 : public:
54 : // Algebra type
55 : typedef TAlgebra algebra_type;
56 :
57 : // Vector type
58 : typedef typename TAlgebra::vector_type vector_type;
59 :
60 : // Matrix type
61 : typedef typename TAlgebra::matrix_type matrix_type;
62 :
63 : /// Matrix Operator type
64 : typedef typename IPreconditioner<TAlgebra>::matrix_operator_type matrix_operator_type;
65 :
66 : /// Base type
67 : typedef IPreconditioner<TAlgebra> base_type;
68 :
69 : /// Ordering type
70 : typedef std::vector<size_t> ordering_container_type;
71 : typedef IOrderingAlgorithm<TAlgebra, ordering_container_type> ordering_algo_type;
72 :
73 : protected:
74 : using base_type::set_debug;
75 : using base_type::debug_writer;
76 : using base_type::write_debug;
77 :
78 : public:
79 : // Constructor
80 0 : GaussSeidelBase() :
81 0 : m_relax(1.0),
82 0 : m_bConsistentInterfaces(false),
83 0 : m_useOverlap(false) {};
84 :
85 : /// clone constructor
86 0 : GaussSeidelBase( const GaussSeidelBase<TAlgebra> &parent )
87 : : base_type(parent),
88 0 : m_bConsistentInterfaces(parent.m_bConsistentInterfaces),
89 0 : m_useOverlap(parent.m_useOverlap),
90 0 : m_spOrderingAlgo(parent.m_spOrderingAlgo)
91 : {
92 0 : set_sor_relax(parent.m_relax);
93 0 : }
94 :
95 : /// set relaxation parameter to define a SOR-method
96 0 : void set_sor_relax(number relaxFactor){ m_relax = relaxFactor;}
97 :
98 : /// activates the new parallelization approach (disabled by default)
99 0 : void enable_consistent_interfaces(bool enable) {m_bConsistentInterfaces = enable;}
100 :
101 0 : void enable_overlap (bool enable) {m_useOverlap = enable;}
102 :
103 : /// sets an ordering algorithm
104 : void set_ordering_algorithm(SmartPtr<ordering_algo_type> ordering_algo){
105 : m_spOrderingAlgo = ordering_algo;
106 : }
107 :
108 : virtual const char* name() const = 0;
109 : protected:
110 :
111 0 : virtual bool supports_parallel() const {return true;}
112 :
113 : // Preprocess routine
114 0 : virtual bool preprocess(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp)
115 : {
116 : PROFILE_BEGIN_GROUP(GaussSeidel_preprocess, "algebra gaussseidel");
117 : matrix_type *pA;
118 : #ifdef UG_PARALLEL
119 : if(pcl::NumProcs() > 1)
120 : {
121 : if(m_useOverlap){
122 : m_A = *pOp;
123 : CreateOverlap(m_A);
124 : m_oD.set_layouts(m_A.layouts());
125 : m_oC.set_layouts(m_A.layouts());
126 : m_oD.resize(m_A.num_rows(), false);
127 : m_oC.resize(m_A.num_rows(), false);
128 : }
129 : else if (m_bConsistentInterfaces)
130 : {
131 : m_A = *pOp;
132 : MatMakeConsistentOverlap0(m_A);
133 : }
134 : else
135 : {
136 : // copy original matrix
137 : MakeConsistent(*pOp, m_A);
138 : // set zero on slaves
139 : std::vector<IndexLayout::Element> vIndex;
140 : CollectUniqueElements(vIndex, m_A.layouts()->slave());
141 : SetDirichletRow(m_A, vIndex);
142 : }
143 : pA = &m_A;
144 : }
145 : else
146 : pA = &(*pOp);
147 : #else
148 : pA = &(*pOp);
149 : #endif
150 0 : THROW_IF_NOT_EQUAL(pA->num_rows(), pA->num_cols());
151 : // UG_ASSERT(CheckDiagonalInvertible(A), "GS: A has noninvertible diagonal");
152 : UG_COND_THROW(CheckDiagonalInvertible(*pA) == false, name() << ": A has noninvertible diagonal");
153 :
154 0 : return true;
155 : }
156 :
157 : // Postprocess routine
158 0 : virtual bool postprocess() {return true;}
159 :
160 : virtual void step(const matrix_type &A, vector_type &c, const vector_type &d, const number relax) = 0;
161 :
162 : // Stepping routine
163 0 : virtual bool step(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp, vector_type& c, const vector_type& d)
164 : {
165 : PROFILE_BEGIN_GROUP(GaussSeidel_step, "algebra gaussseidel");
166 :
167 : #ifdef UG_PARALLEL
168 : if(pcl::NumProcs() > 1)
169 : {
170 : if(m_useOverlap){
171 : for(size_t i = 0; i < d.size(); ++i)
172 : m_oD[i] = d[i];
173 : for(size_t i = d.size(); i < m_oD.size(); ++i)
174 : m_oD[i] = 0;
175 : m_oD.set_storage_type(PST_ADDITIVE);
176 : m_oD.change_storage_type(PST_CONSISTENT);
177 :
178 : step(m_A, m_oC, m_oD, m_relax);
179 :
180 : for(size_t i = 0; i < c.size(); ++i)
181 : c[i] = m_oC[i];
182 : SetLayoutValues(&c, c.layouts()->slave(), 0);
183 : c.set_storage_type(PST_UNIQUE);
184 : }
185 : else if (m_bConsistentInterfaces)
186 : {
187 : // make defect consistent
188 : SmartPtr<vector_type> spDtmp = d.clone();
189 : spDtmp->change_storage_type(PST_CONSISTENT);
190 :
191 : THROW_IF_NOT_EQUAL_3(c.size(), spDtmp->size(), m_A.num_rows());
192 : step(m_A, c, *spDtmp, m_relax);
193 :
194 : // declare c unique to enforce that only master correction is used
195 : // when it is made consistent below
196 : c.set_storage_type(PST_UNIQUE);
197 :
198 : //UG_COND_THROW(!d.has_storage_type(PST_ADDITIVE), "Additive or unique defect expected.");
199 : //step(m_A, c, d, m_relax);
200 : //c.set_storage_type(PST_ADDITIVE);
201 : }
202 : else
203 : {
204 : // todo: do not clone every time
205 : // make defect unique
206 : SmartPtr<vector_type> spDtmp = d.clone();
207 : spDtmp->change_storage_type(PST_UNIQUE);
208 :
209 : THROW_IF_NOT_EQUAL_3(c.size(), spDtmp->size(), m_A.num_rows());
210 : step(m_A, c, *spDtmp, m_relax);
211 : c.set_storage_type(PST_UNIQUE);
212 : }
213 :
214 : // make correction consistent
215 : c.change_storage_type(PST_CONSISTENT);
216 :
217 : return true;
218 : }
219 : else
220 : #endif
221 : {
222 0 : matrix_type &A = *pOp;
223 0 : THROW_IF_NOT_EQUAL_4(c.size(), d.size(), A.num_rows(), A.num_cols());
224 :
225 0 : step(A, c, d, m_relax);
226 : #ifdef UG_PARALLEL
227 : c.set_storage_type(PST_CONSISTENT);
228 : #endif
229 0 : return true;
230 : }
231 : }
232 :
233 : protected:
234 : #ifdef UG_PARALLEL
235 : matrix_type m_A;
236 : #endif
237 : protected:
238 : /// relaxation parameter
239 : number m_relax;
240 :
241 : /// for overlaps only
242 : vector_type m_oD;
243 : vector_type m_oC;
244 :
245 : bool m_bConsistentInterfaces;
246 : bool m_useOverlap;
247 :
248 :
249 : /// for ordering algorithms
250 : SmartPtr<ordering_algo_type> m_spOrderingAlgo;
251 : #ifdef NOT_YET
252 : ordering_container_type m_ordering, m_old_ordering;
253 : std::vector<size_t> m_newIndex, m_oldIndex;
254 : bool m_bSortIsIdentity;
255 : #endif
256 : };
257 :
258 : /// Gauss-Seidel preconditioner for the 'forward' ordering of the dofs
259 : /**
260 : * This class implements the Gauss-Seidel preconditioner (and smoother) for the
261 : * 'forward' ordering of the dofs. When a relaxation parameter is set by the method
262 : * 'set_sor_relax', the resulting preconditioner is better known as (forward) 'SOR'-method.
263 : * References:
264 : * <ul>
265 : * <li> W. Hackbusch. Iterative solution of large sparse systems of equations. New York: Springer, 1994
266 : * </ul>
267 : *
268 : * \tparam TAlgebra Algebra type
269 : */
270 : template <typename TAlgebra>
271 : class GaussSeidel : public GaussSeidelBase<TAlgebra>
272 : {
273 : typedef TAlgebra algebra_type;
274 : typedef typename TAlgebra::vector_type vector_type;
275 : typedef typename TAlgebra::matrix_type matrix_type;
276 : typedef GaussSeidelBase<TAlgebra> base_type;
277 :
278 : public:
279 : // Name of preconditioner
280 0 : virtual const char* name() const {return "Gauss-Seidel";}
281 :
282 : /// constructor
283 0 : GaussSeidel() : base_type() {}
284 :
285 : /// clone constructor
286 0 : GaussSeidel( const GaussSeidel<TAlgebra> &parent )
287 0 : : base_type(parent)
288 : { }
289 :
290 : /// Clone
291 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
292 : {
293 0 : return make_sp(new GaussSeidel<algebra_type>(*this));
294 : }
295 :
296 : // Stepping routine
297 0 : virtual void step(const matrix_type &A, vector_type &c, const vector_type &d, const number relax)
298 : {
299 0 : gs_step_LL(A, c, d, relax);
300 0 : }
301 : };
302 :
303 : /// Gauss-Seidel preconditioner for the 'backward' ordering of the dofs
304 : /**
305 : * This class implements the Gauss-Seidel preconditioner (and smoother) for the
306 : * 'backward' ordering of the dofs. When a relaxation parameter is set by the method
307 : * 'set_sor_relax', the resulting preconditioner is better known as backward 'SOR'-method.
308 : * References:
309 : * <ul>
310 : * <li> W. Hackbusch. Iterative solution of large sparse systems of equations. New York: Springer, 1994
311 : * </ul>
312 : *
313 : * \tparam TAlgebra Algebra type
314 : */
315 : template <typename TAlgebra>
316 : class BackwardGaussSeidel : public GaussSeidelBase<TAlgebra>
317 : {
318 : typedef TAlgebra algebra_type;
319 : typedef typename TAlgebra::vector_type vector_type;
320 : typedef typename TAlgebra::matrix_type matrix_type;
321 : typedef GaussSeidelBase<TAlgebra> base_type;
322 :
323 : public:
324 : // Name of preconditioner
325 0 : virtual const char* name() const {return "Backward Gauss-Seidel";}
326 :
327 : /// constructor
328 0 : BackwardGaussSeidel() : base_type() {}
329 :
330 : /// clone constructor
331 0 : BackwardGaussSeidel( const BackwardGaussSeidel<TAlgebra> &parent )
332 0 : : base_type(parent)
333 : { }
334 :
335 : /// Clone
336 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
337 : {
338 0 : return make_sp(new BackwardGaussSeidel<algebra_type>(*this));
339 : }
340 :
341 : // Stepping routine
342 0 : virtual void step(const matrix_type &A, vector_type &c, const vector_type &d, const number relax)
343 : {
344 0 : gs_step_UR(A, c, d, relax);
345 0 : }
346 : };
347 :
348 :
349 : template <typename TAlgebra>
350 : class SymmetricGaussSeidel : public GaussSeidelBase<TAlgebra>
351 : {
352 : typedef TAlgebra algebra_type;
353 : typedef typename TAlgebra::vector_type vector_type;
354 : typedef typename TAlgebra::matrix_type matrix_type;
355 : typedef GaussSeidelBase<TAlgebra> base_type;
356 :
357 : public:
358 : // Name of preconditioner
359 0 : virtual const char* name() const {return "Symmetric Gauss-Seidel";}
360 :
361 : /// constructor
362 0 : SymmetricGaussSeidel() : base_type() {}
363 :
364 : /// clone constructor
365 0 : SymmetricGaussSeidel( const SymmetricGaussSeidel<TAlgebra> &parent )
366 0 : : base_type(parent)
367 : { }
368 :
369 : /// Clone
370 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
371 : {
372 0 : return make_sp(new SymmetricGaussSeidel<algebra_type>(*this));
373 : }
374 :
375 : // Stepping routine
376 0 : virtual void step(const matrix_type &A, vector_type &c, const vector_type &d, const number relax)
377 : {
378 0 : sgs_step(A, c, d, relax);
379 0 : }
380 : };
381 :
382 : } // end namespace ug
383 :
384 : #endif // __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__GAUSS_SEIDEL__
|
