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__LINEAR_SOLVER__
34 : #define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__LINEAR_SOLVER__
35 : #include <iostream>
36 : #include <string>
37 :
38 : #include "lib_algebra/operator/interface/preconditioned_linear_operator_inverse.h"
39 : #include "lib_algebra/operator/interface/linear_solver_profiling.h"
40 : #include "lib_algebra/operator/interface/pprocess.h"
41 : #ifdef UG_PARALLEL
42 : #include "lib_algebra/parallelization/parallelization.h"
43 : #endif
44 :
45 : namespace ug{
46 :
47 : /// linear solver using abstract preconditioner interface
48 : /**
49 : * This class is a linear iterating scheme, that uses any implementation
50 : * of the ILinearIterator interface to precondition the iteration.
51 : *
52 : * \tparam TAlgebra algebra type
53 : */
54 : template <typename TVector>
55 0 : class LinearSolver
56 : : public IPreconditionedLinearOperatorInverse<TVector>
57 : {
58 : public:
59 : /// Vector type
60 : typedef TVector vector_type;
61 :
62 : /// Base type
63 : typedef IPreconditionedLinearOperatorInverse<vector_type> base_type;
64 :
65 : /// constructors
66 0 : LinearSolver() : base_type() {}
67 :
68 0 : LinearSolver(SmartPtr<ILinearIterator<vector_type,vector_type> > spPrecond)
69 0 : : base_type ( spPrecond ) {}
70 :
71 0 : LinearSolver(SmartPtr<ILinearIterator<vector_type,vector_type> > spPrecond, SmartPtr<IConvergenceCheck<vector_type> > spConvCheck)
72 0 : : base_type ( spPrecond, spConvCheck) {}
73 :
74 : protected:
75 : using base_type::convergence_check;
76 : using base_type::linear_operator;
77 : using base_type::preconditioner;
78 : using base_type::write_debug;
79 :
80 : public:
81 : /// returns the name of the solver
82 0 : virtual const char* name() const {return "Iterative Linear Solver";}
83 :
84 : /// returns if parallel solving is supported
85 0 : virtual bool supports_parallel() const
86 : {
87 0 : if(preconditioner().valid())
88 0 : return preconditioner()->supports_parallel();
89 : else return true;
90 : }
91 :
92 : /**
93 : * Compute a correction c := B*d using one iterative step
94 : * Internally the defect is updated d := d - A*c = b - A*(x+c)
95 : * @param c
96 : * @param d
97 : */
98 0 : bool compute_correction(vector_type &c, vector_type &d)
99 : {
100 0 : if(preconditioner().valid()) {
101 : LS_PROFILE_BEGIN(LS_ApplyPrecond);
102 :
103 0 : if(!preconditioner()->apply_update_defect(c, d))
104 : {
105 : UG_LOG("ERROR in 'LinearSolver': Could not apply preconditioner. Aborting.\n");
106 0 : return false;
107 : }
108 : LS_PROFILE_END(LS_ApplyPrecond);
109 : }
110 : return true;
111 : }
112 :
113 : /// solves the system and returns the last defect
114 0 : virtual bool apply_return_defect(vector_type& x, vector_type& b)
115 : {
116 :
117 : LS_PROFILE_BEGIN(LS_ApplyReturnDefect);
118 :
119 : #ifdef UG_PARALLEL
120 : if(!b.has_storage_type(PST_ADDITIVE) || !x.has_storage_type(PST_CONSISTENT))
121 : UG_THROW("LinearSolver::apply: Inadequate parallel storage format of Vectors: "
122 : << b.get_storage_type() << " for b (expected " << PST_ADDITIVE << "), "
123 : << x.get_storage_type() << " for x (expected " << PST_CONSISTENT << ")");
124 : #endif
125 :
126 : // debug output
127 0 : if(this->vector_debug_writer_valid())
128 0 : write_debug(b, std::string("LS_RHS") + ".vec");
129 :
130 : // rename b as d (for convenience)
131 : vector_type& d = b;
132 :
133 : // build defect: d := b - J*x
134 : LS_PROFILE_BEGIN(LS_BuildDefect);
135 0 : linear_operator()->apply_sub(d, x);
136 : LS_PROFILE_END(LS_BuildDefect);
137 :
138 : // create correction
139 : LS_PROFILE_BEGIN(LS_CreateCorrection);
140 0 : SmartPtr<vector_type> spC = x.clone_without_values();
141 : vector_type& c = *spC;
142 : #ifdef UG_PARALLEL
143 : // this is ok if clone_without_values() inits with zeros
144 : c.set_storage_type(PST_CONSISTENT);
145 : #endif
146 : LS_PROFILE_END(LS_CreateCorrection);
147 :
148 : LS_PROFILE_BEGIN(LS_ComputeStartDefect);
149 0 : prepare_conv_check();
150 0 : convergence_check()->start(d);
151 : LS_PROFILE_END(LS_ComputeStartDefect);
152 :
153 : int loopCnt = 0;
154 0 : write_debugXCD(x, c, d, loopCnt, false);
155 :
156 : // Iteration loop
157 0 : while(!convergence_check()->iteration_ended())
158 : {
159 0 : enter_precond_debug_section(loopCnt);
160 0 : if( !compute_correction(c, d) )
161 : {
162 0 : this->leave_vector_debug_writer_section();
163 0 : return false;
164 : }
165 0 : this->leave_vector_debug_writer_section();
166 :
167 : // post-process the correction
168 0 : m_corr_post_process.apply (c);
169 :
170 : // add correction to solution: x += c
171 : LS_PROFILE_BEGIN(LS_AddCorrection);
172 0 : x += c;
173 : LS_PROFILE_END(LS_AddCorrection);
174 :
175 0 : write_debugXCD(x, c, d, ++loopCnt, true);
176 :
177 : // compute norm of new defect (in parallel)
178 : LS_PROFILE_BEGIN(LS_ComputeNewDefectNorm);
179 0 : convergence_check()->update(d);
180 : LS_PROFILE_END(LS_ComputeNewDefectNorm);
181 : }
182 :
183 : // write some information when ending the iteration
184 0 : if(!convergence_check()->post())
185 : {
186 : UG_LOG("ERROR in 'LinearSolver::apply': post-convergence-check "
187 : "signaled failure. Aborting.\n");
188 : return false;
189 : }
190 :
191 : // end profiling of whole function
192 : LS_PROFILE_END(LS_ApplyReturnDefect);
193 :
194 : // we're done
195 : return true;
196 : }
197 :
198 : /// adds a post-process for the iterates
199 0 : void add_postprocess_corr (SmartPtr<IPProcessVector<vector_type> > p)
200 : {
201 0 : m_corr_post_process.add (p);
202 0 : }
203 :
204 : /// removes a post-process for the iterates
205 0 : void remove_postprocess_corr (SmartPtr<IPProcessVector<vector_type> > p)
206 : {
207 0 : m_corr_post_process.remove (p);
208 0 : }
209 :
210 : protected:
211 : /// prepares the convergence check output
212 0 : void prepare_conv_check()
213 : {
214 0 : convergence_check()->set_name(name());
215 0 : convergence_check()->set_symbol('%');
216 0 : if(preconditioner().valid())
217 : {
218 : std::string s;
219 0 : if(preconditioner().valid())
220 0 : s = std::string(" (Precond: ") + preconditioner()->name() + ")";
221 : else
222 : s = " (No Preconditioner) ";
223 0 : convergence_check()->set_info(s);
224 : }
225 0 : }
226 :
227 : /// debugger output: solution, correction, defect
228 0 : void write_debugXCD(vector_type &x, vector_type &c, vector_type &d, int loopCnt, bool bWriteC)
229 : {
230 0 : if(!this->vector_debug_writer_valid()) return;
231 : char ext[20]; snprintf(ext, sizeof(ext),"_iter%03d", loopCnt);
232 0 : write_debug(d, std::string("LS_Defect") + ext + ".vec");
233 0 : if(bWriteC) write_debug(c, std::string("LS_Correction") + ext + ".vec");
234 0 : write_debug(x, std::string("LS_Solution") + ext + ".vec");
235 : }
236 :
237 : /// debugger section for the preconditioner
238 0 : void enter_precond_debug_section(int loopCnt)
239 : {
240 0 : if(!this->vector_debug_writer_valid()) return;
241 : char ext[20]; snprintf(ext, sizeof(ext),"_iter%03d", loopCnt);
242 0 : this->enter_vector_debug_writer_section(std::string("LS_Precond_") + ext);
243 : }
244 :
245 : protected:
246 : /// postprocessor for the correction in the iterations
247 : PProcessChain<vector_type> m_corr_post_process;
248 : };
249 :
250 : } // end namespace ug
251 :
252 : #endif /* __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__LINEAR_SOLVER__ */
|
