Reference documentation for deal.II version GIT relicensing-233-g802318d791 2024-03-28 20: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
precondition_block_base.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2010 - 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_precondition_block_base_h
16#define dealii_precondition_block_base_h
17
18
19#include <deal.II/base/config.h>
20
25
28
29#include <vector>
30
32
33// Forward declarations
34#ifndef DOXYGEN
35template <typename number>
36class FullMatrix;
37template <typename number>
38class Vector;
39#endif
40
56template <typename number>
58{
59public:
64
85
90 Inversion method = gauss_jordan);
91
96
101 void
103
108 void
109 reinit(unsigned int nblocks,
110 size_type blocksize,
111 bool compress,
112 Inversion method = gauss_jordan);
113
117 void
118 inverses_computed(bool are_they);
119
123 bool
125
129 bool
131
135 bool
137
141 unsigned int
142 size() const;
143
147 template <typename number2>
148 void
150 Vector<number2> &dst,
151 const Vector<number2> &src) const;
152
156 template <typename number2>
157 void
159 Vector<number2> &dst,
160 const Vector<number2> &src) const;
161
167
173
179
183 const FullMatrix<number> &
185
189 const Householder<number> &
191
197
203
207 const FullMatrix<number> &
209
215 void
217
222 std::size_t
224
230
237
238protected:
243
244private:
248 unsigned int n_diagonal_blocks;
249
256 std::vector<FullMatrix<number>> var_inverse_full;
257
264 std::vector<Householder<number>> var_inverse_householder;
265
272 std::vector<LAPACKFullMatrix<number>> var_inverse_svd;
273
279 std::vector<FullMatrix<number>> var_diagonal;
280
281
286
291
297};
298
299//----------------------------------------------------------------------//
300
301template <typename number>
303 Inversion method)
304 : inversion(method)
305 , n_diagonal_blocks(0)
306 , var_store_diagonals(store)
307 , var_same_diagonal(false)
308 , var_inverses_ready(false)
309{}
310
311
312template <typename number>
313inline void
315{
316 if (var_inverse_full.size() != 0)
317 var_inverse_full.erase(var_inverse_full.begin(), var_inverse_full.end());
318 if (var_inverse_householder.size() != 0)
319 var_inverse_householder.erase(var_inverse_householder.begin(),
320 var_inverse_householder.end());
321 if (var_inverse_svd.size() != 0)
322 var_inverse_svd.erase(var_inverse_svd.begin(), var_inverse_svd.end());
323 if (var_diagonal.size() != 0)
324 var_diagonal.erase(var_diagonal.begin(), var_diagonal.end());
325 var_same_diagonal = false;
326 var_inverses_ready = false;
327 n_diagonal_blocks = 0;
328}
329
330template <typename number>
331inline void
333 size_type b,
334 bool compress,
335 Inversion method)
336{
337 inversion = method;
338 var_same_diagonal = compress;
339 var_inverses_ready = false;
340 n_diagonal_blocks = n;
341
342 if (compress)
343 {
344 switch (inversion)
345 {
346 case gauss_jordan:
347 var_inverse_full.resize(1);
348 var_inverse_full[0].reinit(b, b);
349 break;
350 case householder:
351 var_inverse_householder.resize(1);
352 break;
353 case svd:
354 var_inverse_svd.resize(1);
355 var_inverse_svd[0].reinit(b, b);
356 break;
357 default:
359 }
360
361 if (store_diagonals())
362 {
363 var_diagonal.resize(1);
364 var_diagonal[0].reinit(b, b);
365 }
366 }
367 else
368 {
369 // set the arrays to the right
370 // size. we could do it like this:
371 // var_inverse = vector<>(nblocks,FullMatrix<>())
372 // but this would involve copying many
373 // FullMatrix objects.
374 //
375 // the following is a neat trick which
376 // avoids copying
377 if (store_diagonals())
378 {
379 std::vector<FullMatrix<number>> tmp(n, FullMatrix<number>(b));
380 var_diagonal.swap(tmp);
381 }
382
383 switch (inversion)
384 {
385 case gauss_jordan:
386 {
387 std::vector<FullMatrix<number>> tmp(n, FullMatrix<number>(b));
388 var_inverse_full.swap(tmp);
389 break;
390 }
391 case householder:
392 var_inverse_householder.resize(n);
393 break;
394 case svd:
395 {
396 std::vector<LAPACKFullMatrix<number>> tmp(
398 var_inverse_svd.swap(tmp);
399 break;
400 }
401 default:
403 }
404 }
405}
406
407
408template <typename number>
409inline unsigned int
411{
412 return n_diagonal_blocks;
413}
414
415template <typename number>
416template <typename number2>
417inline void
419 Vector<number2> &dst,
420 const Vector<number2> &src) const
421{
422 const size_type ii = same_diagonal() ? 0U : i;
423
424 switch (inversion)
425 {
426 case gauss_jordan:
427 AssertIndexRange(ii, var_inverse_full.size());
428 var_inverse_full[ii].vmult(dst, src);
429 break;
430 case householder:
431 AssertIndexRange(ii, var_inverse_householder.size());
432 var_inverse_householder[ii].vmult(dst, src);
433 break;
434 case svd:
435 AssertIndexRange(ii, var_inverse_svd.size());
436 var_inverse_svd[ii].vmult(dst, src);
437 break;
438 default:
440 }
441}
442
443
444template <typename number>
445template <typename number2>
446inline void
448 Vector<number2> &dst,
449 const Vector<number2> &src) const
450{
451 const size_type ii = same_diagonal() ? 0U : i;
452
453 switch (inversion)
454 {
455 case gauss_jordan:
456 AssertIndexRange(ii, var_inverse_full.size());
457 var_inverse_full[ii].Tvmult(dst, src);
458 break;
459 case householder:
460 AssertIndexRange(ii, var_inverse_householder.size());
461 var_inverse_householder[ii].Tvmult(dst, src);
462 break;
463 case svd:
464 AssertIndexRange(ii, var_inverse_svd.size());
465 var_inverse_svd[ii].Tvmult(dst, src);
466 break;
467 default:
469 }
470}
471
472
473template <typename number>
474inline const FullMatrix<number> &
476{
477 if (same_diagonal())
478 return var_inverse_full[0];
479
480 AssertIndexRange(i, var_inverse_full.size());
481 return var_inverse_full[i];
482}
483
484
485template <typename number>
486inline const Householder<number> &
488{
489 if (same_diagonal())
490 return var_inverse_householder[0];
491
492 AssertIndexRange(i, var_inverse_householder.size());
493 return var_inverse_householder[i];
494}
495
496
497template <typename number>
498inline const LAPACKFullMatrix<number> &
500{
501 if (same_diagonal())
502 return var_inverse_svd[0];
503
504 AssertIndexRange(i, var_inverse_svd.size());
505 return var_inverse_svd[i];
506}
507
508
509template <typename number>
510inline const FullMatrix<number> &
512{
513 Assert(store_diagonals(), ExcDiagonalsNotStored());
514
515 if (same_diagonal())
516 return var_diagonal[0];
517
518 AssertIndexRange(i, var_diagonal.size());
519 return var_diagonal[i];
520}
521
522
523template <typename number>
524inline FullMatrix<number> &
526{
527 Assert(var_inverse_full.size() != 0, ExcInverseNotAvailable());
528
529 if (same_diagonal())
530 return var_inverse_full[0];
531
532 AssertIndexRange(i, var_inverse_full.size());
533 return var_inverse_full[i];
534}
535
536
537template <typename number>
538inline Householder<number> &
540{
541 Assert(var_inverse_householder.size() != 0, ExcInverseNotAvailable());
542
543 if (same_diagonal())
544 return var_inverse_householder[0];
545
546 AssertIndexRange(i, var_inverse_householder.size());
547 return var_inverse_householder[i];
548}
549
550
551template <typename number>
554{
555 Assert(var_inverse_svd.size() != 0, ExcInverseNotAvailable());
556
557 if (same_diagonal())
558 return var_inverse_svd[0];
559
560 AssertIndexRange(i, var_inverse_svd.size());
561 return var_inverse_svd[i];
562}
563
564
565template <typename number>
566inline FullMatrix<number> &
568{
569 Assert(store_diagonals(), ExcDiagonalsNotStored());
570
571 if (same_diagonal())
572 return var_diagonal[0];
573
574 AssertIndexRange(i, var_diagonal.size());
575 return var_diagonal[i];
576}
577
578
579template <typename number>
580inline bool
582{
583 return var_same_diagonal;
584}
585
586
587template <typename number>
588inline bool
590{
591 return var_store_diagonals;
592}
593
594
595template <typename number>
596inline void
598{
599 var_inverses_ready = x;
600}
601
602
603template <typename number>
604inline bool
606{
607 return var_inverses_ready;
608}
609
610
611template <typename number>
612inline void
614{
615 deallog << "PreconditionBlockBase: " << size() << " blocks; ";
616
617 if (inversion == svd)
618 {
619 unsigned int kermin = 100000000, kermax = 0;
620 double sigmin = 1.e300, sigmax = -1.e300;
621 double kappamin = 1.e300, kappamax = -1.e300;
622
623 for (size_type b = 0; b < size(); ++b)
624 {
625 const LAPACKFullMatrix<number> &matrix = inverse_svd(b);
626 size_type k = 1;
627 while (k <= matrix.n_cols() &&
628 matrix.singular_value(matrix.n_cols() - k) == 0)
629 ++k;
630 const double s0 = matrix.singular_value(0);
631 const double sm = matrix.singular_value(matrix.n_cols() - k);
632 const double co = sm / s0;
633
634 if (kermin > k)
635 kermin = k - 1;
636 if (kermax < k)
637 kermax = k - 1;
638 if (s0 < sigmin)
639 sigmin = s0;
640 if (sm > sigmax)
641 sigmax = sm;
642 if (co < kappamin)
643 kappamin = co;
644 if (co > kappamax)
645 kappamax = co;
646 }
647 deallog << "dim ker [" << kermin << ':' << kermax << "] sigma [" << sigmin
648 << ':' << sigmax << "] kappa [" << kappamin << ':' << kappamax
649 << ']' << std::endl;
650 }
651 else if (inversion == householder)
652 {
653 }
654 else if (inversion == gauss_jordan)
655 {
656 }
657 else
658 {
660 }
661}
662
663
664template <typename number>
665inline std::size_t
667{
668 std::size_t mem = sizeof(*this);
669 for (size_type i = 0; i < var_inverse_full.size(); ++i)
670 mem += MemoryConsumption::memory_consumption(var_inverse_full[i]);
671 for (size_type i = 0; i < var_diagonal.size(); ++i)
672 mem += MemoryConsumption::memory_consumption(var_diagonal[i]);
673 return mem;
674}
675
676
678
679#endif
std::vector< FullMatrix< number > > var_diagonal
std::size_t memory_consumption() const
std::vector< LAPACKFullMatrix< number > > var_inverse_svd
void inverses_computed(bool are_they)
void reinit(unsigned int nblocks, size_type blocksize, bool compress, Inversion method=gauss_jordan)
void inverse_vmult(size_type i, Vector< number2 > &dst, const Vector< number2 > &src) const
PreconditionBlockBase(bool store_diagonals=false, Inversion method=gauss_jordan)
Householder< number > & inverse_householder(size_type i)
std::vector< FullMatrix< number > > var_inverse_full
const FullMatrix< number > & diagonal(size_type i) const
std::vector< Householder< number > > var_inverse_householder
FullMatrix< number > & inverse(size_type i)
FullMatrix< number > & diagonal(size_type i)
const Householder< number > & inverse_householder(size_type i) const
const FullMatrix< number > & inverse(size_type i) const
void inverse_Tvmult(size_type i, Vector< number2 > &dst, const Vector< number2 > &src) const
LAPACKFullMatrix< number > & inverse_svd(size_type i)
const LAPACKFullMatrix< number > & inverse_svd(size_type i) const
~PreconditionBlockBase()=default
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define DeclException0(Exception0)
Definition exceptions.h:471
static ::ExceptionBase & ExcDiagonalsNotStored()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInverseNotAvailable()
#define DEAL_II_NOT_IMPLEMENTED()
LogStream deallog
Definition logstream.cc:36
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
unsigned int global_dof_index
Definition types.h:81