Reference documentation for deal.II version GIT 35969cdc9b 2023-12-09 01:10:02+00:00
\(\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 - 2022 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 
28 #include <deal.II/lac/solver.h>
30 
31 #include <cmath>
32 #include <random>
33 
35 
41 namespace internal
42 {
46  namespace SolverIDRImplementation
47  {
52  template <typename VectorType>
53  class TmpVectors
54  {
55  public:
59  TmpVectors(const unsigned int s_param, VectorMemory<VectorType> &vmem);
60 
64  ~TmpVectors() = default;
65 
70  VectorType &
71  operator[](const unsigned int i) const;
72 
79  VectorType &
80  operator()(const unsigned int i, const VectorType &temp);
81 
82  private:
87 
91  std::vector<typename VectorMemory<VectorType>::Pointer> data;
92  };
93  } // namespace SolverIDRImplementation
94 } // namespace internal
95 
118 template <typename VectorType = Vector<double>>
119 class SolverIDR : public SolverBase<VectorType>
120 {
121 public:
126  {
130  explicit AdditionalData(const unsigned int s = 2)
131  : s(s)
132  {}
133 
134  const unsigned int s;
135  };
136 
142  const AdditionalData &data = AdditionalData());
143 
148  explicit SolverIDR(SolverControl &cn,
149  const AdditionalData &data = AdditionalData());
150 
154  virtual ~SolverIDR() override = default;
155 
159  template <typename MatrixType, typename PreconditionerType>
160  void
161  solve(const MatrixType &A,
162  VectorType &x,
163  const VectorType &b,
164  const PreconditionerType &preconditioner);
165 
166 protected:
172  virtual void
173  print_vectors(const unsigned int step,
174  const VectorType &x,
175  const VectorType &r,
176  const VectorType &d) const;
177 
178 private:
183 };
184 
186 /*------------------------- Implementation ----------------------------*/
187 
188 #ifndef DOXYGEN
189 
190 
191 namespace internal
192 {
193  namespace SolverIDRImplementation
194  {
195  template <typename VectorType>
196  inline TmpVectors<VectorType>::TmpVectors(const unsigned int s_param,
198  : mem(vmem)
199  , data(s_param)
200  {}
201 
202 
203 
204  template <typename VectorType>
205  inline VectorType &
206  TmpVectors<VectorType>::operator[](const unsigned int i) const
207  {
208  AssertIndexRange(i, data.size());
209 
210  Assert(data[i] != nullptr, ExcNotInitialized());
211  return *data[i];
212  }
213 
214 
215 
216  template <typename VectorType>
217  inline VectorType &
218  TmpVectors<VectorType>::operator()(const unsigned int i,
219  const VectorType &temp)
220  {
221  AssertIndexRange(i, data.size());
222  if (data[i] == nullptr)
223  {
224  data[i] = std::move(typename VectorMemory<VectorType>::Pointer(mem));
225  data[i]->reinit(temp, true);
226  }
227  return *data[i];
228  }
229 
230 
231 
232  template <typename VectorType,
233  std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
234  * = nullptr>
235  unsigned int
236  n_blocks(const VectorType &)
237  {
238  return 1;
239  }
240 
241 
242 
243  template <typename VectorType,
244  std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
245  nullptr>
246  unsigned int
247  n_blocks(const VectorType &vector)
248  {
249  return vector.n_blocks();
250  }
251 
252 
253 
254  template <typename VectorType,
255  std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
256  * = nullptr>
257  VectorType &
258  block(VectorType &vector, const unsigned int b)
259  {
260  AssertDimension(b, 0);
261  (void)b;
262  return vector;
263  }
264 
265 
266 
267  template <typename VectorType,
268  std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
269  nullptr>
270  typename VectorType::BlockType &
271  block(VectorType &vector, const unsigned int b)
272  {
273  return vector.block(b);
274  }
275 
276  } // namespace SolverIDRImplementation
277 } // namespace internal
278 
279 
280 
281 template <typename VectorType>
284  const AdditionalData &data)
285  : SolverBase<VectorType>(cn, mem)
286  , additional_data(data)
287 {}
288 
289 
290 
291 template <typename VectorType>
292 SolverIDR<VectorType>::SolverIDR(SolverControl &cn, const AdditionalData &data)
293  : SolverBase<VectorType>(cn)
294  , additional_data(data)
295 {}
296 
297 
298 
299 template <typename VectorType>
300 void
301 SolverIDR<VectorType>::print_vectors(const unsigned int,
302  const VectorType &,
303  const VectorType &,
304  const VectorType &) const
305 {}
306 
307 
308 
309 template <typename VectorType>
310 template <typename MatrixType, typename PreconditionerType>
311 void
312 SolverIDR<VectorType>::solve(const MatrixType &A,
313  VectorType &x,
314  const VectorType &b,
315  const PreconditionerType &preconditioner)
316 {
317  LogStream::Prefix prefix("IDR(s)");
318 
320  unsigned int step = 0;
321 
322  const unsigned int s = additional_data.s;
323 
324  // Define temporary vectors which do not do not depend on s
325  typename VectorMemory<VectorType>::Pointer r_pointer(this->memory);
326  typename VectorMemory<VectorType>::Pointer v_pointer(this->memory);
327  typename VectorMemory<VectorType>::Pointer uhat_pointer(this->memory);
328 
329  VectorType &r = *r_pointer;
330  VectorType &v = *v_pointer;
331  VectorType &uhat = *uhat_pointer;
332 
333  r.reinit(x, true);
334  v.reinit(x, true);
335  uhat.reinit(x, true);
336 
337  // Initial residual
338  A.vmult(r, x);
339  r.sadd(-1.0, 1.0, b);
340 
341  using value_type = typename VectorType::value_type;
342  using real_type = typename numbers::NumberTraits<value_type>::real_type;
343 
344  // Check for convergent initial guess
345  real_type res = r.l2_norm();
346  iteration_state = this->iteration_status(step, res, x);
347  if (iteration_state == SolverControl::success)
348  return;
349 
350  // Initialize sets of vectors/matrices whose size dependent on s
354  FullMatrix<value_type> M(s, s);
355 
356  // Random number generator for vector entries of
357  // Q (normal distribution, mean=0 sigma=1)
358  std::mt19937 rng;
359  std::normal_distribution<> normal_distribution(0.0, 1.0);
360  for (unsigned int i = 0; i < s; ++i)
361  {
362  // Initialize vectors
363  G(i, x);
364  U(i, x);
365 
366  // Compute random set of s orthonormalized vectors Q
367  // Note: the first vector is chosen to be the initial
368  // residual to match BiCGStab (as is done in comparisons
369  // with BiCGStab in the papers listed in the documentation
370  // of this function)
371  VectorType &tmp_q = Q(i, x);
372  if (i != 0)
373  {
374  for (unsigned int b = 0;
376  ++b)
377  for (auto indx : internal::SolverIDRImplementation::block(tmp_q, b)
378  .locally_owned_elements())
379  internal::SolverIDRImplementation::block(tmp_q, b)(indx) =
380  normal_distribution(rng);
381  tmp_q.compress(VectorOperation::insert);
382  }
383  else
384  tmp_q = r;
385 
386  for (unsigned int j = 0; j < i; ++j)
387  {
388  v = Q[j];
389  v *= (v * tmp_q) / (tmp_q * tmp_q);
390  tmp_q.add(-1.0, v);
391  }
392 
393  if (i != 0)
394  tmp_q *= 1.0 / tmp_q.l2_norm();
395 
396  M(i, i) = 1.;
397  }
398 
399  value_type omega = 1.;
400 
401  bool early_exit = false;
402 
403  // Outer iteration
404  while (iteration_state == SolverControl::iterate)
405  {
406  ++step;
407 
408  // Compute phi
409  Vector<value_type> phi(s);
410  for (unsigned int i = 0; i < s; ++i)
411  phi(i) = Q[i] * r;
412 
413  // Inner iteration over s
414  for (unsigned int k = 0; k < s; ++k)
415  {
416  // Solve M(k:s)*gamma = phi(k:s)
417  Vector<value_type> gamma(s - k);
418  {
419  Vector<value_type> phik(s - k);
420  FullMatrix<value_type> Mk(s - k, s - k);
421  std::vector<unsigned int> indices;
422  unsigned int j = 0;
423  for (unsigned int i = k; i < s; ++i, ++j)
424  {
425  indices.push_back(i);
426  phik(j) = phi(i);
427  }
428  Mk.extract_submatrix_from(M, indices, indices);
429 
430  FullMatrix<value_type> Mk_inv(s - k, s - k);
431  Mk_inv.invert(Mk);
432  Mk_inv.vmult(gamma, phik);
433  }
434 
435  v = r;
436 
437  if (step > 1)
438  {
439  for (unsigned int i = k, j = 0; i < s; ++i, ++j)
440  v.add(-gamma(j), G[i]);
441  }
442 
443  preconditioner.vmult(uhat, v);
444 
445  if (step > 1)
446  {
447  uhat.sadd(omega, gamma(0), U[k]);
448  for (unsigned int i = k + 1, j = 1; i < s; ++i, ++j)
449  uhat.add(gamma(j), U[i]);
450  }
451  else
452  uhat *= omega;
453 
454  A.vmult(G[k], uhat);
455 
456  // Update G and U
457  // Orthogonalize G[k] to Q0,..,Q_{k-1} and update uhat
458  if (k > 0)
459  {
460  value_type alpha = Q[0] * G[k] / M(0, 0);
461  for (unsigned int i = 1; i < k; ++i)
462  {
463  const value_type alpha_old = alpha;
464  alpha = G[k].add_and_dot(-alpha, G[i - 1], Q[i]) / M(i, i);
465 
466  // update uhat every other iteration to reduce vector access
467  if (i % 2 == 1)
468  uhat.add(-alpha_old, U[i - 1], -alpha, U[i]);
469  }
470  M(k, k) = G[k].add_and_dot(-alpha, G[k - 1], Q[k]);
471  if (k % 2 == 1)
472  uhat.add(-alpha, U[k - 1]);
473  }
474  else
475  M(k, k) = G[k] * Q[k];
476 
477  U[k].swap(uhat);
478 
479  // Update kth column of M
480  for (unsigned int i = k + 1; i < s; ++i)
481  M(i, k) = Q[i] * G[k];
482 
483  // Orthogonalize r to Q0,...,Qk, update x
484  {
485  const value_type beta = phi(k) / M(k, k);
486  r.add(-beta, G[k]);
487  x.add(beta, U[k]);
488 
489  print_vectors(step, x, r, U[k]);
490 
491  // Check for early convergence. If so, store
492  // information in early_exit so that outer iteration
493  // is broken before recomputing the residual
494  res = r.l2_norm();
495  iteration_state = this->iteration_status(step, res, x);
496  if (iteration_state != SolverControl::iterate)
497  {
498  early_exit = true;
499  break;
500  }
501 
502  // Update phi
503  if (k + 1 < s)
504  {
505  for (unsigned int i = 0; i < k + 1; ++i)
506  phi(i) = 0.0;
507  for (unsigned int i = k + 1; i < s; ++i)
508  phi(i) -= beta * M(i, k);
509  }
510  }
511  }
512  if (early_exit == true)
513  break;
514 
515  // Update r and x
516  preconditioner.vmult(uhat, r);
517  A.vmult(v, uhat);
518 
519  omega = (v * r) / (v * v);
520 
521  res = std::sqrt(r.add_and_dot(-1.0 * omega, v, r));
522  x.add(omega, uhat);
523 
524  print_vectors(step, x, r, uhat);
525 
526  // Check for convergence
527  iteration_state = this->iteration_status(step, res, x);
528  if (iteration_state != SolverControl::iterate)
529  break;
530  }
531 
532  if (iteration_state != SolverControl::success)
533  AssertThrow(false, SolverControl::NoConvergence(step, res));
534 }
535 
536 
537 #endif // DOXYGEN
538 
540 
541 #endif
@ iterate
Continue iteration.
@ success
Stop iteration, goal reached.
virtual void print_vectors(const unsigned int step, const VectorType &x, const VectorType &r, const VectorType &d) const
SolverIDR(SolverControl &cn, const AdditionalData &data=AdditionalData())
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
SolverIDR(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
AdditionalData additional_data
Definition: solver_idr.h:182
virtual ~SolverIDR() override=default
Definition: vector.h:110
TmpVectors(const unsigned int s_param, VectorMemory< VectorType > &vmem)
VectorType & operator()(const unsigned int i, const VectorType &temp)
VectorMemory< VectorType > & mem
Definition: solver_idr.h:86
VectorType & operator[](const unsigned int i) const
std::vector< typename VectorMemory< VectorType >::Pointer > data
Definition: solver_idr.h:91
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1631
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
#define AssertThrow(cond, exc)
Definition: exceptions.h:1732
static const char U
static const char A
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition: operators.h:50
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)
long double gamma(const unsigned int n)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const unsigned int s
Definition: solver_idr.h:134
AdditionalData(const unsigned int s=2)
Definition: solver_idr.h:130