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__PRECONDITIONED_LINEAR_OPERATOR_INVERSE__
34 : #define __H__LIB_ALGEBRA__OPERATOR__INTERFACE__PRECONDITIONED_LINEAR_OPERATOR_INVERSE__
35 :
36 : #include "linear_operator_inverse.h"
37 : #include "linear_iterator.h"
38 : #include "lib_algebra/operator/debug_writer.h"
39 : #include "common/util/smart_pointer.h"
40 : #include "lib_algebra/operator/convergence_check.h"
41 : #include "common/log.h"
42 : #include "linear_solver_profiling.h"
43 :
44 : #undef DEBUG_FOR_AMG
45 :
46 : namespace ug{
47 : ///////////////////////////////////////////////////////////////////////////////
48 : // Inverse of a Linear Operator using a ILinearIterator as preconditioner
49 : ///////////////////////////////////////////////////////////////////////////////
50 :
51 : /// describes an inverse linear mapping X->X
52 : /**
53 : * This a useful derived class from ILinearOperatorInverse, that uses a
54 : * ILinearIterator in order to precondition the solution process. This is
55 : * used e.g. in LinearSolver, CG and BiCGStab.
56 : *
57 : * \tparam X domain and range space
58 : */
59 : template <typename X>
60 : class IPreconditionedLinearOperatorInverse
61 : : public ILinearOperatorInverse<X>,
62 : public VectorDebugWritingObject<X>
63 : {
64 : public:
65 : /// Domain space
66 : typedef X domain_function_type;
67 :
68 : /// Range space
69 : typedef X codomain_function_type;
70 :
71 : /// Base class
72 : typedef ILinearOperatorInverse<X,X> base_type;
73 :
74 : protected:
75 : using base_type::linear_operator;
76 : using base_type::m_spConvCheck;
77 :
78 : public:
79 : using VectorDebugWritingObject<X>::write_debug;
80 : using base_type::name;
81 : using base_type::apply_return_defect;
82 :
83 : public:
84 : /// Empty constructor
85 0 : IPreconditionedLinearOperatorInverse()
86 0 : : m_bRecompute(false), m_spPrecond(NULL)
87 : #ifdef DEBUG_FOR_AMG
88 : , m_amgDebug(0)
89 : #endif
90 : {}
91 :
92 : /// constructor setting the preconditioner
93 0 : IPreconditionedLinearOperatorInverse(SmartPtr<ILinearIterator<X,X> > spPrecond)
94 0 : : m_bRecompute(false), m_spPrecond(spPrecond)
95 : #ifdef DEBUG_FOR_AMG
96 : , m_amgDebug(0)
97 : #endif
98 0 : {}
99 :
100 : /// constructor setting the preconditioner
101 0 : IPreconditionedLinearOperatorInverse(SmartPtr<ILinearIterator<X,X> > spPrecond,
102 : SmartPtr<IConvergenceCheck<X> > spConvCheck)
103 : : base_type(spConvCheck),
104 0 : m_bRecompute(false), m_spPrecond(spPrecond)
105 : #ifdef DEBUG_FOR_AMG
106 : , m_amgDebug(0)
107 : #endif
108 0 : {}
109 :
110 : /// sets the preconditioner
111 0 : void set_preconditioner(SmartPtr<ILinearIterator<X, X> > spPrecond)
112 : {
113 0 : m_spPrecond = spPrecond;
114 0 : }
115 :
116 : /// returns the preconditioner
117 : /// \{
118 : SmartPtr<ILinearIterator<X, X> > preconditioner(){return m_spPrecond;}
119 : ConstSmartPtr<ILinearIterator<X, X> > preconditioner() const {return m_spPrecond;}
120 : /// \}
121 :
122 : /// initializes the solver for an operator
123 0 : virtual bool init(SmartPtr<ILinearOperator<X,X> > J, const X& u)
124 : {
125 0 : if(!base_type::init(J, u)) return false;
126 :
127 : LS_PROFILE_BEGIN(LS_InitPrecond);
128 0 : if(m_spPrecond.valid())
129 0 : if(!m_spPrecond->init(J, u))
130 0 : UG_THROW(name() << "::init: Cannot init Preconditioner "
131 : "Operator for Operator J.");
132 : LS_PROFILE_END(LS_InitPrecond);
133 :
134 : return true;
135 : }
136 :
137 : /// initializes the solver for an operator
138 0 : virtual bool init(SmartPtr<ILinearOperator<X,X> > L)
139 : {
140 0 : if(!base_type::init(L)) return false;
141 :
142 : LS_PROFILE_BEGIN(LS_InitPrecond);
143 0 : if(m_spPrecond.valid())
144 0 : if(!m_spPrecond->init(L))
145 0 : UG_THROW(name() <<"::prepare: Cannot init Preconditioner "
146 : "Operator for Operator L.");
147 : LS_PROFILE_END(LS_InitPrecond);
148 :
149 : return true;
150 : }
151 :
152 0 : virtual bool apply(X& x, const X& b)
153 : {
154 : // copy defect
155 0 : SmartPtr<X> spB = b.clone(); X& bTmp = *spB;
156 : // X bTmp; bTmp.resize(b.size()); bTmp = b;
157 :
158 : // solve on copy of defect
159 0 : bool bRes = apply_return_defect(x, bTmp);
160 :
161 : // write updated defect
162 0 : write_debug(bTmp, "LS_UpdatedDefectEnd.vec");
163 :
164 : // compute defect again, for debug purpose
165 0 : if(m_bRecompute)
166 : {
167 : // recompute defect
168 0 : bTmp = b; linear_operator()->apply_sub(bTmp, x);
169 0 : number norm = bTmp.norm();
170 :
171 : // print norm of recomputed defect
172 0 : UG_LOG("%%%% DEBUG "<<name()<<": (Re)computed defect has norm: "
173 : <<norm<<"\n");
174 :
175 : // write true end defect
176 0 : write_debug(bTmp, "LS_TrueDefectEnd.vec");
177 : }
178 :
179 : #ifdef DEBUG_FOR_AMG
180 : if (m_amgDebug>0)
181 : {
182 : // convergence post-check
183 : X myError(x.size());
184 : myError.set_random(-1.0, 1.0);
185 :
186 : bTmp.set(0.0);
187 :
188 : this->write_debug(myError, "AMGDebugPre");
189 : apply_return_defect(myError, bTmp);
190 : this->write_debug(myError, "AMGDebugPost");
191 : }
192 : #endif
193 :
194 :
195 : // return
196 0 : return bRes;
197 : }
198 :
199 : /// returns config information of convergence check and preconditioner
200 0 : std::string config_string_preconditioner_convergence_check() const
201 : {
202 0 : std::stringstream ss;
203 0 : ss << " Convergence Check: ";
204 0 : if(m_spConvCheck.valid()) ss << ConfigShift(m_spConvCheck->config_string()) << "\n";
205 0 : else ss << " NOT SET!\n";
206 0 : ss << " Preconditioner: ";
207 0 : if(m_spPrecond.valid()) ss << ConfigShift(m_spPrecond->config_string()) << "\n";
208 0 : else ss << " NOT SET!\n";
209 0 : return ss.str();
210 0 : }
211 :
212 : /// returns information about configuration parameters
213 0 : virtual std::string config_string() const
214 : {
215 0 : std::stringstream ss;
216 0 : ss << name() << "\n" << config_string_preconditioner_convergence_check();
217 0 : return ss.str();
218 0 : }
219 :
220 : /// for debug: computes norm again after whole calculation of apply
221 0 : void set_compute_fresh_defect_when_finished(bool bRecompute)
222 : {
223 0 : m_bRecompute = bRecompute;
224 0 : }
225 :
226 : protected:
227 : /// flag if fresh defect should be computed when finish for debug purpose
228 : bool m_bRecompute;
229 : /// Iterator used in the iterative scheme to compute the correction and update the defect
230 : SmartPtr<ILinearIterator<X,X> > m_spPrecond;
231 :
232 : #ifdef DEBUG_FOR_AMG
233 : public:
234 : void set_debug_amg(int b) {m_amgDebug = b;}
235 : int m_amgDebug;
236 : #endif
237 : };
238 :
239 : }
240 : #endif /* __H__LIB_ALGEBRA__OPERATOR__INTERFACE__PRECONDITIONED_LINEAR_OPERATOR_INVERSE__ */
|
