Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
solver_richardson.h
1 // ---------------------------------------------------------------------
2 // @f$Id: solver_richardson.h 31349 2013-10-20 19:07:06Z maier @f$
3 //
4 // Copyright (C) 1999 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__solver_richardson_h
18 #define __deal2__solver_richardson_h
19 
20 
21 #include <deal.II/base/config.h>
22 #include <deal.II/base/logstream.h>
23 #include <deal.II/lac/solver.h>
24 #include <deal.II/lac/solver_control.h>
25 #include <deal.II/base/subscriptor.h>
26 
28 
31 
54 template <class VECTOR = Vector<double> >
55 class SolverRichardson : public Solver<VECTOR>
56 {
57 public:
64  {
70  AdditionalData (const double omega = 1,
71  const bool use_preconditioned_residual = false);
72 
76  double omega;
77 
82 
83  };
84 
90  const AdditionalData &data=AdditionalData());
91 
98  const AdditionalData &data=AdditionalData());
99 
103  virtual ~SolverRichardson ();
104 
109  template<class MATRIX, class PRECONDITIONER>
110  void
111  solve (const MATRIX &A,
112  VECTOR &x,
113  const VECTOR &b,
114  const PRECONDITIONER &precondition);
115 
119  template<class MATRIX, class PRECONDITIONER>
120  void
121  Tsolve (const MATRIX &A,
122  VECTOR &x,
123  const VECTOR &b,
124  const PRECONDITIONER &precondition);
125 
130  void set_omega (const double om=1.);
131 
141  virtual void print_vectors(const unsigned int step,
142  const VECTOR &x,
143  const VECTOR &r,
144  const VECTOR &d) const;
145 
146 protected:
151  virtual typename VECTOR::value_type criterion();
152 
159  VECTOR *Vr;
169  VECTOR *Vd;
170 
175 
186  typename VECTOR::value_type res;
187 };
188 
190 /*----------------- Implementation of the Richardson Method ------------------*/
191 
192 #ifndef DOXYGEN
193 
194 template <class VECTOR>
195 inline
197 AdditionalData (const double omega,
198  const bool use_preconditioned_residual)
199  :
200  omega(omega),
201  use_preconditioned_residual(use_preconditioned_residual)
202 {}
203 
204 
205 template <class VECTOR>
208  const AdditionalData &data)
209  :
210  Solver<VECTOR> (cn,mem),
211  additional_data(data)
212 {}
213 
214 
215 
216 template <class VECTOR>
218  const AdditionalData &data)
219  :
220  Solver<VECTOR> (cn),
221  additional_data(data)
222 {}
223 
224 
225 
226 template <class VECTOR>
228 {}
229 
230 
231 template <class VECTOR>
232 template <class MATRIX, class PRECONDITIONER>
233 void
235  VECTOR &x,
236  const VECTOR &b,
237  const PRECONDITIONER &precondition)
238 {
240 
241  // Memory allocation
242  Vr = this->memory.alloc();
243  VECTOR &r = *Vr;
244  r.reinit(x);
245  Vd = this->memory.alloc();
246  VECTOR &d = *Vd;
247  d.reinit(x);
248 
249  deallog.push("Richardson");
250 
251  try
252  {
253  // Main loop
254  for (int iter=0; conv==SolverControl::iterate; iter++)
255  {
256  // Do not use residual,
257  // but do it in 2 steps
258  A.vmult(r,x);
259  r.sadd(-1.,1.,b);
260  precondition.vmult(d,r);
261 
262  // The required norm of the
263  // (preconditioned)
264  // residual is computed in
265  // criterion() and stored
266  // in res.
267  conv = this->control().check (iter, criterion());
268 // conv = this->control().check (iter, std::sqrt(A.matrix_norm_square(r)));
269  if (conv != SolverControl::iterate)
270  break;
271 
272  x.add(additional_data.omega,d);
273  print_vectors(iter,x,r,d);
274  }
275  }
276  catch (...)
277  {
278  this->memory.free(Vr);
279  this->memory.free(Vd);
280  deallog.pop();
281  throw;
282  }
283  // Deallocate Memory
284  this->memory.free(Vr);
285  this->memory.free(Vd);
286  deallog.pop();
287 
288  // in case of failure: throw exception
289  if (this->control().last_check() != SolverControl::success)
290  AssertThrow(false, SolverControl::NoConvergence (this->control().last_step(),
291  this->control().last_value()));
292  // otherwise exit as normal
293 }
294 
295 
296 template <class VECTOR>
297 template <class MATRIX, class PRECONDITIONER>
298 void
300  VECTOR &x,
301  const VECTOR &b,
302  const PRECONDITIONER &precondition)
303 {
305 
306  // Memory allocation
307  Vr = this->memory.alloc();
308  VECTOR &r = *Vr;
309  r.reinit(x);
310  Vd =this-> memory.alloc();
311  VECTOR &d = *Vd;
312  d.reinit(x);
313 
314  deallog.push("RichardsonT");
315 
316  try
317  {
318  // Main loop
319  for (int iter=0; conv==SolverControl::iterate; iter++)
320  {
321  // Do not use Tresidual,
322  // but do it in 2 steps
323  A.Tvmult(r,x);
324  r.sadd(-1.,1.,b);
325  precondition.Tvmult(d,r);
326 
327  conv = this->control().check (iter, criterion());
328  if (conv != SolverControl::iterate)
329  break;
330 
331  x.add(additional_data.omega,d);
332  print_vectors(iter,x,r,d);
333  }
334  }
335  catch (...)
336  {
337  this->memory.free(Vr);
338  this->memory.free(Vd);
339  deallog.pop();
340  throw;
341  }
342 
343  // Deallocate Memory
344  this->memory.free(Vr);
345  this->memory.free(Vd);
346  deallog.pop();
347  // in case of failure: throw exception
348  if (this->control().last_check() != SolverControl::success)
349  AssertThrow(false, SolverControl::NoConvergence (this->control().last_step(),
350  this->control().last_value()));
351  // otherwise exit as normal
352 }
353 
354 
355 template <class VECTOR>
356 void
357 SolverRichardson<VECTOR>::print_vectors(const unsigned int,
358  const VECTOR &,
359  const VECTOR &,
360  const VECTOR &) const
361 {}
362 
363 
364 
365 template <class VECTOR>
366 inline typename VECTOR::value_type
368 {
369  if (!additional_data.use_preconditioned_residual)
370  res = Vr->l2_norm();
371  else
372  res = Vd->l2_norm();
373  return res;
374 }
375 
376 
377 template <class VECTOR>
378 inline void
379 SolverRichardson<VECTOR>::set_omega (const double om)
380 {
381  additional_data.omega=om;
382 }
383 
384 #endif // DOXYGEN
385 
386 DEAL_II_NAMESPACE_CLOSE
387 
388 #endif
VECTOR::value_type res
void Tvmult(VECTOR &u, const VECTOR &v) const
void pop()
SolverRichardson(SolverControl &cn, VectorMemory< VECTOR > &mem, const AdditionalData &data=AdditionalData())
void vmult(VECTOR &u, const VECTOR &v) const
AdditionalData(const double omega=1, const bool use_preconditioned_residual=false)
Continue iteration.
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
void solve(const MATRIX &A, VECTOR &x, const VECTOR &b, const PRECONDITIONER &precondition)
void Tsolve(const MATRIX &A, VECTOR &x, const VECTOR &b, const PRECONDITIONER &precondition)
AdditionalData additional_data
Definition: solver.h:147
void push(const std::string &text)
virtual ~SolverRichardson()
void set_omega(const double om=1.)
virtual VECTOR::value_type criterion()
virtual void print_vectors(const unsigned int step, const VECTOR &x, const VECTOR &r, const VECTOR &d) const
Stop iteration, goal reached.