ug4
linear_solver.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__LINEAR_SOLVER__
34 #define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__LINEAR_SOLVER__
35 #include <iostream>
36 #include <string>
37 
41 #ifdef UG_PARALLEL
43 #endif
44 
45 namespace ug{
46 
48 
54 template <typename TVector>
56  : public IPreconditionedLinearOperatorInverse<TVector>
57 {
58  public:
60  typedef TVector vector_type;
61 
64 
67 
69  : base_type ( spPrecond ) {}
70 
72  : base_type ( spPrecond, spConvCheck) {}
73 
74  protected:
79 
80  public:
82  virtual const char* name() const {return "Iterative Linear Solver";}
83 
85  virtual bool supports_parallel() const
86  {
87  if(preconditioner().valid())
88  return preconditioner()->supports_parallel();
89  else return true;
90  }
91 
99  {
100  if(preconditioner().valid()) {
101  LS_PROFILE_BEGIN(LS_ApplyPrecond);
102 
103  if(!preconditioner()->apply_update_defect(c, d))
104  {
105  UG_LOG("ERROR in 'LinearSolver': Could not apply preconditioner. Aborting.\n");
106  return false;
107  }
108  LS_PROFILE_END(LS_ApplyPrecond);
109  }
110  return true;
111  }
112 
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  if(this->vector_debug_writer_valid())
128  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  linear_operator()->apply_sub(d, x);
136  LS_PROFILE_END(LS_BuildDefect);
137 
138  // create correction
139  LS_PROFILE_BEGIN(LS_CreateCorrection);
140  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);
150  convergence_check()->start(d);
151  LS_PROFILE_END(LS_ComputeStartDefect);
152 
153  int loopCnt = 0;
154  write_debugXCD(x, c, d, loopCnt, false);
155 
156  // Iteration loop
157  while(!convergence_check()->iteration_ended())
158  {
160  if( !compute_correction(c, d) )
161  {
163  return false;
164  }
166 
167  // post-process the correction
169 
170  // add correction to solution: x += c
171  LS_PROFILE_BEGIN(LS_AddCorrection);
172  x += c;
173  LS_PROFILE_END(LS_AddCorrection);
174 
175  write_debugXCD(x, c, d, ++loopCnt, true);
176 
177  // compute norm of new defect (in parallel)
178  LS_PROFILE_BEGIN(LS_ComputeNewDefectNorm);
179  convergence_check()->update(d);
180  LS_PROFILE_END(LS_ComputeNewDefectNorm);
181  }
182 
183  // write some information when ending the iteration
184  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 
200  {
202  }
203 
206  {
208  }
209 
210  protected:
213  {
214  convergence_check()->set_name(name());
215  convergence_check()->set_symbol('%');
216  if(preconditioner().valid())
217  {
218  std::string s;
219  if(preconditioner().valid())
220  s = std::string(" (Precond: ") + preconditioner()->name() + ")";
221  else
222  s = " (No Preconditioner) ";
223  convergence_check()->set_info(s);
224  }
225  }
226 
228  void write_debugXCD(vector_type &x, vector_type &c, vector_type &d, int loopCnt, bool bWriteC)
229  {
230  if(!this->vector_debug_writer_valid()) return;
231  char ext[20]; snprintf(ext, sizeof(ext),"_iter%03d", loopCnt);
232  write_debug(d, std::string("LS_Defect") + ext + ".vec");
233  if(bWriteC) write_debug(c, std::string("LS_Correction") + ext + ".vec");
234  write_debug(x, std::string("LS_Solution") + ext + ".vec");
235  }
236 
238  void enter_precond_debug_section(int loopCnt)
239  {
240  if(!this->vector_debug_writer_valid()) return;
241  char ext[20]; snprintf(ext, sizeof(ext),"_iter%03d", loopCnt);
242  this->enter_vector_debug_writer_section(std::string("LS_Precond_") + ext);
243  }
244 
245  protected:
248 };
249 
250 } // end namespace ug
251 
252 #endif /* __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__LINEAR_SOLVER__ */
parameterString p
parameterString s
Definition: smart_pointer.h:108
Definition: convergence_check.h:72
describes a linear iterator
Definition: linear_iterator.h:81
ConstSmartPtr< IConvergenceCheck< X > > convergence_check() const
returns the convergence check
Definition: linear_operator_inverse.h:202
SmartPtr< ILinearOperator< Y, X > > linear_operator()
returns the current Operator this Inverse Operator is initialized for
Definition: linear_operator_inverse.h:227
virtual bool apply_update_defect(Y &u, X &f)
Definition: linear_operator_inverse.h:190
interface for pre- and postprocess functions
Definition: pprocess.h:59
describes an inverse linear mapping X->X
Definition: preconditioned_linear_operator_inverse.h:63
SmartPtr< ILinearIterator< X, X > > preconditioner()
Definition: preconditioned_linear_operator_inverse.h:118
linear solver using abstract preconditioner interface
Definition: linear_solver.h:57
TVector vector_type
Vector type.
Definition: linear_solver.h:60
bool compute_correction(vector_type &c, vector_type &d)
Definition: linear_solver.h:98
virtual const char * name() const
returns the name of the solver
Definition: linear_solver.h:82
void write_debugXCD(vector_type &x, vector_type &c, vector_type &d, int loopCnt, bool bWriteC)
debugger output: solution, correction, defect
Definition: linear_solver.h:228
IPreconditionedLinearOperatorInverse< vector_type > base_type
Base type.
Definition: linear_solver.h:63
LinearSolver(SmartPtr< ILinearIterator< vector_type, vector_type > > spPrecond)
Definition: linear_solver.h:68
virtual bool supports_parallel() const
returns if parallel solving is supported
Definition: linear_solver.h:85
LinearSolver(SmartPtr< ILinearIterator< vector_type, vector_type > > spPrecond, SmartPtr< IConvergenceCheck< vector_type > > spConvCheck)
Definition: linear_solver.h:71
virtual bool apply_return_defect(vector_type &x, vector_type &b)
solves the system and returns the last defect
Definition: linear_solver.h:114
SmartPtr< ILinearIterator< X, X > > preconditioner()
Definition: preconditioned_linear_operator_inverse.h:118
void add_postprocess_corr(SmartPtr< IPProcessVector< vector_type > > p)
adds a post-process for the iterates
Definition: linear_solver.h:199
PProcessChain< vector_type > m_corr_post_process
postprocessor for the correction in the iterations
Definition: linear_solver.h:247
void enter_precond_debug_section(int loopCnt)
debugger section for the preconditioner
Definition: linear_solver.h:238
void remove_postprocess_corr(SmartPtr< IPProcessVector< vector_type > > p)
removes a post-process for the iterates
Definition: linear_solver.h:205
LinearSolver()
constructors
Definition: linear_solver.h:66
void prepare_conv_check()
prepares the convergence check output
Definition: linear_solver.h:212
void add(SmartPtr< p_process_type > p)
adds an operation at the end of the chain
Definition: pprocess.h:96
void apply(vector_type &v)
performs all the operations
Definition: pprocess.h:86
void remove(SmartPtr< p_process_type > p)
removes an operation from the chain
Definition: pprocess.h:108
bool vector_debug_writer_valid() const
returns true if the debug writer is set
Definition: debug_writer.h:290
void enter_vector_debug_writer_section(std::string secDir)
enters a debugging section
Definition: debug_writer.h:331
void leave_vector_debug_writer_section()
leaves a debugging section
Definition: debug_writer.h:342
void write_debug(const vector_type &vec, const char *filename)
writing debug output for a vector (if debug writer set)
Definition: debug_writer.h:293
TVector vector_type
type of vector
Definition: debug_writer.h:269
@ PST_CONSISTENT
Definition: parallel_storage_type.h:68
@ PST_ADDITIVE
Definition: parallel_storage_type.h:69
#define UG_THROW(msg)
Definition: error.h:57
#define UG_LOG(msg)
Definition: log.h:367
#define LS_PROFILE_BEGIN(name)
Definition: linear_solver_profiling.h:40
#define LS_PROFILE_END(name)
Definition: linear_solver_profiling.h:41
the ug namespace