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_idr.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 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_idr_h
17 #define dealii_solver_idr_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/logstream.h>
25 #include <deal.II/base/utilities.h>
26 
28 #include <deal.II/lac/solver.h>
30 
31 #include <cmath>
32 #include <random>
33 
35 
38 
39 namespace internal
40 {
44  namespace SolverIDRImplementation
45  {
50  template <typename VectorType>
51  class TmpVectors
52  {
53  public:
57  TmpVectors(const unsigned int s_param, VectorMemory<VectorType> &vmem);
58 
62  ~TmpVectors() = default;
63 
68  VectorType &operator[](const unsigned int i) const;
69 
76  VectorType &
77  operator()(const unsigned int i, const VectorType &temp);
78 
79  private:
84 
88  std::vector<typename VectorMemory<VectorType>::Pointer> data;
89  };
90  } // namespace SolverIDRImplementation
91 } // namespace internal
92 
117 template <class VectorType = Vector<double>>
118 class SolverIDR : public SolverBase<VectorType>
119 {
120 public:
125  {
129  explicit AdditionalData(const unsigned int s = 2)
130  : s(s)
131  {}
132 
133  const unsigned int s;
134  };
135 
141  const AdditionalData & data = AdditionalData());
142 
147  explicit SolverIDR(SolverControl & cn,
148  const AdditionalData &data = AdditionalData());
149 
153  virtual ~SolverIDR() override = default;
154 
158  template <typename MatrixType, typename PreconditionerType>
159  void
160  solve(const MatrixType & A,
161  VectorType & x,
162  const VectorType & b,
163  const PreconditionerType &preconditioner);
164 
165 protected:
171  virtual void
172  print_vectors(const unsigned int step,
173  const VectorType & x,
174  const VectorType & r,
175  const VectorType & d) const;
176 
177 private:
182 };
183 
185 /*------------------------- Implementation ----------------------------*/
186 
187 #ifndef DOXYGEN
188 
189 
190 namespace internal
191 {
192  namespace SolverIDRImplementation
193  {
194  template <class VectorType>
195  inline TmpVectors<VectorType>::TmpVectors(const unsigned int s_param,
197  : mem(vmem)
198  , data(s_param)
199  {}
200 
201 
202 
203  template <class VectorType>
205  operator[](const unsigned int i) const
206  {
207  AssertIndexRange(i, data.size());
208 
209  Assert(data[i] != nullptr, ExcNotInitialized());
210  return *data[i];
211  }
212 
213 
214 
215  template <class VectorType>
216  inline VectorType &
217  TmpVectors<VectorType>::operator()(const unsigned int i,
218  const VectorType & temp)
219  {
220  AssertIndexRange(i, data.size());
221  if (data[i] == nullptr)
222  {
223  data[i] = std::move(typename VectorMemory<VectorType>::Pointer(mem));
224  data[i]->reinit(temp);
225  }
226  return *data[i];
227  }
228  } // namespace SolverIDRImplementation
229 } // namespace internal
230 
231 
232 
233 template <class VectorType>
236  const AdditionalData & data)
238  , additional_data(data)
239 {}
240 
241 
242 
243 template <class VectorType>
246  , additional_data(data)
247 {}
248 
249 
250 
251 template <class VectorType>
252 void
253 SolverIDR<VectorType>::print_vectors(const unsigned int,
254  const VectorType &,
255  const VectorType &,
256  const VectorType &) const
257 {}
258 
259 
260 
261 template <class VectorType>
262 template <typename MatrixType, typename PreconditionerType>
263 void
264 SolverIDR<VectorType>::solve(const MatrixType & A,
265  VectorType & x,
266  const VectorType & b,
267  const PreconditionerType &preconditioner)
268 {
269  LogStream::Prefix prefix("IDR(s)");
270 
272  unsigned int step = 0;
273 
274  const unsigned int s = additional_data.s;
275 
276  // Define temporary vectors which do not do not depend on s
277  typename VectorMemory<VectorType>::Pointer r_pointer(this->memory);
278  typename VectorMemory<VectorType>::Pointer v_pointer(this->memory);
279  typename VectorMemory<VectorType>::Pointer vhat_pointer(this->memory);
280  typename VectorMemory<VectorType>::Pointer uhat_pointer(this->memory);
281  typename VectorMemory<VectorType>::Pointer ghat_pointer(this->memory);
282 
283  VectorType &r = *r_pointer;
284  VectorType &v = *v_pointer;
285  VectorType &vhat = *vhat_pointer;
286  VectorType &uhat = *uhat_pointer;
287  VectorType &ghat = *ghat_pointer;
288 
289  r.reinit(x, true);
290  v.reinit(x, true);
291  vhat.reinit(x, true);
292  uhat.reinit(x, true);
293  ghat.reinit(x, true);
294 
295  // Initial residual
296  A.vmult(r, x);
297  r.sadd(-1.0, 1.0, b);
298 
299  // Check for convergent initial guess
300  double res = r.l2_norm();
301  iteration_state = this->iteration_status(step, res, x);
302  if (iteration_state == SolverControl::success)
303  return;
304 
305  // Initialize sets of vectors/matrices whose size dependent on s
309  FullMatrix<double> M(s, s);
310 
311  // Random number generator for vector entries of
312  // Q (normal distribution, mean=0 sigma=1)
313  std::mt19937 rng;
314  std::normal_distribution<> normal_distribution(0.0, 1.0);
315  for (unsigned int i = 0; i < s; ++i)
316  {
317  VectorType &tmp_g = G(i, x);
318  VectorType &tmp_u = U(i, x);
319  tmp_g = 0;
320  tmp_u = 0;
321 
322  // Compute random set of s orthonormalized vectors Q
323  // Note: the first vector is chosen to be the initial
324  // residual to match BiCGStab (as is done in comparisons
325  // with BiCGStab in the papers listed in the documentation
326  // of this function)
327  VectorType &tmp_q = Q(i, x);
328  if (i != 0)
329  {
330  for (auto indx : tmp_q.locally_owned_elements())
331  tmp_q(indx) = normal_distribution(rng);
332  tmp_q.compress(VectorOperation::insert);
333  }
334  else
335  tmp_q = r;
336 
337  for (unsigned int j = 0; j < i; ++j)
338  {
339  v = Q[j];
340  v *= (v * tmp_q) / (tmp_q * tmp_q);
341  tmp_q.add(-1.0, v);
342  }
343 
344  if (i != 0)
345  tmp_q *= 1.0 / tmp_q.l2_norm();
346 
347  M(i, i) = 1.;
348  }
349 
350  double omega = 1.;
351 
352  bool early_exit = false;
353 
354  // Outer iteration
355  while (iteration_state == SolverControl::iterate)
356  {
357  ++step;
358 
359  // Compute phi
360  Vector<double> phi(s);
361  for (unsigned int i = 0; i < s; ++i)
362  phi(i) = Q[i] * r;
363 
364  // Inner iteration over s
365  for (unsigned int k = 0; k < s; ++k)
366  {
367  // Solve M(k:s)*gamma = phi(k:s)
368  Vector<double> gamma(s - k);
369  {
370  Vector<double> phik(s - k);
371  FullMatrix<double> Mk(s - k, s - k);
372  std::vector<unsigned int> indices;
373  unsigned int j = 0;
374  for (unsigned int i = k; i < s; ++i, ++j)
375  {
376  indices.push_back(i);
377  phik(j) = phi(i);
378  }
379  Mk.extract_submatrix_from(M, indices, indices);
380 
381  FullMatrix<double> Mk_inv(s - k, s - k);
382  Mk_inv.invert(Mk);
383  Mk_inv.vmult(gamma, phik);
384  }
385 
386  {
387  v = r;
388 
389  unsigned int j = 0;
390  for (unsigned int i = k; i < s; ++i, ++j)
391  v.add(-1.0 * gamma(j), G[i]);
392  preconditioner.vmult(vhat, v);
393 
394  uhat = vhat;
395  uhat *= omega;
396  j = 0;
397  for (unsigned int i = k; i < s; ++i, ++j)
398  uhat.add(gamma(j), U[i]);
399  A.vmult(ghat, uhat);
400  }
401 
402  // Update G and U
403  // Orthogonalize ghat to Q0,..,Q_{k-1}
404  // and update uhat
405  for (unsigned int i = 0; i < k; ++i)
406  {
407  double alpha = (Q[i] * ghat) / M(i, i);
408  ghat.add(-alpha, G[i]);
409  uhat.add(-alpha, U[i]);
410  }
411  G[k] = ghat;
412  U[k] = uhat;
413 
414  // Update kth column of M
415  for (unsigned int i = k; i < s; ++i)
416  M(i, k) = Q[i] * G[k];
417 
418  // Orthogonalize r to Q0,...,Qk,
419  // update x
420  {
421  double beta = phi(k) / M(k, k);
422  r.add(-1.0 * beta, G[k]);
423  x.add(beta, U[k]);
424 
425  print_vectors(step, x, r, U[k]);
426 
427  // Check for early convergence. If so, store
428  // information in early_exit so that outer iteration
429  // is broken before recomputing the residual
430  res = r.l2_norm();
431  iteration_state = this->iteration_status(step, res, x);
432  if (iteration_state != SolverControl::iterate)
433  {
434  early_exit = true;
435  break;
436  }
437 
438  // Update phi
439  if (k + 1 < s)
440  {
441  for (unsigned int i = 0; i < k + 1; ++i)
442  phi(i) = 0.0;
443  for (unsigned int i = k + 1; i < s; ++i)
444  phi(i) -= beta * M(i, k);
445  }
446  }
447  }
448  if (early_exit == true)
449  break;
450 
451  // Update r and x
452  preconditioner.vmult(vhat, r);
453  A.vmult(v, vhat);
454 
455  omega = (v * r) / (v * v);
456 
457  r.add(-1.0 * omega, v);
458  x.add(omega, vhat);
459 
460  print_vectors(step, x, r, vhat);
461 
462  // Check for convergence
463  res = r.l2_norm();
464  iteration_state = this->iteration_status(step, res, x);
465  if (iteration_state != SolverControl::iterate)
466  break;
467  }
468 
469  if (iteration_state != SolverControl::success)
470  AssertThrow(false, SolverControl::NoConvergence(step, res));
471 }
472 
473 
474 #endif // DOXYGEN
475 
477 
478 #endif
const unsigned int s
Definition: solver_idr.h:133
TmpVectors(const unsigned int s_param, VectorMemory< VectorType > &vmem)
Continue iteration.
AdditionalData(const unsigned int s=2)
Definition: solver_idr.h:129
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
std::vector< typename VectorMemory< VectorType >::Pointer > data
Definition: solver_idr.h:88
void invert(const FullMatrix< number2 > &M)
static const char U
static ::ExceptionBase & ExcNotInitialized()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
VectorType & operator()(const unsigned int i, const VectorType &temp)
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
AdditionalData additional_data
Definition: solver_idr.h:181
Stop iteration, goal reached.
SolverIDR(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
#define Assert(cond, exc)
Definition: exceptions.h:1419
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
void extract_submatrix_from(const MatrixType &matrix, const std::vector< index_type > &row_index_set, const std::vector< index_type > &column_index_set)
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)
virtual void print_vectors(const unsigned int step, const VectorType &x, const VectorType &r, const VectorType &d) const
static const char A
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
VectorMemory< VectorType > & mem
Definition: solver_idr.h:83
VectorType & operator[](const unsigned int i) const
long double gamma(const unsigned int n)