Line data Source code
1 : /*
2 : * Copyright (c) 2013-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__LIB_ALGEBRA__OPERATOR__INTERFACE__PRECONDITIONER__
34 : #define __H__LIB_ALGEBRA__OPERATOR__INTERFACE__PRECONDITIONER__
35 :
36 : #include "linear_iterator.h"
37 : #include "matrix_operator.h"
38 : #include "lib_algebra/operator/damping.h"
39 : #include "lib_algebra/operator/debug_writer.h"
40 :
41 : #ifdef UG_PARALLEL
42 : #include "lib_algebra/parallelization/parallelization.h"
43 : #endif
44 :
45 : namespace ug{
46 : ///////////////////////////////////////////////////////////////////////////////
47 : // Matrix Iterator Operator (Preconditioner)
48 : ///////////////////////////////////////////////////////////////////////////////
49 :
50 : /// describes a linear iterator that is based on a matrix operator
51 : /**
52 : * This class is the base class for all linear iterators, that act on matrix
53 : * based linear operators. We call the matrix based iterator a 'Preconditioner'.
54 : * They are used in iterative schemes when solving a linear system.
55 : * Usually, a linear problem like L*u = f is intended to be solved. This is
56 : * done in an iterative way by performing an iteration of
57 : *
58 : * start: compute d := f - L*u
59 : * iterate: - c := B*d (compute correction)
60 : * - u := u + c (update solution)
61 : * - d := d - L*c (update defect)
62 : *
63 : * This iterator class describes the application of B in the scheme above, where
64 : * L is internally represented by a matrix.
65 : *
66 : * This method derives from the ILinearIterator-interface and thus must
67 : * implement these to steps:
68 : *
69 : * 1. init(L, u) or init(L):
70 : * These methods initialize the iterator and one of these methods has to
71 : * be called before any of the apply methods can be used. Passing the
72 : * linear operator indicates that this operator is used as underlying
73 : * for the iterator.
74 : *
75 : * 2. apply or apply_return_defect:
76 : * These methods are used to compute the correction (and to update the
77 : * defect at the same time). Note, that these methods can only be called
78 : * when the iterator has been initialized.
79 : *
80 : * This splitting has been made, since initialization may be computationally
81 : * expansive. Thus, the user of this class has the choice when to call this
82 : * initialization. E.g. when the operator is applied several times the init of
83 : * the iterator is only needed once.
84 : *
85 : * In order to facilitate the implementation of derived classes a default
86 : * implementation of these virtual method is given. Thus, the implementational
87 : * part for a matrix based iterator is to implement the three functions:
88 : *
89 : * 1. preprocess: Initializes the preconditioner for a matrix, that is used as
90 : * the underlying matrix. Calling these method again with a
91 : * different matrix results in a different preconditioner.
92 : *
93 : * 2. step: Computes the new correction.
94 : *
95 : * 3. postprocess: Clean up (if needed)
96 : *
97 : * \tparam TAlgebra Type of Algebra
98 : */
99 : template <typename TAlgebra>
100 : class IPreconditioner :
101 : public ILinearIterator<typename TAlgebra::vector_type>,
102 : public DebugWritingObject<TAlgebra>
103 : {
104 : public:
105 : /// Algebra type
106 : typedef TAlgebra algebra_type;
107 :
108 : /// Vector type
109 : typedef typename TAlgebra::vector_type vector_type;
110 :
111 : /// Matrix type
112 : typedef typename TAlgebra::matrix_type matrix_type;
113 :
114 : /// Matrix Operator type
115 : typedef MatrixOperator<matrix_type, vector_type> matrix_operator_type;
116 : using DebugWritingObject<TAlgebra>::set_debug;
117 : protected:
118 : using ILinearIterator<vector_type>::damping;
119 :
120 : using DebugWritingObject<TAlgebra>::debug_writer;
121 : using DebugWritingObject<TAlgebra>::write_debug;
122 :
123 : public:
124 : /// default constructor
125 0 : IPreconditioner() :
126 0 : m_spDefectOperator(NULL), m_spApproxOperator(NULL), m_bInit(false), m_bOtherApproxOperator(false)
127 : {};
128 :
129 : /// constructor setting debug writer
130 : IPreconditioner(SmartPtr<IDebugWriter<algebra_type> > spDebugWriter) :
131 : DebugWritingObject<TAlgebra>(spDebugWriter),
132 : m_spDefectOperator(NULL), m_spApproxOperator(NULL), m_bInit(false), m_bOtherApproxOperator(false)
133 : {};
134 :
135 : /// clone constructor
136 0 : IPreconditioner( const IPreconditioner<TAlgebra> &parent ) :
137 : ILinearIterator<vector_type>(parent),
138 : DebugWritingObject<TAlgebra>(parent),
139 0 : m_spDefectOperator(NULL), m_spApproxOperator(NULL), m_bInit(false), m_bOtherApproxOperator(false)
140 : {
141 0 : }
142 : protected:
143 : /// returns the name of iterator
144 : /**
145 : * This method returns the name of the iterator operator. This function is
146 : * typically needed, when the iterator operator is used inside of another
147 : * operator and some debug output should be printed
148 : *
149 : * \returns const char* name of inverse operator
150 : */
151 : virtual const char* name() const = 0;
152 :
153 : /// initializes the preconditioner
154 : /**
155 : * This method is used to initialize the preconditioner. Usually, here
156 : * are performed computationally expensive operations, that should only be
157 : * computed once for an underlying matrix (e.g. LU factorization), while
158 : * the preconditioner will by applied (using 'step'-method) several times.
159 : *
160 : * \param[in] mat underlying matrix (i.e. L in L*u = f)
161 : * \returns bool success flag
162 : */
163 : virtual bool preprocess(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp) = 0;
164 :
165 : /// computes a new correction c = B*d
166 : /**
167 : * This method computes a new correction c = B*d. It can only be called,
168 : * when the preprocess has been done.
169 : *
170 : * \param[in] mat underlying matrix (i.e. L in L*u = f)
171 : * \param[out] c correction
172 : * \param[in] d defect
173 : * \returns bool success flag
174 : */
175 : virtual bool step(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp, vector_type& c, const vector_type& d) = 0;
176 :
177 : /// cleans the operator
178 : virtual bool postprocess() = 0;
179 :
180 : public:
181 : /// implements the ILinearIterator-interface for matrix based preconditioner
182 : /**
183 : * This method implements the ILinearIterator interface. It check if the
184 : * passed linear operator is matrix based (otherwise this preconditioner can
185 : * not be used for the linear operator). Then the request is forwarded to
186 : * the implementation of matrix based operators.
187 : *
188 : * \param[in] J linear operator
189 : * \param[in] u linearization point
190 : * \returns bool success flag
191 : */
192 0 : virtual bool init(SmartPtr<ILinearOperator<vector_type> > J,
193 : const vector_type& u)
194 : {
195 : // cast to matrix based operator
196 0 : SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp =
197 : J.template cast_dynamic<MatrixOperator<matrix_type, vector_type> >();
198 :
199 : // Check that matrix if of correct type
200 0 : if(pOp.invalid())
201 0 : UG_THROW(name() << "::init': Passed Operator is "
202 : "not based on matrix. This Preconditioner can only "
203 : "handle matrix-based operators.");
204 :
205 : // forward request to matrix based implementation
206 0 : return init(pOp);
207 : }
208 :
209 : /// implements the ILinearIterator-interface for matrix based preconditioner
210 : /**
211 : * This method implements the ILinearIterator interface. It check if the
212 : * passed linear operator is matrix based (otherwise this preconditioner can
213 : * not be used for the linear operator). Then the request is forwarded to
214 : * the implementation of matrix based operators.
215 : *
216 : * \param[in] L linear operator
217 : * \returns bool success flag
218 : */
219 0 : bool init(SmartPtr<ILinearOperator<vector_type> > L)
220 : {
221 : // Remember operator
222 0 : m_spDefectOperator = L;
223 0 : if(m_bOtherApproxOperator) return true;
224 :
225 : // cast to matrix based operator
226 0 : SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp =
227 : L.template cast_dynamic<MatrixOperator<matrix_type, vector_type> >();
228 :
229 : // Check that matrix if of correct type
230 0 : if(pOp.invalid())
231 0 : UG_THROW(name() << "::init': Passed Operator is "
232 : "not based on matrix. This Preconditioner can only "
233 : "handle matrix-based operators.");
234 :
235 : // forward request to matrix based implementation
236 0 : return init(pOp);
237 :
238 : }
239 :
240 :
241 : /// initializes the preconditioner for a matrix based operator
242 : /**
243 : * This method initializes the preconditioner for matrix based operators.
244 : * It performs some default checks and then forwards internally the
245 : * initialization to the (virtual) 'preprocess'-method
246 : *
247 : * \param[in] Op matrix based operator
248 : * \returns bool success flag
249 : */
250 0 : bool init(SmartPtr<MatrixOperator<matrix_type, vector_type> > Op)
251 : {
252 : // Remember operator
253 0 : m_spApproxOperator = Op;
254 0 : m_spDefectOperator = Op;
255 :
256 : // Check that matrix exists
257 0 : if(m_spApproxOperator.invalid())
258 0 : UG_THROW(name() << "::init': Passed Operator is invalid.");
259 :
260 : // Preprocess
261 0 : if(!preprocess(m_spApproxOperator))
262 : {
263 0 : UG_LOG("ERROR in '"<<name()<<"::init': Preprocess failed.\n");
264 0 : return false;
265 : }
266 :
267 : // Remember, that operator has been initialized
268 0 : m_bInit = true;
269 :
270 : // we're done
271 0 : return true;
272 : }
273 :
274 : /// compute new correction c = B*d
275 : /**
276 : * This method implements the virtual method of the ILinearIterator-interface.
277 : * Basically, besides some common checks the request is forwarded to the
278 : * (virtual) 'step'-method.
279 : *
280 : * \param[out] c correction
281 : * \param[in] d defect
282 : * \returns bool success flag
283 : */
284 0 : virtual bool apply(vector_type& c, const vector_type& d)
285 : {
286 : // Check that operator is initialized
287 0 : if(!m_bInit)
288 : {
289 0 : UG_LOG("ERROR in '"<<name()<<"::apply': Iterator not initialized.\n");
290 0 : return false;
291 : }
292 :
293 : // Check parallel status
294 : #ifdef UG_PARALLEL
295 : if(!d.has_storage_type(PST_ADDITIVE))
296 : UG_THROW(name() << "::apply: Wrong parallel "
297 : "storage format. Defect must be additive.");
298 : #endif
299 :
300 : // Check sizes
301 0 : THROW_IF_NOT_EQUAL_4(c.size(), d.size(),
302 : m_spApproxOperator->num_rows(), m_spApproxOperator->num_cols());
303 :
304 : // apply iterator: c = B*d
305 0 : if(!step(m_spApproxOperator, c, d))
306 : {
307 0 : UG_LOG("ERROR in '"<<name()<<"::apply': Step Routine failed.\n");
308 0 : return false;
309 : }
310 :
311 : // apply scaling
312 0 : const number kappa = damping()->damping(c, d, m_spApproxOperator);
313 0 : if(kappa != 1.0){
314 : c *= kappa;
315 : }
316 :
317 : // Correction is always consistent
318 : #ifdef UG_PARALLEL
319 : if(!c.change_storage_type(PST_CONSISTENT))
320 : UG_THROW(name() << "::apply': Cannot change "
321 : "parallel storage type of correction to consistent.");
322 : #endif
323 :
324 : // we're done
325 : return true;
326 : }
327 :
328 : /// compute new correction c = B*d and update defect d:= d - L*c
329 : /**
330 : * This method implements the virtual method of the ILinearIterator-interface.
331 : * Basically, the request is forwarded to the 'apply'-method and then the
332 : * update of the defect is computed afterwards.
333 : *
334 : * \param[out] c correction
335 : * \param[in, out] d defect on entry, updated defect on exit
336 : * \returns bool success flag
337 : */
338 0 : virtual bool apply_update_defect(vector_type& c, vector_type& d)
339 : {
340 : // compute new correction
341 0 : if(!apply(c, d)) return false;
342 :
343 : // update defect d := d - A*c
344 0 : m_spDefectOperator->apply_sub(d, c);
345 :
346 : // we're done
347 0 : return true;
348 : }
349 :
350 0 : virtual void set_approximation(SmartPtr<MatrixOperator<matrix_type,vector_type> > approx)
351 : {
352 0 : UG_COND_THROW(!approx.valid(), "");
353 0 : m_spApproxOperator = approx;
354 0 : if(!preprocess(m_spApproxOperator))
355 : {
356 0 : UG_THROW("ERROR in '"<<name()<<"::init': Preprocess failed.\n");
357 : return;
358 : }
359 0 : m_bOtherApproxOperator = true;
360 : }
361 :
362 : /// virtual destructor
363 0 : virtual ~IPreconditioner() {};
364 :
365 : /// underlying matrix based operator for calculation of defect
366 : SmartPtr<MatrixOperator<matrix_type, vector_type> > defect_operator()
367 : {
368 : return m_spDefectOperator;
369 : }
370 :
371 : /// underlying matrix based operator used for the preconditioner
372 : SmartPtr<MatrixOperator<matrix_type, vector_type> > approx_operator()
373 : {
374 : return m_spApproxOperator;
375 : }
376 :
377 : protected:
378 : /// underlying matrix based operator for calculation of defect
379 : SmartPtr<ILinearOperator<vector_type> > m_spDefectOperator;
380 : /// underlying matrix based operator used for the preconditioner
381 : SmartPtr<MatrixOperator<matrix_type, vector_type> > m_spApproxOperator;
382 :
383 :
384 : /// init flag indicating if init has been called
385 : bool m_bInit;
386 :
387 : bool m_bOtherApproxOperator;
388 : };
389 :
390 :
391 :
392 : } // end namespace ug
393 : #endif /* __H__LIB_ALGEBRA__OPERATOR__INTERFACE__PRECONDITIONER__ */
|
