Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
solver_cg.h
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 
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/logstream.h>
24 #include <deal.II/base/subscriptor.h>
25 
26 #include <deal.II/lac/solver.h>
27 #include <deal.II/lac/solver_control.h>
28 #include <deal.II/lac/tridiagonal_matrix.h>
29 
30 #include <cmath>
31 
32 DEAL_II_NAMESPACE_OPEN
33 
34 // forward declaration
36 
37 
40 
95 template <typename VectorType = Vector<double>>
96 class SolverCG : public SolverBase<VectorType>
97 {
98 public:
103 
110  {};
111 
115  SolverCG(SolverControl & cn,
117  const AdditionalData & data = AdditionalData());
118 
123  SolverCG(SolverControl &cn, const AdditionalData &data = AdditionalData());
124 
128  virtual ~SolverCG() override = default;
129 
133  template <typename MatrixType, typename PreconditionerType>
134  void
135  solve(const MatrixType & A,
136  VectorType & x,
137  const VectorType & b,
138  const PreconditionerType &preconditioner);
139 
146  boost::signals2::connection
148  const std::function<void(typename VectorType::value_type,
149  typename VectorType::value_type)> &slot);
150 
157  boost::signals2::connection
158  connect_condition_number_slot(const std::function<void(double)> &slot,
159  const bool every_iteration = false);
160 
167  boost::signals2::connection
169  const std::function<void(const std::vector<double> &)> &slot,
170  const bool every_iteration = false);
171 
172 protected:
178  virtual void
179  print_vectors(const unsigned int step,
180  const VectorType & x,
181  const VectorType & r,
182  const VectorType & d) const;
183 
189  static void
191  const std::vector<typename VectorType::value_type> &diagonal,
192  const std::vector<typename VectorType::value_type> &offdiagonal,
193  const boost::signals2::signal<void(const std::vector<double> &)>
195  const boost::signals2::signal<void(double)> &cond_signal);
196 
201 
205  boost::signals2::signal<void(typename VectorType::value_type,
206  typename VectorType::value_type)>
208 
213  boost::signals2::signal<void(double)> condition_number_signal;
214 
219  boost::signals2::signal<void(double)> all_condition_numbers_signal;
220 
225  boost::signals2::signal<void(const std::vector<double> &)> eigenvalues_signal;
226 
231  boost::signals2::signal<void(const std::vector<double> &)>
233 };
234 
237 /*------------------------- Implementation ----------------------------*/
238 
239 #ifndef DOXYGEN
240 
241 template <typename VectorType>
244  const AdditionalData & data)
245  : SolverBase<VectorType>(cn, mem)
246  , additional_data(data)
247 {}
248 
249 
250 
251 template <typename VectorType>
252 SolverCG<VectorType>::SolverCG(SolverControl &cn, const AdditionalData &data)
253  : SolverBase<VectorType>(cn)
254  , additional_data(data)
255 {}
256 
257 
258 
259 template <typename VectorType>
260 void
261 SolverCG<VectorType>::print_vectors(const unsigned int,
262  const VectorType &,
263  const VectorType &,
264  const VectorType &) const
265 {}
266 
267 
268 
269 template <typename VectorType>
270 inline void
272  const std::vector<typename VectorType::value_type> &diagonal,
273  const std::vector<typename VectorType::value_type> &offdiagonal,
274  const boost::signals2::signal<void(const std::vector<double> &)>
275  & eigenvalues_signal,
276  const boost::signals2::signal<void(double)> &cond_signal)
277 {
278  // Avoid computing eigenvalues unless they are needed.
279  if (!cond_signal.empty() || !eigenvalues_signal.empty())
280  {
282  true);
283  for (size_type i = 0; i < diagonal.size(); ++i)
284  {
285  T(i, i) = diagonal[i];
286  if (i < diagonal.size() - 1)
287  T(i, i + 1) = offdiagonal[i];
288  }
289  T.compute_eigenvalues();
290  // Need two eigenvalues to estimate the condition number.
291  if (diagonal.size() > 1)
292  {
293  auto condition_number = T.eigenvalue(T.n() - 1) / T.eigenvalue(0);
294  // Condition number is real valued and nonnegative; simply take
295  // the absolute value:
296  cond_signal(std::abs(condition_number));
297  }
298  // Avoid copying the eigenvalues of T to a vector unless a signal is
299  // connected.
300  if (!eigenvalues_signal.empty())
301  {
302  std::vector<double> eigenvalues(T.n());
303  for (unsigned int j = 0; j < T.n(); ++j)
304  {
305  // for a hermitian matrix, all eigenvalues are real-valued
306  // and non-negative, simply return the absolute value:
307  eigenvalues[j] = std::abs(T.eigenvalue(j));
308  }
309  eigenvalues_signal(eigenvalues);
310  }
311  }
312 }
313 
314 
315 
316 template <typename VectorType>
317 template <typename MatrixType, typename PreconditionerType>
318 void
319 SolverCG<VectorType>::solve(const MatrixType & A,
320  VectorType & x,
321  const VectorType & b,
322  const PreconditionerType &preconditioner)
323 {
324  using number = typename VectorType::value_type;
325 
327 
328  LogStream::Prefix prefix("cg");
329 
330  // Memory allocation
331  typename VectorMemory<VectorType>::Pointer g_pointer(this->memory);
332  typename VectorMemory<VectorType>::Pointer d_pointer(this->memory);
333  typename VectorMemory<VectorType>::Pointer h_pointer(this->memory);
334 
335  // define some aliases for simpler access
336  VectorType &g = *g_pointer;
337  VectorType &d = *d_pointer;
338  VectorType &h = *h_pointer;
339 
340  // Should we build the matrix for eigenvalue computations?
341  const bool do_eigenvalues =
342  !condition_number_signal.empty() || !all_condition_numbers_signal.empty() ||
343  !eigenvalues_signal.empty() || !all_eigenvalues_signal.empty();
344 
345  // vectors used for eigenvalue
346  // computations
347  std::vector<typename VectorType::value_type> diagonal;
348  std::vector<typename VectorType::value_type> offdiagonal;
349 
350  int it = 0;
351  double res = -std::numeric_limits<double>::max();
352 
353  typename VectorType::value_type eigen_beta_alpha = 0;
354 
355  // resize the vectors, but do not set
356  // the values since they'd be overwritten
357  // soon anyway.
358  g.reinit(x, true);
359  d.reinit(x, true);
360  h.reinit(x, true);
361 
362  number gh, beta;
363 
364  // compute residual. if vector is
365  // zero, then short-circuit the
366  // full computation
367  if (!x.all_zero())
368  {
369  A.vmult(g, x);
370  g.add(-1., b);
371  }
372  else
373  g.equ(-1., b);
374  res = g.l2_norm();
375 
376  conv = this->iteration_status(0, res, x);
377  if (conv != SolverControl::iterate)
378  return;
379 
380  if (std::is_same<PreconditionerType, PreconditionIdentity>::value == false)
381  {
382  preconditioner.vmult(h, g);
383 
384  d.equ(-1., h);
385 
386  gh = g * h;
387  }
388  else
389  {
390  d.equ(-1., g);
391  gh = res * res;
392  }
393 
394  while (conv == SolverControl::iterate)
395  {
396  it++;
397  A.vmult(h, d);
398 
399  number alpha = d * h;
400  Assert(std::abs(alpha) != 0., ExcDivideByZero());
401  alpha = gh / alpha;
402 
403  x.add(alpha, d);
404  res = std::sqrt(std::abs(g.add_and_dot(alpha, h, g)));
405 
406  print_vectors(it, x, g, d);
407 
408  conv = this->iteration_status(it, res, x);
409  if (conv != SolverControl::iterate)
410  break;
411 
412  if (std::is_same<PreconditionerType, PreconditionIdentity>::value ==
413  false)
414  {
415  preconditioner.vmult(h, g);
416 
417  beta = gh;
418  Assert(std::abs(beta) != 0., ExcDivideByZero());
419  gh = g * h;
420  beta = gh / beta;
421  d.sadd(beta, -1., h);
422  }
423  else
424  {
425  beta = gh;
426  gh = res * res;
427  beta = gh / beta;
428  d.sadd(beta, -1., g);
429  }
430 
431  this->coefficients_signal(alpha, beta);
432  // set up the vectors
433  // containing the diagonal
434  // and the off diagonal of
435  // the projected matrix.
436  if (do_eigenvalues)
437  {
438  diagonal.push_back(number(1.) / alpha + eigen_beta_alpha);
439  eigen_beta_alpha = beta / alpha;
440  offdiagonal.push_back(std::sqrt(beta) / alpha);
441  }
442  compute_eigs_and_cond(diagonal,
443  offdiagonal,
444  all_eigenvalues_signal,
445  all_condition_numbers_signal);
446  }
447 
448  compute_eigs_and_cond(diagonal,
449  offdiagonal,
450  eigenvalues_signal,
451  condition_number_signal);
452 
453  // in case of failure: throw exception
454  if (conv != SolverControl::success)
455  AssertThrow(false, SolverControl::NoConvergence(it, res));
456  // otherwise exit as normal
457 }
458 
459 
460 
461 template <typename VectorType>
462 boost::signals2::connection
464  const std::function<void(typename VectorType::value_type,
465  typename VectorType::value_type)> &slot)
466 {
467  return coefficients_signal.connect(slot);
468 }
469 
470 
471 
472 template <typename VectorType>
473 boost::signals2::connection
475  const std::function<void(double)> &slot,
476  const bool every_iteration)
477 {
478  if (every_iteration)
479  {
480  return all_condition_numbers_signal.connect(slot);
481  }
482  else
483  {
484  return condition_number_signal.connect(slot);
485  }
486 }
487 
488 
489 
490 template <typename VectorType>
491 boost::signals2::connection
493  const std::function<void(const std::vector<double> &)> &slot,
494  const bool every_iteration)
495 {
496  if (every_iteration)
497  {
498  return all_eigenvalues_signal.connect(slot);
499  }
500  else
501  {
502  return eigenvalues_signal.connect(slot);
503  }
504 }
505 
506 
507 
508 #endif // DOXYGEN
509 
510 DEAL_II_NAMESPACE_CLOSE
511 
512 #endif
virtual ~SolverCG() override=default
types::global_dof_index size_type
Definition: solver_cg.h:102
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:225
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:219
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
static ::ExceptionBase & ExcDivideByZero()
Matrix is diagonal.
boost::signals2::signal< void(double)> condition_number_signal
Definition: solver_cg.h:213
boost::signals2::connection connect_condition_number_slot(const std::function< void(double)> &slot, const bool every_iteration=false)
Stop iteration, goal reached.
#define Assert(cond, exc)
Definition: exceptions.h:1407
AdditionalData additional_data
Definition: solver_cg.h:200
boost::signals2::signal< void(typename VectorType::value_type, typename VectorType::value_type)> coefficients_signal
Definition: solver_cg.h:207
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
unsigned int global_dof_index
Definition: types.h:89
boost::signals2::signal< void(const std::vector< double > &)> all_eigenvalues_signal
Definition: solver_cg.h:232
boost::signals2::connection connect_eigenvalues_slot(const std::function< void(const std::vector< double > &)> &slot, const bool every_iteration=false)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
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)