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