Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel
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__JACOBI__
34 : #define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__JACOBI__
35 :
36 : #include "lib_algebra/operator/interface/preconditioner.h"
37 : #include "lib_algebra/small_algebra/additional_math.h"
38 : #include "lib_algebra/cpu_algebra/vector.h"
39 :
40 : #ifdef UG_PARALLEL
41 : #include "lib_algebra/parallelization/parallelization.h"
42 : #endif
43 :
44 : namespace ug{
45 :
46 : /////////////////////////////////////////////////////////////////////////////////////////////
47 : /// Jacobi-Iteration
48 : /**
49 : * Here, the Jacobi-iteration is described for solving the linear equation
50 : *
51 : * \f$ A * x = b. A \in R^{nxn}, x \in R^n, b \in R^n \f$.
52 : *
53 : * Most of the common linear iteration-methods base on the decomposition of A into
54 : * its diagonal (D) and strict-upper(-U) and strict-lower part (-L),
55 : *
56 : * \f$ A = D - L - U \f$.
57 : *
58 : * Among others, W. Hackbusch ('Iterative Loesung grosser Gleichungssysteme'),
59 : * distinguishes three different forms for describing a linear iteration scheme.
60 : * The general 'first normal-form' of a linear iteration scheme takes the form
61 : *
62 : * \f$ x^{m+1} = M * x^m + N * b \f$,
63 : *
64 : * with some Matrices \f$ M \f$ and \f$ N \in R^{nxn} \f$. m denotes the iteration index.
65 : * The general 'second normal-form' of a linear iteration scheme takes the form
66 : *
67 : * \f$ x^{m+1} = x^m - N * (A * x^m - b) \f$.
68 : *
69 : * Those linear iteration schemes, which can be represented by the second normal-form
70 : * are the linear, consistent iteration schemes.
71 : *
72 : * Introducing the correction \f$ c{m+1} := x^{m+1} - x^m \f$ and the defect
73 : * \f$ d^m := b - A * x^m \f$ the second normal-form can be rewritten as
74 : *
75 : * \f$ c = N * d \f$.
76 : *
77 : * The matrix of the second normal-form for the Jacobi-method takes the simple form
78 : *
79 : * \f$ N = D^{-1} \f$. .
80 : *
81 : * References:
82 : * <ul>
83 : * <li> W. Hackbusch. Iterative Loesung grosser Gleichungssysteme
84 : * </ul>
85 : */
86 :
87 :
88 : /// Jacobi Preconditioner
89 : template <typename TAlgebra>
90 : class Jacobi : public IPreconditioner<TAlgebra>
91 : {
92 : public:
93 : /// Algebra type
94 : typedef TAlgebra algebra_type;
95 :
96 : /// Vector type
97 : typedef typename TAlgebra::vector_type vector_type;
98 :
99 : /// Matrix type
100 : typedef typename TAlgebra::matrix_type matrix_type;
101 :
102 : /// Matrix Operator type
103 : typedef typename IPreconditioner<TAlgebra>::matrix_operator_type matrix_operator_type;
104 :
105 : /// Base type
106 : typedef IPreconditioner<TAlgebra> base_type;
107 :
108 : protected:
109 : using base_type::set_debug;
110 : using base_type::debug_writer;
111 : using base_type::write_debug;
112 : using base_type::damping;
113 : using base_type::approx_operator;
114 :
115 : public:
116 : /// default constructor
117 0 : Jacobi() {this->set_damp(1.0); m_bBlock = true;};
118 :
119 : /// constructor setting the damping parameter
120 0 : Jacobi(number damp) {this->set_damp(damp); m_bBlock = true;};
121 :
122 : /// clone constructor
123 0 : Jacobi( const Jacobi<TAlgebra> &parent )
124 0 : : base_type(parent)
125 : {
126 0 : set_block(parent.m_bBlock);
127 : }
128 :
129 : /// Clone
130 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
131 : {
132 0 : return make_sp(new Jacobi<algebra_type>(*this));
133 : }
134 :
135 :
136 : /// returns if parallel solving is supported
137 0 : virtual bool supports_parallel() const {return true;}
138 :
139 :
140 : /// Destructor
141 0 : virtual ~Jacobi()
142 0 : {};
143 :
144 : /// sets if blocked jacobi is used (inverting block-diagonal), or plain (scalar) diagonal if false
145 : void set_block(bool b)
146 : {
147 0 : m_bBlock = b;
148 : }
149 :
150 : protected:
151 : /// Name of preconditioner
152 0 : virtual const char* name() const {return "Jacobi";}
153 :
154 : /// Preprocess routine
155 0 : virtual bool preprocess(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp)
156 : {
157 :
158 : PROFILE_BEGIN_GROUP(Jacobi_preprocess, "algebra Jacobi");
159 :
160 0 : matrix_type &mat = *pOp;
161 : // create help vector to apply diagonal
162 : size_t size = mat.num_rows();
163 0 : if(size != mat.num_cols())
164 : {
165 : UG_LOG("Square Matrix needed for Jacobi Iteration.\n");
166 0 : return false;
167 : }
168 :
169 : // resize
170 0 : m_diagInv.resize(size);
171 : #ifdef UG_PARALLEL
172 : // temporary vector for the diagonal
173 : ParallelVector<Vector< typename matrix_type::value_type > > diag;
174 : diag.resize(size);
175 :
176 : // copy the layouts+communicator into the vector
177 : diag.set_layouts(mat.layouts());
178 :
179 : // copy diagonal
180 : for(size_t i = 0; i < diag.size(); ++i){
181 : diag[i] = mat(i, i);
182 : }
183 :
184 : // make diagonal consistent
185 : diag.set_storage_type(PST_ADDITIVE);
186 : diag.change_storage_type(PST_CONSISTENT);
187 :
188 :
189 : if(diag.size() > 0)
190 : if(CheckVectorInvertible(diag) == false)
191 : return false;
192 : // UG_ASSERT(CheckVectorInvertible(diag), "Jacobi: A has noninvertible diagonal");
193 :
194 : #endif
195 :
196 : // get damping in constant case to damp at once
197 : number damp = 1.0;
198 0 : if(damping()->constant_damping())
199 0 : damp = damping()->damping();
200 :
201 : typename matrix_type::value_type m;
202 : // invert diagonal and multiply by damping
203 0 : for(size_t i = 0; i < mat.num_rows(); ++i)
204 : {
205 : #ifdef UG_PARALLEL
206 : typename matrix_type::value_type &d = diag[i];
207 : #else
208 : typename matrix_type::value_type &d = mat(i,i);
209 : #endif
210 0 : if(!m_bBlock)
211 0 : GetDiag(m, d);
212 : else
213 0 : m = d;
214 0 : m *= 1./damp;
215 : GetInverse(m_diagInv[i], m);
216 : }
217 :
218 : // done
219 : return true;
220 : }
221 :
222 0 : virtual bool step(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp, vector_type& c, const vector_type& d)
223 : {
224 : PROFILE_BEGIN_GROUP(Jacobi_step, "algebra Jacobi");
225 :
226 : // multiply defect with diagonal, c = damp * D^{-1} * d
227 : // note, that the damping is already included in the inverse diagonal
228 0 : for(size_t i = 0; i < m_diagInv.size(); ++i)
229 : {
230 : // c[i] = m_diagInv[i] * d[i];
231 0 : MatMult(c[i], 1.0, m_diagInv[i], d[i]);
232 : }
233 :
234 : #ifdef UG_PARALLEL
235 :
236 : // the computed correction is additive
237 : c.set_storage_type(PST_ADDITIVE);
238 :
239 : // we make it consistent
240 : if(!c.change_storage_type(PST_CONSISTENT))
241 : {
242 : UG_LOG("ERROR in 'JacobiPreconditioner::apply': "
243 : "Cannot change parallel status of correction to consistent.\n");
244 : return false;
245 : }
246 : #endif
247 : // done
248 0 : return true;
249 : }
250 :
251 : /// Postprocess routine
252 0 : virtual bool postprocess() {return true;}
253 :
254 :
255 : // overwrite function in order to specially treat constant damping
256 0 : virtual bool apply(vector_type& c, const vector_type& d)
257 : {
258 : PROFILE_BEGIN_GROUP(Jacobi_apply, "algebra Jacobi");
259 : // Check that operator is initialized
260 0 : if(!this->m_bInit)
261 : {
262 0 : UG_LOG("ERROR in '"<<name()<<"::apply': Iterator not initialized.\n");
263 0 : return false;
264 : }
265 :
266 : // Check parallel status
267 : #ifdef UG_PARALLEL
268 : if(!d.has_storage_type(PST_ADDITIVE))
269 : UG_THROW(name() << "::apply: Wrong parallel "
270 : "storage format. Defect must be additive.");
271 : #endif
272 :
273 : // Check sizes
274 0 : THROW_IF_NOT_EQUAL_4(c.size(), d.size(), approx_operator()->num_rows(), approx_operator()->num_cols());
275 :
276 : // apply iterator: c = B*d
277 0 : if(!step(approx_operator(), c, d))
278 : {
279 0 : UG_LOG("ERROR in '"<<name()<<"::apply': Step Routine failed.\n");
280 0 : return false;
281 : }
282 :
283 : // apply scaling
284 0 : if(!damping()->constant_damping()){
285 0 : const number kappa = damping()->damping(c, d, approx_operator());
286 0 : if(kappa != 1.0){
287 : c *= kappa;
288 : }
289 : }
290 :
291 : // Correction is always consistent
292 : #ifdef UG_PARALLEL
293 : if(!c.change_storage_type(PST_CONSISTENT))
294 : UG_THROW(name() << "::apply': Cannot change "
295 : "parallel storage type of correction to consistent.");
296 : #endif
297 :
298 : // we're done
299 : return true;
300 : }
301 :
302 : protected:
303 : /// type of block-inverse
304 : typedef typename block_traits<typename matrix_type::value_type>::inverse_type inverse_type;
305 :
306 : /// storage of the inverse diagonal in parallel
307 : std::vector<inverse_type> m_diagInv;
308 : bool m_bBlock;
309 :
310 :
311 : };
312 :
313 : } // end namespace ug
314 :
315 : //#include "gpujacobi.h"
316 :
317 : #endif
|
