deal.II version GIT relicensing-2206-gaa53ff9447 2024-12-02 09:10:00+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
26
29
30#include <vector>
31
33
34// Forward declarations
35#ifndef DOXYGEN
36template <typename number>
37class FullMatrix;
38template <typename number>
39class Vector;
40#endif
41
57template <typename number>
59{
60public:
65
86
91 Inversion method = gauss_jordan);
92
97
102 void
104
109 void
110 reinit(unsigned int nblocks,
111 size_type blocksize,
112 bool compress,
113 Inversion method = gauss_jordan);
114
118 void
119 inverses_computed(bool are_they);
120
124 bool
126
130 bool
132
136 bool
138
142 unsigned int
143 size() const;
144
148 template <typename number2>
149 void
151 Vector<number2> &dst,
152 const Vector<number2> &src) const;
153
157 template <typename number2>
158 void
160 Vector<number2> &dst,
161 const Vector<number2> &src) const;
162
168
174
180
184 const FullMatrix<number> &
186
190 const Householder<number> &
192
198
204
208 const FullMatrix<number> &
210
216 void
218
223 std::size_t
225
231
238
239protected:
244
245private:
249 unsigned int n_diagonal_blocks;
250
257 std::vector<FullMatrix<number>> var_inverse_full;
258
265 std::vector<Householder<number>> var_inverse_householder;
266
273 std::vector<LAPACKFullMatrix<number>> var_inverse_svd;
274
280 std::vector<FullMatrix<number>> var_diagonal;
281
282
287
292
298};
299
300//----------------------------------------------------------------------//
301
302template <typename number>
304 Inversion method)
305 : inversion(method)
306 , n_diagonal_blocks(0)
307 , var_store_diagonals(store)
308 , var_same_diagonal(false)
309 , var_inverses_ready(false)
310{}
311
312
313template <typename number>
314inline void
316{
317 if (var_inverse_full.size() != 0)
318 var_inverse_full.erase(var_inverse_full.begin(), var_inverse_full.end());
319 if (var_inverse_householder.size() != 0)
320 var_inverse_householder.erase(var_inverse_householder.begin(),
321 var_inverse_householder.end());
322 if (var_inverse_svd.size() != 0)
323 var_inverse_svd.erase(var_inverse_svd.begin(), var_inverse_svd.end());
324 if (var_diagonal.size() != 0)
325 var_diagonal.erase(var_diagonal.begin(), var_diagonal.end());
326 var_same_diagonal = false;
327 var_inverses_ready = false;
328 n_diagonal_blocks = 0;
329}
330
331template <typename number>
332inline void
334 size_type b,
335 bool compress,
336 Inversion method)
337{
338 inversion = method;
339 var_same_diagonal = compress;
340 var_inverses_ready = false;
341 n_diagonal_blocks = n;
342
343 if (compress)
344 {
345 switch (inversion)
346 {
347 case gauss_jordan:
348 var_inverse_full.resize(1);
349 var_inverse_full[0].reinit(b, b);
350 break;
351 case householder:
352 var_inverse_householder.resize(1);
353 break;
354 case svd:
355 var_inverse_svd.resize(1);
356 var_inverse_svd[0].reinit(b, b);
357 break;
358 default:
360 }
361
362 if (store_diagonals())
363 {
364 var_diagonal.resize(1);
365 var_diagonal[0].reinit(b, b);
366 }
367 }
368 else
369 {
370 // set the arrays to the right
371 // size. we could do it like this:
372 // var_inverse = vector<>(nblocks,FullMatrix<>())
373 // but this would involve copying many
374 // FullMatrix objects.
375 //
376 // the following is a neat trick which
377 // avoids copying
378 if (store_diagonals())
379 {
380 std::vector<FullMatrix<number>> tmp(n, FullMatrix<number>(b));
381 var_diagonal.swap(tmp);
382 }
383
384 switch (inversion)
385 {
386 case gauss_jordan:
387 {
388 std::vector<FullMatrix<number>> tmp(n, FullMatrix<number>(b));
389 var_inverse_full.swap(tmp);
390 break;
391 }
392 case householder:
393 var_inverse_householder.resize(n);
394 break;
395 case svd:
396 {
397 std::vector<LAPACKFullMatrix<number>> tmp(
399 var_inverse_svd.swap(tmp);
400 break;
401 }
402 default:
404 }
405 }
406}
407
408
409template <typename number>
410inline unsigned int
412{
413 return n_diagonal_blocks;
414}
415
416template <typename number>
417template <typename number2>
418inline void
420 Vector<number2> &dst,
421 const Vector<number2> &src) const
422{
423 const size_type ii = same_diagonal() ? 0U : i;
424
425 switch (inversion)
426 {
427 case gauss_jordan:
428 AssertIndexRange(ii, var_inverse_full.size());
429 var_inverse_full[ii].vmult(dst, src);
430 break;
431 case householder:
432 AssertIndexRange(ii, var_inverse_householder.size());
433 var_inverse_householder[ii].vmult(dst, src);
434 break;
435 case svd:
436 AssertIndexRange(ii, var_inverse_svd.size());
437 var_inverse_svd[ii].vmult(dst, src);
438 break;
439 default:
441 }
442}
443
444
445template <typename number>
446template <typename number2>
447inline void
449 Vector<number2> &dst,
450 const Vector<number2> &src) const
451{
452 const size_type ii = same_diagonal() ? 0U : i;
453
454 switch (inversion)
455 {
456 case gauss_jordan:
457 AssertIndexRange(ii, var_inverse_full.size());
458 var_inverse_full[ii].Tvmult(dst, src);
459 break;
460 case householder:
461 AssertIndexRange(ii, var_inverse_householder.size());
462 var_inverse_householder[ii].Tvmult(dst, src);
463 break;
464 case svd:
465 AssertIndexRange(ii, var_inverse_svd.size());
466 var_inverse_svd[ii].Tvmult(dst, src);
467 break;
468 default:
470 }
471}
472
473
474template <typename number>
475inline const FullMatrix<number> &
477{
478 if (same_diagonal())
479 return var_inverse_full[0];
480
481 AssertIndexRange(i, var_inverse_full.size());
482 return var_inverse_full[i];
483}
484
485
486template <typename number>
487inline const Householder<number> &
489{
490 if (same_diagonal())
491 return var_inverse_householder[0];
492
493 AssertIndexRange(i, var_inverse_householder.size());
494 return var_inverse_householder[i];
495}
496
497
498template <typename number>
499inline const LAPACKFullMatrix<number> &
501{
502 if (same_diagonal())
503 return var_inverse_svd[0];
504
505 AssertIndexRange(i, var_inverse_svd.size());
506 return var_inverse_svd[i];
507}
508
509
510template <typename number>
511inline const FullMatrix<number> &
513{
514 Assert(store_diagonals(), ExcDiagonalsNotStored());
515
516 if (same_diagonal())
517 return var_diagonal[0];
518
519 AssertIndexRange(i, var_diagonal.size());
520 return var_diagonal[i];
521}
522
523
524template <typename number>
525inline FullMatrix<number> &
527{
528 Assert(var_inverse_full.size() != 0, ExcInverseNotAvailable());
529
530 if (same_diagonal())
531 return var_inverse_full[0];
532
533 AssertIndexRange(i, var_inverse_full.size());
534 return var_inverse_full[i];
535}
536
537
538template <typename number>
539inline Householder<number> &
541{
542 Assert(var_inverse_householder.size() != 0, ExcInverseNotAvailable());
543
544 if (same_diagonal())
545 return var_inverse_householder[0];
546
547 AssertIndexRange(i, var_inverse_householder.size());
548 return var_inverse_householder[i];
549}
550
551
552template <typename number>
555{
556 Assert(var_inverse_svd.size() != 0, ExcInverseNotAvailable());
557
558 if (same_diagonal())
559 return var_inverse_svd[0];
560
561 AssertIndexRange(i, var_inverse_svd.size());
562 return var_inverse_svd[i];
563}
564
565
566template <typename number>
567inline FullMatrix<number> &
569{
570 Assert(store_diagonals(), ExcDiagonalsNotStored());
571
572 if (same_diagonal())
573 return var_diagonal[0];
574
575 AssertIndexRange(i, var_diagonal.size());
576 return var_diagonal[i];
577}
578
579
580template <typename number>
581inline bool
583{
584 return var_same_diagonal;
585}
586
587
588template <typename number>
589inline bool
591{
592 return var_store_diagonals;
593}
594
595
596template <typename number>
597inline void
599{
600 var_inverses_ready = x;
601}
602
603
604template <typename number>
605inline bool
607{
608 return var_inverses_ready;
609}
610
611
612template <typename number>
613inline void
615{
616 deallog << "PreconditionBlockBase: " << size() << " blocks; ";
617
618 if (inversion == svd)
619 {
620 unsigned int kermin = 100000000, kermax = 0;
621 double sigmin = 1.e300, sigmax = -1.e300;
622 double kappamin = 1.e300, kappamax = -1.e300;
623
624 for (size_type b = 0; b < size(); ++b)
625 {
626 const LAPACKFullMatrix<number> &matrix = inverse_svd(b);
627 size_type k = 1;
628 while (k <= matrix.n_cols() &&
629 matrix.singular_value(matrix.n_cols() - k) == 0)
630 ++k;
631 const double s0 = matrix.singular_value(0);
632 const double sm = matrix.singular_value(matrix.n_cols() - k);
633 const double co = sm / s0;
634
635 if (kermin > k)
636 kermin = k - 1;
637 if (kermax < k)
638 kermax = k - 1;
639 if (s0 < sigmin)
640 sigmin = s0;
641 if (sm > sigmax)
642 sigmax = sm;
643 if (co < kappamin)
644 kappamin = co;
645 if (co > kappamax)
646 kappamax = co;
647 }
648 deallog << "dim ker [" << kermin << ':' << kermax << "] sigma [" << sigmin
649 << ':' << sigmax << "] kappa [" << kappamin << ':' << kappamax
650 << ']' << std::endl;
651 }
652 else if (inversion == householder)
653 {
654 }
655 else if (inversion == gauss_jordan)
656 {
657 }
658 else
659 {
661 }
662}
663
664
665template <typename number>
666inline std::size_t
668{
669 std::size_t mem = sizeof(*this);
670 for (size_type i = 0; i < var_inverse_full.size(); ++i)
671 mem += MemoryConsumption::memory_consumption(var_inverse_full[i]);
672 for (size_type i = 0; i < var_diagonal.size(); ++i)
673 mem += MemoryConsumption::memory_consumption(var_diagonal[i]);
674 return mem;
675}
676
677
679
680#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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
#define DeclException0(Exception0)
Definition exceptions.h:466
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::size_t size
Definition mpi.cc:734
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