Reference documentation for deal.II version 9.5.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// 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
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 <class VectorType = Vector<double>>
119class SolverIDR : public SolverBase<VectorType>
120{
121public:
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
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
166protected:
172 virtual void
173 print_vectors(const unsigned int step,
174 const VectorType & x,
175 const VectorType & r,
176 const VectorType & d) const;
177
178private:
183};
184
186/*------------------------- Implementation ----------------------------*/
187
188#ifndef DOXYGEN
189
190
191namespace internal
192{
193 namespace SolverIDRImplementation
194 {
195 template <class VectorType>
196 inline TmpVectors<VectorType>::TmpVectors(const unsigned int s_param,
198 : mem(vmem)
199 , data(s_param)
200 {}
201
202
203
204 template <class 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 <class 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
281template <class VectorType>
284 const AdditionalData & data)
285 : SolverBase<VectorType>(cn, mem)
286 , additional_data(data)
287{}
288
289
290
291template <class VectorType>
292SolverIDR<VectorType>::SolverIDR(SolverControl &cn, const AdditionalData &data)
293 : SolverBase<VectorType>(cn)
294 , additional_data(data)
295{}
296
297
298
299template <class VectorType>
300void
301SolverIDR<VectorType>::print_vectors(const unsigned int,
302 const VectorType &,
303 const VectorType &,
304 const VectorType &) const
305{}
306
307
308
309template <class VectorType>
310template <typename MatrixType, typename PreconditionerType>
311void
312SolverIDR<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
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;
375 b < internal::SolverIDRImplementation::n_blocks(tmp_q);
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)
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
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#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:134
AdditionalData(const unsigned int s=2)
Definition solver_idr.h:130