Reference documentation for deal.II version 9.6.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\}}\)
Loading...
Searching...
No Matches
solver_idr.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2019 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_solver_idr_h
16#define dealii_solver_idr_h
17
18
19#include <deal.II/base/config.h>
20
25
28#include <deal.II/lac/solver.h>
30
31#include <cmath>
32#include <random>
33
35
41namespace internal
42{
46 namespace SolverIDRImplementation
47 {
52 template <typename VectorType>
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
118template <typename VectorType = Vector<double>>
120class SolverIDR : public SolverBase<VectorType>
121{
122public:
127 {
131 explicit AdditionalData(const unsigned int s = 2)
132 : s(s)
133 {}
134
135 const unsigned int s;
136 };
137
143 const AdditionalData &data = AdditionalData());
144
150 const AdditionalData &data = AdditionalData());
151
155 virtual ~SolverIDR() override = default;
156
160 template <typename MatrixType, typename PreconditionerType>
164 void solve(const MatrixType &A,
165 VectorType &x,
166 const VectorType &b,
167 const PreconditionerType &preconditioner);
168
169protected:
175 virtual void
176 print_vectors(const unsigned int step,
177 const VectorType &x,
178 const VectorType &r,
179 const VectorType &d) const;
180
181private:
185 AdditionalData additional_data;
186};
187
189/*------------------------- Implementation ----------------------------*/
190
191#ifndef DOXYGEN
192
193
194namespace internal
195{
196 namespace SolverIDRImplementation
197 {
198 template <typename VectorType>
199 inline TmpVectors<VectorType>::TmpVectors(const unsigned int s_param,
201 : mem(vmem)
202 , data(s_param)
203 {}
204
205
206
207 template <typename VectorType>
208 inline VectorType &
209 TmpVectors<VectorType>::operator[](const unsigned int i) const
210 {
211 AssertIndexRange(i, data.size());
212
213 Assert(data[i] != nullptr, ExcNotInitialized());
214 return *data[i];
215 }
216
217
218
219 template <typename VectorType>
220 inline VectorType &
221 TmpVectors<VectorType>::operator()(const unsigned int i,
222 const VectorType &temp)
223 {
224 AssertIndexRange(i, data.size());
225 if (data[i] == nullptr)
226 {
227 data[i] = std::move(typename VectorMemory<VectorType>::Pointer(mem));
228 data[i]->reinit(temp, true);
229 }
230 return *data[i];
231 }
232
233
234
235 template <typename VectorType,
236 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
237 * = nullptr>
238 unsigned int
239 n_blocks(const VectorType &)
240 {
241 return 1;
242 }
243
244
245
246 template <typename VectorType,
247 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
248 nullptr>
249 unsigned int
250 n_blocks(const VectorType &vector)
251 {
252 return vector.n_blocks();
253 }
254
255
256
257 template <typename VectorType,
258 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
259 * = nullptr>
260 VectorType &
261 block(VectorType &vector, const unsigned int b)
262 {
263 AssertDimension(b, 0);
264 (void)b;
265 return vector;
266 }
267
268
269
270 template <typename VectorType,
271 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
272 nullptr>
273 typename VectorType::BlockType &
274 block(VectorType &vector, const unsigned int b)
275 {
276 return vector.block(b);
277 }
278
279 } // namespace SolverIDRImplementation
280} // namespace internal
281
282
283
284template <typename VectorType>
288 const AdditionalData &data)
289 : SolverBase<VectorType>(cn, mem)
290 , additional_data(data)
291{}
292
293
294
295template <typename VectorType>
297SolverIDR<VectorType>::SolverIDR(SolverControl &cn, const AdditionalData &data)
298 : SolverBase<VectorType>(cn)
299 , additional_data(data)
300{}
301
302
303
304template <typename VectorType>
306void SolverIDR<VectorType>::print_vectors(const unsigned int,
307 const VectorType &,
308 const VectorType &,
309 const VectorType &) const
310{}
311
312
313
314template <typename VectorType>
316template <typename MatrixType, typename PreconditionerType>
320void SolverIDR<VectorType>::solve(const MatrixType &A,
321 VectorType &x,
322 const VectorType &b,
323 const PreconditionerType &preconditioner)
324{
325 LogStream::Prefix prefix("IDR(s)");
326
328 unsigned int step = 0;
329
330 const unsigned int s = additional_data.s;
331
332 // Define temporary vectors which do not do not depend on s
333 typename VectorMemory<VectorType>::Pointer r_pointer(this->memory);
334 typename VectorMemory<VectorType>::Pointer v_pointer(this->memory);
335 typename VectorMemory<VectorType>::Pointer uhat_pointer(this->memory);
336
337 VectorType &r = *r_pointer;
338 VectorType &v = *v_pointer;
339 VectorType &uhat = *uhat_pointer;
340
341 r.reinit(x, true);
342 v.reinit(x, true);
343 uhat.reinit(x, true);
344
345 // Initial residual
346 A.vmult(r, x);
347 r.sadd(-1.0, 1.0, b);
348
349 using value_type = typename VectorType::value_type;
350 using real_type = typename numbers::NumberTraits<value_type>::real_type;
351
352 // Check for convergent initial guess
353 real_type res = r.l2_norm();
354 iteration_state = this->iteration_status(step, res, x);
355 if (iteration_state == SolverControl::success)
356 return;
357
358 // Initialize sets of vectors/matrices whose size dependent on s
363
364 // Random number generator for vector entries of
365 // Q (normal distribution, mean=0 sigma=1)
366 std::mt19937 rng;
367 std::normal_distribution<> normal_distribution(0.0, 1.0);
368 for (unsigned int i = 0; i < s; ++i)
369 {
370 // Initialize vectors
371 G(i, x);
372 U(i, x);
373
374 // Compute random set of s orthonormalized vectors Q
375 // Note: the first vector is chosen to be the initial
376 // residual to match BiCGStab (as is done in comparisons
377 // with BiCGStab in the papers listed in the documentation
378 // of this function)
379 VectorType &tmp_q = Q(i, x);
380 if (i != 0)
381 {
382 for (unsigned int b = 0;
383 b < internal::SolverIDRImplementation::n_blocks(tmp_q);
384 ++b)
385 for (auto indx : internal::SolverIDRImplementation::block(tmp_q, b)
386 .locally_owned_elements())
387 internal::SolverIDRImplementation::block(tmp_q, b)(indx) =
388 normal_distribution(rng);
389 tmp_q.compress(VectorOperation::insert);
390 }
391 else
392 tmp_q = r;
393
394 for (unsigned int j = 0; j < i; ++j)
395 {
396 v = Q[j];
397 v *= (v * tmp_q) / (tmp_q * tmp_q);
398 tmp_q.add(-1.0, v);
399 }
400
401 if (i != 0)
402 tmp_q *= 1.0 / tmp_q.l2_norm();
403
404 M(i, i) = 1.;
405 }
406
407 value_type omega = 1.;
408
409 bool early_exit = false;
410
411 // Outer iteration
412 while (iteration_state == SolverControl::iterate)
413 {
414 ++step;
415
416 // Compute phi
417 Vector<value_type> phi(s);
418 for (unsigned int i = 0; i < s; ++i)
419 phi(i) = Q[i] * r;
420
421 // Inner iteration over s
422 for (unsigned int k = 0; k < s; ++k)
423 {
424 // Solve M(k:s)*gamma = phi(k:s)
426 {
427 Vector<value_type> phik(s - k);
428 FullMatrix<value_type> Mk(s - k, s - k);
429 std::vector<unsigned int> indices;
430 unsigned int j = 0;
431 for (unsigned int i = k; i < s; ++i, ++j)
432 {
433 indices.push_back(i);
434 phik(j) = phi(i);
435 }
436 Mk.extract_submatrix_from(M, indices, indices);
437
438 FullMatrix<value_type> Mk_inv(s - k, s - k);
439 Mk_inv.invert(Mk);
440 Mk_inv.vmult(gamma, phik);
441 }
442
443 v = r;
444
445 if (step > 1)
446 {
447 for (unsigned int i = k, j = 0; i < s; ++i, ++j)
448 v.add(-gamma(j), G[i]);
449 }
450
451 preconditioner.vmult(uhat, v);
452
453 if (step > 1)
454 {
455 uhat.sadd(omega, gamma(0), U[k]);
456 for (unsigned int i = k + 1, j = 1; i < s; ++i, ++j)
457 uhat.add(gamma(j), U[i]);
458 }
459 else
460 uhat *= omega;
461
462 A.vmult(G[k], uhat);
463
464 // Update G and U
465 // Orthogonalize G[k] to Q0,..,Q_{k-1} and update uhat
466 if (k > 0)
467 {
468 value_type alpha = Q[0] * G[k] / M(0, 0);
469 for (unsigned int i = 1; i < k; ++i)
470 {
471 const value_type alpha_old = alpha;
472 alpha = G[k].add_and_dot(-alpha, G[i - 1], Q[i]) / M(i, i);
473
474 // update uhat every other iteration to reduce vector access
475 if (i % 2 == 1)
476 uhat.add(-alpha_old, U[i - 1], -alpha, U[i]);
477 }
478 M(k, k) = G[k].add_and_dot(-alpha, G[k - 1], Q[k]);
479 if (k % 2 == 1)
480 uhat.add(-alpha, U[k - 1]);
481 }
482 else
483 M(k, k) = G[k] * Q[k];
484
485 U[k].swap(uhat);
486
487 // Update kth column of M
488 for (unsigned int i = k + 1; i < s; ++i)
489 M(i, k) = Q[i] * G[k];
490
491 // Orthogonalize r to Q0,...,Qk, update x
492 {
493 const value_type beta = phi(k) / M(k, k);
494 r.add(-beta, G[k]);
495 x.add(beta, U[k]);
496
497 print_vectors(step, x, r, U[k]);
498
499 // Check for early convergence. If so, store
500 // information in early_exit so that outer iteration
501 // is broken before recomputing the residual
502 res = r.l2_norm();
503 iteration_state = this->iteration_status(step, res, x);
504 if (iteration_state != SolverControl::iterate)
505 {
506 early_exit = true;
507 break;
508 }
509
510 // Update phi
511 if (k + 1 < s)
512 {
513 for (unsigned int i = 0; i < k + 1; ++i)
514 phi(i) = 0.0;
515 for (unsigned int i = k + 1; i < s; ++i)
516 phi(i) -= beta * M(i, k);
517 }
518 }
519 }
520 if (early_exit == true)
521 break;
522
523 // Update r and x
524 preconditioner.vmult(uhat, r);
525 A.vmult(v, uhat);
526
527 omega = (v * r) / (v * v);
528
529 res = std::sqrt(r.add_and_dot(-1.0 * omega, v, r));
530 x.add(omega, uhat);
531
532 print_vectors(step, x, r, uhat);
533
534 // Check for convergence
535 iteration_state = this->iteration_status(step, res, x);
536 if (iteration_state != SolverControl::iterate)
537 break;
538 }
539
540 if (iteration_state != SolverControl::success)
541 AssertThrow(false, SolverControl::NoConvergence(step, res));
542}
543
544
545#endif // DOXYGEN
546
548
549#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())
virtual ~SolverIDR() override=default
TmpVectors(const unsigned int s_param, VectorMemory< VectorType > &vmem)
VectorType & operator()(const unsigned int i, const VectorType &temp)
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:503
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:177
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcNotInitialized()
#define AssertThrow(cond, exc)
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:135
AdditionalData(const unsigned int s=2)
Definition solver_idr.h:131