Reference documentation for deal.II version 9.3.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 - 2020 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 computations
345  std::vector<typename VectorType::value_type> diagonal;
346  std::vector<typename VectorType::value_type> offdiagonal;
347 
348  typename VectorType::value_type eigen_beta_alpha = 0;
349 
350  // resize the vectors, but do not set the values since they'd be overwritten
351  // soon anyway.
352  g.reinit(x, true);
353  d.reinit(x, true);
354  h.reinit(x, true);
355 
356  int it = 0;
357  number gh = number();
358  number beta = number();
359  number alpha = number();
360  number old_alpha = number();
361 
362  // compute residual. if vector is zero, then short-circuit the full
363  // computation
364  if (!x.all_zero())
365  {
366  A.vmult(g, x);
367  g.add(-1., b);
368  }
369  else
370  g.equ(-1., b);
371 
372  double res = g.l2_norm();
373  conv = this->iteration_status(0, res, x);
374  if (conv != SolverControl::iterate)
375  return;
376 
377  while (conv == SolverControl::iterate)
378  {
379  it++;
380  old_alpha = alpha;
381 
382  if (it > 1)
383  {
384  if (std::is_same<PreconditionerType, PreconditionIdentity>::value ==
385  false)
386  {
387  preconditioner.vmult(h, g);
388  beta = gh;
389  Assert(std::abs(beta) != 0., ExcDivideByZero());
390  gh = g * h;
391  beta = gh / beta;
392  d.sadd(beta, -1., h);
393  }
394  else
395  {
396  beta = gh;
397  gh = res * res;
398  beta = gh / beta;
399  d.sadd(beta, -1., g);
400  }
401  }
402  else
403  {
404  if (std::is_same<PreconditionerType, PreconditionIdentity>::value ==
405  false)
406  {
407  preconditioner.vmult(h, g);
408  d.equ(-1., h);
409  gh = g * h;
410  }
411  else
412  {
413  d.equ(-1., g);
414  gh = res * res;
415  }
416  }
417 
418  A.vmult(h, d);
419 
420  alpha = d * h;
421  Assert(std::abs(alpha) != 0., ExcDivideByZero());
422  alpha = gh / alpha;
423 
424  x.add(alpha, d);
425  res = std::sqrt(std::abs(g.add_and_dot(alpha, h, g)));
426 
427  print_vectors(it, x, g, d);
428 
429  if (it > 1)
430  {
431  this->coefficients_signal(old_alpha, beta);
432  // set up the vectors containing the diagonal and the off diagonal of
433  // the projected matrix.
434  if (do_eigenvalues)
435  {
436  diagonal.push_back(number(1.) / old_alpha + eigen_beta_alpha);
437  eigen_beta_alpha = beta / old_alpha;
438  offdiagonal.push_back(std::sqrt(beta) / old_alpha);
439  }
440  compute_eigs_and_cond(diagonal,
441  offdiagonal,
444  }
445 
446  conv = this->iteration_status(it, res, x);
447  }
448 
449  compute_eigs_and_cond(diagonal,
450  offdiagonal,
453 
454  // in case of failure: throw exception
455  if (conv != SolverControl::success)
456  AssertThrow(false, SolverControl::NoConvergence(it, res));
457  // otherwise exit as normal
458 }
459 
460 
461 
462 template <typename VectorType>
463 boost::signals2::connection
465  const std::function<void(typename VectorType::value_type,
466  typename VectorType::value_type)> &slot)
467 {
468  return coefficients_signal.connect(slot);
469 }
470 
471 
472 
473 template <typename VectorType>
474 boost::signals2::connection
476  const std::function<void(double)> &slot,
477  const bool every_iteration)
478 {
479  if (every_iteration)
480  {
481  return all_condition_numbers_signal.connect(slot);
482  }
483  else
484  {
485  return condition_number_signal.connect(slot);
486  }
487 }
488 
489 
490 
491 template <typename VectorType>
492 boost::signals2::connection
494  const std::function<void(const std::vector<double> &)> &slot,
495  const bool every_iteration)
496 {
497  if (every_iteration)
498  {
499  return all_eigenvalues_signal.connect(slot);
500  }
501  else
502  {
503  return eigenvalues_signal.connect(slot);
504  }
505 }
506 
507 
508 
509 #endif // DOXYGEN
510 
512 
513 #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:1575
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:1465
AdditionalData additional_data
Definition: solver_cg.h:199
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
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:394
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)
VectorMemory< VectorType > & memory
Definition: solver.h:420