Reference documentation for deal.II version Git e3a3ec7800 2020-08-07 14:08:19 +0200
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
solver_cg.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_solver_cg_h
17 #define dealii_solver_cg_h
18 
19 
20 #include <deal.II/base/config.h>
21 
23 #include <deal.II/base/logstream.h>
25 
26 #include <deal.II/lac/solver.h>
29 
30 #include <cmath>
31 
33 
34 // forward declaration
35 #ifndef DOXYGEN
37 #endif
38 
39 
42 
94 template <typename VectorType = Vector<double>>
95 class SolverCG : public SolverBase<VectorType>
96 {
97 public:
102 
109  {};
110 
114  SolverCG(SolverControl & cn,
116  const AdditionalData & data = AdditionalData());
117 
122  SolverCG(SolverControl &cn, const AdditionalData &data = AdditionalData());
123 
127  virtual ~SolverCG() override = default;
128 
132  template <typename MatrixType, typename PreconditionerType>
133  void
134  solve(const MatrixType & A,
135  VectorType & x,
136  const VectorType & b,
137  const PreconditionerType &preconditioner);
138 
145  boost::signals2::connection
147  const std::function<void(typename VectorType::value_type,
148  typename VectorType::value_type)> &slot);
149 
156  boost::signals2::connection
157  connect_condition_number_slot(const std::function<void(double)> &slot,
158  const bool every_iteration = false);
159 
166  boost::signals2::connection
168  const std::function<void(const std::vector<double> &)> &slot,
169  const bool every_iteration = false);
170 
171 protected:
177  virtual void
178  print_vectors(const unsigned int step,
179  const VectorType & x,
180  const VectorType & r,
181  const VectorType & d) const;
182 
188  static void
190  const std::vector<typename VectorType::value_type> &diagonal,
191  const std::vector<typename VectorType::value_type> &offdiagonal,
192  const boost::signals2::signal<void(const std::vector<double> &)>
194  const boost::signals2::signal<void(double)> &cond_signal);
195 
200 
204  boost::signals2::signal<void(typename VectorType::value_type,
205  typename VectorType::value_type)>
207 
212  boost::signals2::signal<void(double)> condition_number_signal;
213 
218  boost::signals2::signal<void(double)> all_condition_numbers_signal;
219 
224  boost::signals2::signal<void(const std::vector<double> &)> eigenvalues_signal;
225 
230  boost::signals2::signal<void(const std::vector<double> &)>
232 };
233 
236 /*------------------------- Implementation ----------------------------*/
237 
238 #ifndef DOXYGEN
239 
240 template <typename VectorType>
243  const AdditionalData & data)
244  : SolverBase<VectorType>(cn, mem)
245  , additional_data(data)
246 {}
247 
248 
249 
250 template <typename VectorType>
253  , additional_data(data)
254 {}
255 
256 
257 
258 template <typename VectorType>
259 void
260 SolverCG<VectorType>::print_vectors(const unsigned int,
261  const VectorType &,
262  const VectorType &,
263  const VectorType &) const
264 {}
265 
266 
267 
268 template <typename VectorType>
269 inline void
271  const std::vector<typename VectorType::value_type> &diagonal,
272  const std::vector<typename VectorType::value_type> &offdiagonal,
273  const boost::signals2::signal<void(const std::vector<double> &)>
275  const boost::signals2::signal<void(double)> &cond_signal)
276 {
277  // Avoid computing eigenvalues unless they are needed.
278  if (!cond_signal.empty() || !eigenvalues_signal.empty())
279  {
281  true);
282  for (size_type i = 0; i < diagonal.size(); ++i)
283  {
284  T(i, i) = diagonal[i];
285  if (i < diagonal.size() - 1)
286  T(i, i + 1) = offdiagonal[i];
287  }
288  T.compute_eigenvalues();
289  // Need two eigenvalues to estimate the condition number.
290  if (diagonal.size() > 1)
291  {
292  auto condition_number = T.eigenvalue(T.n() - 1) / T.eigenvalue(0);
293  // Condition number is real valued and nonnegative; simply take
294  // the absolute value:
295  cond_signal(std::abs(condition_number));
296  }
297  // Avoid copying the eigenvalues of T to a vector unless a signal is
298  // connected.
299  if (!eigenvalues_signal.empty())
300  {
301  std::vector<double> eigenvalues(T.n());
302  for (unsigned int j = 0; j < T.n(); ++j)
303  {
304  // for a hermitian matrix, all eigenvalues are real-valued
305  // and non-negative, simply return the absolute value:
306  eigenvalues[j] = std::abs(T.eigenvalue(j));
307  }
309  }
310  }
311 }
312 
313 
314 
315 template <typename VectorType>
316 template <typename MatrixType, typename PreconditionerType>
317 void
318 SolverCG<VectorType>::solve(const MatrixType & A,
319  VectorType & x,
320  const VectorType & b,
321  const PreconditionerType &preconditioner)
322 {
323  using number = typename VectorType::value_type;
324 
326 
327  LogStream::Prefix prefix("cg");
328 
329  // Memory allocation
330  typename VectorMemory<VectorType>::Pointer g_pointer(this->memory);
331  typename VectorMemory<VectorType>::Pointer d_pointer(this->memory);
332  typename VectorMemory<VectorType>::Pointer h_pointer(this->memory);
333 
334  // define some aliases for simpler access
335  VectorType &g = *g_pointer;
336  VectorType &d = *d_pointer;
337  VectorType &h = *h_pointer;
338 
339  // Should we build the matrix for eigenvalue computations?
340  const bool do_eigenvalues =
342  !eigenvalues_signal.empty() || !all_eigenvalues_signal.empty();
343 
344  // vectors used for eigenvalue
345  // computations
346  std::vector<typename VectorType::value_type> diagonal;
347  std::vector<typename VectorType::value_type> offdiagonal;
348 
349  int it = 0;
350  double res = -std::numeric_limits<double>::max();
351 
352  typename VectorType::value_type eigen_beta_alpha = 0;
353 
354  // resize the vectors, but do not set
355  // the values since they'd be overwritten
356  // soon anyway.
357  g.reinit(x, true);
358  d.reinit(x, true);
359  h.reinit(x, true);
360 
361  number gh, beta;
362 
363  // compute residual. if vector is
364  // zero, then short-circuit the
365  // full computation
366  if (!x.all_zero())
367  {
368  A.vmult(g, x);
369  g.add(-1., b);
370  }
371  else
372  g.equ(-1., b);
373  res = g.l2_norm();
374 
375  conv = this->iteration_status(0, res, x);
376  if (conv != SolverControl::iterate)
377  return;
378 
379  if (std::is_same<PreconditionerType, PreconditionIdentity>::value == false)
380  {
381  preconditioner.vmult(h, g);
382 
383  d.equ(-1., h);
384 
385  gh = g * h;
386  }
387  else
388  {
389  d.equ(-1., g);
390  gh = res * res;
391  }
392 
393  while (conv == SolverControl::iterate)
394  {
395  it++;
396  A.vmult(h, d);
397 
398  number alpha = d * h;
399  Assert(std::abs(alpha) != 0., ExcDivideByZero());
400  alpha = gh / alpha;
401 
402  x.add(alpha, d);
403  res = std::sqrt(std::abs(g.add_and_dot(alpha, h, g)));
404 
405  print_vectors(it, x, g, d);
406 
407  conv = this->iteration_status(it, res, x);
408  if (conv != SolverControl::iterate)
409  break;
410 
411  if (std::is_same<PreconditionerType, PreconditionIdentity>::value ==
412  false)
413  {
414  preconditioner.vmult(h, g);
415 
416  beta = gh;
417  Assert(std::abs(beta) != 0., ExcDivideByZero());
418  gh = g * h;
419  beta = gh / beta;
420  d.sadd(beta, -1., h);
421  }
422  else
423  {
424  beta = gh;
425  gh = res * res;
426  beta = gh / beta;
427  d.sadd(beta, -1., g);
428  }
429 
430  this->coefficients_signal(alpha, beta);
431  // set up the vectors
432  // containing the diagonal
433  // and the off diagonal of
434  // the projected matrix.
435  if (do_eigenvalues)
436  {
437  diagonal.push_back(number(1.) / alpha + eigen_beta_alpha);
438  eigen_beta_alpha = beta / alpha;
439  offdiagonal.push_back(std::sqrt(beta) / alpha);
440  }
441  compute_eigs_and_cond(diagonal,
442  offdiagonal,
445  }
446 
447  compute_eigs_and_cond(diagonal,
448  offdiagonal,
451 
452  // in case of failure: throw exception
453  if (conv != SolverControl::success)
454  AssertThrow(false, SolverControl::NoConvergence(it, res));
455  // otherwise exit as normal
456 }
457 
458 
459 
460 template <typename VectorType>
461 boost::signals2::connection
463  const std::function<void(typename VectorType::value_type,
464  typename VectorType::value_type)> &slot)
465 {
466  return coefficients_signal.connect(slot);
467 }
468 
469 
470 
471 template <typename VectorType>
472 boost::signals2::connection
474  const std::function<void(double)> &slot,
475  const bool every_iteration)
476 {
477  if (every_iteration)
478  {
479  return all_condition_numbers_signal.connect(slot);
480  }
481  else
482  {
483  return condition_number_signal.connect(slot);
484  }
485 }
486 
487 
488 
489 template <typename VectorType>
490 boost::signals2::connection
492  const std::function<void(const std::vector<double> &)> &slot,
493  const bool every_iteration)
494 {
495  if (every_iteration)
496  {
497  return all_eigenvalues_signal.connect(slot);
498  }
499  else
500  {
501  return eigenvalues_signal.connect(slot);
502  }
503 }
504 
505 
506 
507 #endif // DOXYGEN
508 
510 
511 #endif
virtual ~SolverCG() override=default
SolverCG(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
Continue iteration.
boost::signals2::signal< void(const std::vector< double > &)> eigenvalues_signal
Definition: solver_cg.h:224
virtual void print_vectors(const unsigned int step, const VectorType &x, const VectorType &r, const VectorType &d) const
boost::signals2::signal< void(double)> all_condition_numbers_signal
Definition: solver_cg.h:218
boost::signals2::signal< SolverControl::State(const unsigned int iteration, const double check_value, const VectorType &current_iterate), StateCombiner > iteration_status
Definition: solver.h:471
#define AssertThrow(cond, exc)
Definition: exceptions.h:1521
static ::ExceptionBase & ExcDivideByZero()
Matrix is diagonal.
boost::signals2::signal< void(double)> condition_number_signal
Definition: solver_cg.h:212
boost::signals2::connection connect_condition_number_slot(const std::function< void(double)> &slot, const bool every_iteration=false)
static const char T
Stop iteration, goal reached.
#define Assert(cond, exc)
Definition: exceptions.h:1411
AdditionalData additional_data
Definition: solver_cg.h:199
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
boost::signals2::signal< void(typename VectorType::value_type, typename VectorType::value_type)> coefficients_signal
Definition: solver_cg.h:206
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
static const char A
unsigned int global_dof_index
Definition: types.h:76
boost::signals2::signal< void(const std::vector< double > &)> all_eigenvalues_signal
Definition: solver_cg.h:231
boost::signals2::connection connect_eigenvalues_slot(const std::function< void(const std::vector< double > &)> &slot, const bool every_iteration=false)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
Eigenvalue vector is filled.
static void compute_eigs_and_cond(const std::vector< typename VectorType::value_type > &diagonal, const std::vector< typename VectorType::value_type > &offdiagonal, const boost::signals2::signal< void(const std::vector< double > &)> &eigenvalues_signal, const boost::signals2::signal< void(double)> &cond_signal)
boost::signals2::connection connect_coefficients_slot(const std::function< void(typename VectorType::value_type, typename VectorType::value_type)> &slot)
T max(const T &t, const MPI_Comm &mpi_communicator)
VectorMemory< VectorType > & memory
Definition: solver.h:420