Reference documentation for deal.II version 9.2.0
\(\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 
97 template <typename VectorType = Vector<double>>
98 class SolverCG : public SolverBase<VectorType>
99 {
100 public:
105 
112  {};
113 
117  SolverCG(SolverControl & cn,
119  const AdditionalData & data = AdditionalData());
120 
125  SolverCG(SolverControl &cn, const AdditionalData &data = AdditionalData());
126 
130  virtual ~SolverCG() override = default;
131 
135  template <typename MatrixType, typename PreconditionerType>
136  void
137  solve(const MatrixType & A,
138  VectorType & x,
139  const VectorType & b,
140  const PreconditionerType &preconditioner);
141 
148  boost::signals2::connection
150  const std::function<void(typename VectorType::value_type,
151  typename VectorType::value_type)> &slot);
152 
159  boost::signals2::connection
160  connect_condition_number_slot(const std::function<void(double)> &slot,
161  const bool every_iteration = false);
162 
169  boost::signals2::connection
171  const std::function<void(const std::vector<double> &)> &slot,
172  const bool every_iteration = false);
173 
174 protected:
180  virtual void
181  print_vectors(const unsigned int step,
182  const VectorType & x,
183  const VectorType & r,
184  const VectorType & d) const;
185 
191  static void
193  const std::vector<typename VectorType::value_type> &diagonal,
194  const std::vector<typename VectorType::value_type> &offdiagonal,
195  const boost::signals2::signal<void(const std::vector<double> &)>
197  const boost::signals2::signal<void(double)> &cond_signal);
198 
203 
207  boost::signals2::signal<void(typename VectorType::value_type,
208  typename VectorType::value_type)>
210 
215  boost::signals2::signal<void(double)> condition_number_signal;
216 
221  boost::signals2::signal<void(double)> all_condition_numbers_signal;
222 
227  boost::signals2::signal<void(const std::vector<double> &)> eigenvalues_signal;
228 
233  boost::signals2::signal<void(const std::vector<double> &)>
235 };
236 
239 /*------------------------- Implementation ----------------------------*/
240 
241 #ifndef DOXYGEN
242 
243 template <typename VectorType>
246  const AdditionalData & data)
247  : SolverBase<VectorType>(cn, mem)
248  , additional_data(data)
249 {}
250 
251 
252 
253 template <typename VectorType>
256  , additional_data(data)
257 {}
258 
259 
260 
261 template <typename VectorType>
262 void
263 SolverCG<VectorType>::print_vectors(const unsigned int,
264  const VectorType &,
265  const VectorType &,
266  const VectorType &) const
267 {}
268 
269 
270 
271 template <typename VectorType>
272 inline void
274  const std::vector<typename VectorType::value_type> &diagonal,
275  const std::vector<typename VectorType::value_type> &offdiagonal,
276  const boost::signals2::signal<void(const std::vector<double> &)>
278  const boost::signals2::signal<void(double)> &cond_signal)
279 {
280  // Avoid computing eigenvalues unless they are needed.
281  if (!cond_signal.empty() || !eigenvalues_signal.empty())
282  {
284  true);
285  for (size_type i = 0; i < diagonal.size(); ++i)
286  {
287  T(i, i) = diagonal[i];
288  if (i < diagonal.size() - 1)
289  T(i, i + 1) = offdiagonal[i];
290  }
291  T.compute_eigenvalues();
292  // Need two eigenvalues to estimate the condition number.
293  if (diagonal.size() > 1)
294  {
295  auto condition_number = T.eigenvalue(T.n() - 1) / T.eigenvalue(0);
296  // Condition number is real valued and nonnegative; simply take
297  // the absolute value:
298  cond_signal(std::abs(condition_number));
299  }
300  // Avoid copying the eigenvalues of T to a vector unless a signal is
301  // connected.
302  if (!eigenvalues_signal.empty())
303  {
304  std::vector<double> eigenvalues(T.n());
305  for (unsigned int j = 0; j < T.n(); ++j)
306  {
307  // for a hermitian matrix, all eigenvalues are real-valued
308  // and non-negative, simply return the absolute value:
309  eigenvalues[j] = std::abs(T.eigenvalue(j));
310  }
312  }
313  }
314 }
315 
316 
317 
318 template <typename VectorType>
319 template <typename MatrixType, typename PreconditionerType>
320 void
321 SolverCG<VectorType>::solve(const MatrixType & A,
322  VectorType & x,
323  const VectorType & b,
324  const PreconditionerType &preconditioner)
325 {
326  using number = typename VectorType::value_type;
327 
329 
330  LogStream::Prefix prefix("cg");
331 
332  // Memory allocation
333  typename VectorMemory<VectorType>::Pointer g_pointer(this->memory);
334  typename VectorMemory<VectorType>::Pointer d_pointer(this->memory);
335  typename VectorMemory<VectorType>::Pointer h_pointer(this->memory);
336 
337  // define some aliases for simpler access
338  VectorType &g = *g_pointer;
339  VectorType &d = *d_pointer;
340  VectorType &h = *h_pointer;
341 
342  // Should we build the matrix for eigenvalue computations?
343  const bool do_eigenvalues =
345  !eigenvalues_signal.empty() || !all_eigenvalues_signal.empty();
346 
347  // vectors used for eigenvalue
348  // computations
349  std::vector<typename VectorType::value_type> diagonal;
350  std::vector<typename VectorType::value_type> offdiagonal;
351 
352  int it = 0;
353  double res = -std::numeric_limits<double>::max();
354 
355  typename VectorType::value_type eigen_beta_alpha = 0;
356 
357  // resize the vectors, but do not set
358  // the values since they'd be overwritten
359  // soon anyway.
360  g.reinit(x, true);
361  d.reinit(x, true);
362  h.reinit(x, true);
363 
364  number gh, beta;
365 
366  // compute residual. if vector is
367  // zero, then short-circuit the
368  // full computation
369  if (!x.all_zero())
370  {
371  A.vmult(g, x);
372  g.add(-1., b);
373  }
374  else
375  g.equ(-1., b);
376  res = g.l2_norm();
377 
378  conv = this->iteration_status(0, res, x);
379  if (conv != SolverControl::iterate)
380  return;
381 
383  {
384  preconditioner.vmult(h, g);
385 
386  d.equ(-1., h);
387 
388  gh = g * h;
389  }
390  else
391  {
392  d.equ(-1., g);
393  gh = res * res;
394  }
395 
396  while (conv == SolverControl::iterate)
397  {
398  it++;
399  A.vmult(h, d);
400 
401  number alpha = d * h;
402  Assert(std::abs(alpha) != 0., ExcDivideByZero());
403  alpha = gh / alpha;
404 
405  x.add(alpha, d);
406  res = std::sqrt(std::abs(g.add_and_dot(alpha, h, g)));
407 
408  print_vectors(it, x, g, d);
409 
410  conv = this->iteration_status(it, res, x);
411  if (conv != SolverControl::iterate)
412  break;
413 
415  false)
416  {
417  preconditioner.vmult(h, g);
418 
419  beta = gh;
420  Assert(std::abs(beta) != 0., ExcDivideByZero());
421  gh = g * h;
422  beta = gh / beta;
423  d.sadd(beta, -1., h);
424  }
425  else
426  {
427  beta = gh;
428  gh = res * res;
429  beta = gh / beta;
430  d.sadd(beta, -1., g);
431  }
432 
433  this->coefficients_signal(alpha, beta);
434  // set up the vectors
435  // containing the diagonal
436  // and the off diagonal of
437  // the projected matrix.
438  if (do_eigenvalues)
439  {
440  diagonal.push_back(number(1.) / alpha + eigen_beta_alpha);
441  eigen_beta_alpha = beta / alpha;
442  offdiagonal.push_back(std::sqrt(beta) / alpha);
443  }
444  compute_eigs_and_cond(diagonal,
445  offdiagonal,
448  }
449 
450  compute_eigs_and_cond(diagonal,
451  offdiagonal,
454 
455  // in case of failure: throw exception
456  if (conv != SolverControl::success)
457  AssertThrow(false, SolverControl::NoConvergence(it, res));
458  // otherwise exit as normal
459 }
460 
461 
462 
463 template <typename VectorType>
464 boost::signals2::connection
466  const std::function<void(typename VectorType::value_type,
467  typename VectorType::value_type)> &slot)
468 {
469  return coefficients_signal.connect(slot);
470 }
471 
472 
473 
474 template <typename VectorType>
475 boost::signals2::connection
477  const std::function<void(double)> &slot,
478  const bool every_iteration)
479 {
480  if (every_iteration)
481  {
482  return all_condition_numbers_signal.connect(slot);
483  }
484  else
485  {
486  return condition_number_signal.connect(slot);
487  }
488 }
489 
490 
491 
492 template <typename VectorType>
493 boost::signals2::connection
495  const std::function<void(const std::vector<double> &)> &slot,
496  const bool every_iteration)
497 {
498  if (every_iteration)
499  {
500  return all_eigenvalues_signal.connect(slot);
501  }
502  else
503  {
504  return eigenvalues_signal.connect(slot);
505  }
506 }
507 
508 
509 
510 #endif // DOXYGEN
511 
513 
514 #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:227
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:221
boost::signals2::signal< SolverControl::State(const unsigned int iteration, const double check_value, const VectorType &current_iterate), StateCombiner > iteration_status
Definition: solver.h:463
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
static ::ExceptionBase & ExcDivideByZero()
Matrix is diagonal.
boost::signals2::signal< void(double)> condition_number_signal
Definition: solver_cg.h:215
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:1419
AdditionalData additional_data
Definition: solver_cg.h:202
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
boost::signals2::signal< void(typename VectorType::value_type, typename VectorType::value_type)> coefficients_signal
Definition: solver_cg.h:209
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:234
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:358
Eigenvalue vector is filled.
static const bool value
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:412