ug4
Loading...
Searching...
No Matches
jacobi.h
Go to the documentation of this file.
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
39
40#ifdef UG_PARALLEL
42#endif
43
44namespace ug{
45
48
89template <typename TAlgebra>
90class Jacobi : public IPreconditioner<TAlgebra>
91{
92 public:
94 typedef TAlgebra algebra_type;
95
97 typedef typename TAlgebra::vector_type vector_type;
98
100 typedef typename TAlgebra::matrix_type matrix_type;
101
104
107
108 protected:
112 using base_type::damping;
114
115 public:
117 Jacobi() {this->set_damp(1.0); m_bBlock = true;};
118
120 Jacobi(number damp) {this->set_damp(damp); m_bBlock = true;};
121
123 Jacobi( const Jacobi<TAlgebra> &parent )
124 : base_type(parent)
125 {
126 set_block(parent.m_bBlock);
127 }
128
131 {
132 return make_sp(new Jacobi<algebra_type>(*this));
133 }
134
135
137 virtual bool supports_parallel() const {return true;}
138
139
141 virtual ~Jacobi()
142 {};
143
145 void set_block(bool b)
146 {
147 m_bBlock = b;
148 }
149
150 protected:
152 virtual const char* name() const {return "Jacobi";}
153
156 {
157
158 PROFILE_BEGIN_GROUP(Jacobi_preprocess, "algebra Jacobi");
159
160 matrix_type &mat = *pOp;
161 // create help vector to apply diagonal
162 size_t size = mat.num_rows();
163 if(size != mat.num_cols())
164 {
165 UG_LOG("Square Matrix needed for Jacobi Iteration.\n");
166 return false;
167 }
168
169 // resize
170 m_diagInv.resize(size);
171#ifdef UG_PARALLEL
172 // temporary vector for the diagonal
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
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 if(damping()->constant_damping())
199 damp = damping()->damping();
200
201 typename matrix_type::value_type m;
202 // invert diagonal and multiply by damping
203 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 if(!m_bBlock)
211 GetDiag(m, d);
212 else
213 m = d;
214 m *= 1./damp;
215 GetInverse(m_diagInv[i], m);
216 }
217
218 // done
219 return true;
220 }
221
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 for(size_t i = 0; i < m_diagInv.size(); ++i)
229 {
230 // c[i] = m_diagInv[i] * d[i];
231 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 return true;
249 }
250
252 virtual bool postprocess() {return true;}
253
254
255 // overwrite function in order to specially treat constant damping
256 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 if(!this->m_bInit)
261 {
262 UG_LOG("ERROR in '"<<name()<<"::apply': Iterator not initialized.\n");
263 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
275
276 // apply iterator: c = B*d
277 if(!step(approx_operator(), c, d))
278 {
279 UG_LOG("ERROR in '"<<name()<<"::apply': Step Routine failed.\n");
280 return false;
281 }
282
283 // apply scaling
284 if(!damping()->constant_damping()){
285 const number kappa = damping()->damping(c, d, approx_operator());
286 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:
305
307 std::vector<inverse_type> m_diagInv;
309
310
311};
312
313} // end namespace ug
314
315//#include "gpujacobi.h"
316
317#endif
Definition smart_pointer.h:108
virtual void set_debug(SmartPtr< IDebugWriter< algebra_type > > spDebugWriter)
set debug writer
Definition debug_writer.h:384
void write_debug(const matrix_type &mat, const char *filename)
write debug output for a matrix (if debug writer set)
Definition debug_writer.h:399
SmartPtr< IDebugWriter< algebra_type > > debug_writer()
returns the debug writer
Definition debug_writer.h:391
void set_damp(SmartPtr< IDamping< TAlgebra::vector_type, TAlgebra::vector_type > > spScaling)
sets a scaling for the correction
Definition linear_iterator.h:163
SmartPtr< IDamping< TAlgebra::vector_type, TAlgebra::vector_type > > damping()
returns the scaling
Definition linear_iterator.h:173
describes a linear iterator that is based on a matrix operator
Definition preconditioner.h:103
SmartPtr< MatrixOperator< matrix_type, vector_type > > approx_operator()
underlying matrix based operator used for the preconditioner
Definition preconditioner.h:372
bool m_bInit
init flag indicating if init has been called
Definition preconditioner.h:385
Jacobi-Iteration.
Definition jacobi.h:91
SmartPtr< MatrixOperator< matrix_type, vector_type > > approx_operator()
underlying matrix based operator used for the preconditioner
Definition preconditioner.h:372
bool m_bBlock
Definition jacobi.h:308
virtual bool postprocess()
Postprocess routine.
Definition jacobi.h:252
TAlgebra::matrix_type matrix_type
Matrix type.
Definition jacobi.h:100
virtual bool preprocess(SmartPtr< MatrixOperator< matrix_type, vector_type > > pOp)
Preprocess routine.
Definition jacobi.h:155
virtual bool step(SmartPtr< MatrixOperator< matrix_type, vector_type > > pOp, vector_type &c, const vector_type &d)
computes a new correction c = B*d
Definition jacobi.h:222
virtual const char * name() const
Name of preconditioner.
Definition jacobi.h:152
TAlgebra::vector_type vector_type
Vector type.
Definition jacobi.h:97
virtual bool apply(vector_type &c, const vector_type &d)
compute new correction c = B*d
Definition jacobi.h:256
virtual SmartPtr< ILinearIterator< vector_type > > clone()
Clone.
Definition jacobi.h:130
IPreconditioner< TAlgebra > base_type
Base type.
Definition jacobi.h:106
Jacobi(number damp)
constructor setting the damping parameter
Definition jacobi.h:120
std::vector< inverse_type > m_diagInv
storage of the inverse diagonal in parallel
Definition jacobi.h:307
virtual bool supports_parallel() const
returns if parallel solving is supported
Definition jacobi.h:137
IPreconditioner< TAlgebra >::matrix_operator_type matrix_operator_type
Matrix Operator type.
Definition jacobi.h:103
Jacobi(const Jacobi< TAlgebra > &parent)
clone constructor
Definition jacobi.h:123
TAlgebra algebra_type
Algebra type.
Definition jacobi.h:94
block_traits< typenamematrix_type::value_type >::inverse_type inverse_type
type of block-inverse
Definition jacobi.h:304
virtual ~Jacobi()
Destructor.
Definition jacobi.h:141
Jacobi()
default constructor
Definition jacobi.h:117
void set_block(bool b)
sets if blocked jacobi is used (inverting block-diagonal), or plain (scalar) diagonal if false
Definition jacobi.h:145
Definition parallel_vector.h:60
bool change_storage_type(ParallelStorageType type)
changes to the requested storage type if possible
Definition parallel_vector_impl.h:118
void set_storage_type(uint type)
sets the storage type
Definition parallel_vector.h:104
void set_layouts(ConstSmartPtr< AlgebraLayouts > layouts)
sets the algebra layouts
Definition parallel_vector.h:100
@ PST_CONSISTENT
Definition parallel_storage_type.h:68
@ PST_ADDITIVE
Definition parallel_storage_type.h:69
bool CheckVectorInvertible(const TVector &v)
Definition sparsematrix_util.h:1116
#define THROW_IF_NOT_EQUAL_4(s1, s2, s3, s4)
Definition error.h:183
#define UG_THROW(msg)
Definition error.h:57
#define UG_LOG(msg)
Definition log.h:367
double number
Definition types.h:124
the ug namespace
bool GetInverse(block_traits< T >::inverse_type &inv, const T &m)
void GetDiag(double &a, double b)
Definition additional_math.h:144
bool MatMult(vector_t &dest, const number &beta1, const matrix_t &A1, const vector_t &w1)
calculates dest = beta1 * A1;
Definition operations_mat.h:59
#define PROFILE_BEGIN_GROUP(name, groups)
Definition profiler.h:255
SmartPtr< T, FreePolicy > make_sp(T *inst)
returns a SmartPtr for the passed raw pointer
Definition smart_pointer.h:836
size_t num_rows() const
Definition sparsematrix_interface.h:38
size_t num_cols() const
Definition sparsematrix_interface.h:39
Definition communication_policies.h:58