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__LINEAR_OPERATOR_INVERSE__
34 : #define __H__LIB_ALGEBRA__OPERATOR__INTERFACE__LINEAR_OPERATOR_INVERSE__
35 :
36 : #include "linear_operator.h"
37 : #include "linear_iterator.h"
38 : #include "lib_algebra/operator/convergence_check.h"
39 : #include "common/error.h"
40 : #include "common/util/smart_pointer.h"
41 :
42 : namespace ug{
43 :
44 :
45 : ///////////////////////////////////////////////////////////////////////////////
46 : // Inverse of a Linear Operator
47 : ///////////////////////////////////////////////////////////////////////////////
48 :
49 : /// describes an inverse linear mapping X->Y
50 : /**
51 : * This class is the base class for the inversion of linear operator given in
52 : * form of the ILinearOperator interface class. Given a operator L, the basic
53 : * usage of this class is to invert this operator, i.e. to compute the solution
54 : * u of
55 : *
56 : * L*u = f i.e. u := L^{-1} f
57 : *
58 : * This application has been split up into three steps:
59 : *
60 : * 1. init(): This method initializes the inverse operator. The inverse operator
61 : * is initialized the way that, its application will be the inverse
62 : * application of the operator L passed in by this function. The
63 : * prepare method can only be called, when this method has been
64 : * called once.
65 : *
66 : * 3. apply(): This method performs the inversion. Before this method is called
67 : * the init and prepare methods have to be called.
68 : *
69 : * This splitting has been made, since initialization and preparation may be
70 : * computationally expansive. Thus, the user of this class has the choice
71 : * when to call this initialization/preparation. E.g. when the operator is
72 : * applied several times on the same vectors, those have only to be prepared
73 : * once and the init of the operator is only needed once.
74 : *
75 : * \tparam X domain space
76 : * \tparam Y range space
77 : */
78 : template <typename X, typename Y = X>
79 : class ILinearOperatorInverse : public ILinearIterator<X,Y>
80 : {
81 : public:
82 : /// Domain space
83 : typedef X domain_function_type;
84 :
85 : /// Range space
86 : typedef Y codomain_function_type;
87 :
88 : public:
89 : /// constructor setting convergence check to (100, 1e-12, 1e-12, true)
90 0 : ILinearOperatorInverse()
91 : : m_spLinearOperator(NULL),
92 0 : m_spConvCheck(new StdConvCheck<X>(100, 1e-12, 1e-12, true))
93 0 : {}
94 :
95 : /// Default constructor
96 0 : ILinearOperatorInverse(SmartPtr<IConvergenceCheck<X> > spConvCheck)
97 : : m_spLinearOperator(NULL),
98 0 : m_spConvCheck(spConvCheck)
99 0 : {}
100 :
101 : /// virtual destructor
102 0 : virtual ~ILinearOperatorInverse() {};
103 :
104 : /// returns the name of the operator inverse
105 : /**
106 : * This method returns the name of the inverse operator. This function is
107 : * typically needed, when the inverse operator is used inside of another and
108 : * some debug output should be printed
109 : *
110 : * \returns const char* name of inverse operator
111 : */
112 : virtual const char* name() const = 0;
113 :
114 :
115 : /// returns information about configuration parameters
116 : /**
117 : * this should return necessary information about parameters and possibly
118 : * calling config_string of subcomponents.
119 : *
120 : * \returns std::string necessary information about configuration parameters
121 : */
122 0 : virtual std::string config_string() const { return name(); }
123 :
124 : /// returns if parallel solving is supported
125 : virtual bool supports_parallel() const = 0;
126 :
127 : /// initializes for the inverse for a linear operator
128 : /**
129 : * This method passes the operator L that is inverted by this operator. In
130 : * addition some preparation step can be made.
131 : *
132 : * \param[in] L linear operator to invert
133 : * \returns bool success flag
134 : */
135 0 : virtual bool init(SmartPtr<ILinearOperator<Y,X> > L)
136 : {
137 : // remember operator
138 0 : m_spLinearOperator = L;
139 0 : return true;
140 : }
141 :
142 : /// initializes for the inverse for a linearized operator at linearization point u
143 : /**
144 : * This method passes the linear operator J(u) that should be inverted by
145 : * this operator. As second argument the linearization point is passed.
146 : * This is needed e.g. for the geometric multigrid method, that inverts
147 : * a linearized operator based on coarser grid operators, that have to be
148 : * initialized based on the linearization point.
149 : *
150 : * \param[in] J linearized operator to invert
151 : * \param[in] u linearization point
152 : * \returns bool success flag
153 : */
154 0 : virtual bool init(SmartPtr<ILinearOperator<Y,X> > J, const Y& u)
155 : {
156 : // remember operator
157 0 : m_spLinearOperator = J;
158 0 : return true;
159 : }
160 :
161 : /// applies inverse operator, i.e. returns u = A^{-1} f
162 : /**
163 : * This method applies the inverse operator, i.e. u = A^{-1} f. The
164 : * domain function f remains unchanged.
165 : * Note, that this method can always be implemented by creating a copy of
166 : * f and calling apply_return_defect with this copy.
167 : *
168 : * \param[in] f right-hand side
169 : * \param[out] u solution
170 : * \returns bool success flag
171 : */
172 : virtual bool apply(Y& u, const X& f) = 0;
173 :
174 : /// applies inverse operator, i.e. returns u = A^{-1} f and returns defect d := f - A*u
175 : /**
176 : * This method applies the inverse operator, i.e. u = A^{-1} f. The
177 : * domain function f is changed in the way, that the defect d := f - A*u
178 : * is returned in the function. This is always useful, when the inverting
179 : * algorithm can (or must) update the defect during computation (this is
180 : * e.g. the case for the geometric multigrid method).
181 : * Note, that this method can always be implemented by calling apply and
182 : * then computing d := f - A*u.
183 : *
184 : * \param[in,out] f right-hand side
185 : * \param[out] u solution
186 : * \returns bool success flag
187 : */
188 : virtual bool apply_return_defect(Y& u, X& f) = 0;
189 :
190 0 : virtual bool apply_update_defect(Y& u, X& f)
191 : {
192 0 : return apply_return_defect(u,f);
193 : }
194 :
195 0 : virtual SmartPtr<ILinearIterator<X,Y> > clone()
196 : {
197 0 : UG_THROW("No cloning implemented.");
198 : return SPNULL;
199 : }
200 :
201 : /// returns the convergence check
202 0 : ConstSmartPtr<IConvergenceCheck<X> > convergence_check() const {return m_spConvCheck;}
203 :
204 : /// returns the convergence check
205 : SmartPtr<IConvergenceCheck<X> > convergence_check() {return m_spConvCheck;}
206 :
207 : /// returns the current defect
208 0 : number defect() const {return convergence_check()->defect();}
209 :
210 : /// returns the current number of steps
211 0 : int step() const {return convergence_check()->step();}
212 :
213 : /// returns the current relative reduction
214 0 : number reduction() const {return convergence_check()->reduction();}
215 :
216 : /// returns the standard offset for output
217 0 : virtual int standard_offset() const {return 3;}
218 :
219 : /// set the convergence check
220 0 : void set_convergence_check(SmartPtr<IConvergenceCheck<X> > spConvCheck)
221 : {
222 0 : m_spConvCheck = spConvCheck;
223 0 : m_spConvCheck->set_offset(standard_offset());
224 0 : };
225 :
226 : /// returns the current Operator this Inverse Operator is initialized for
227 0 : SmartPtr<ILinearOperator<Y,X> > linear_operator()
228 : {
229 0 : if(m_spLinearOperator.invalid())
230 0 : UG_THROW(name() << ": Linear Operator that should be "
231 : "inverted has not been set.");
232 :
233 0 : return m_spLinearOperator;
234 : }
235 :
236 : protected:
237 : /// Operator that is inverted by this Inverse Operator
238 : SmartPtr<ILinearOperator<Y,X> > m_spLinearOperator;
239 :
240 : /// smart pointer holding the convergence check
241 : SmartPtr<IConvergenceCheck<X> > m_spConvCheck;
242 : };
243 :
244 : }
245 : #endif /* __H__LIB_ALGEBRA__OPERATOR__INTERFACE__LINEAR_OPERATOR_INVERSE__ */
|
