Reference documentation for deal.II version GIT 5983d193e2 2023-05-27 16:50: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\}}\)
precondition_block_base.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2020 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_precondition_block_base_h
17 #define dealii_precondition_block_base_h
18 
19 
20 #include <deal.II/base/config.h>
21 
26 
29 
30 #include <vector>
31 
33 
34 // Forward declarations
35 #ifndef DOXYGEN
36 template <typename number>
37 class FullMatrix;
38 template <typename number>
39 class Vector;
40 #endif
41 
57 template <typename number>
59 {
60 public:
65 
70  enum Inversion
71  {
84  svd
85  };
86 
91  Inversion method = gauss_jordan);
92 
97 
102  void
103  clear();
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
125  same_diagonal() const;
126 
130  bool
132 
136  bool
137  inverses_ready() const;
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> &
185  inverse(size_type i) const;
186 
190  const Householder<number> &
192 
198 
204 
208  const FullMatrix<number> &
209  diagonal(size_type i) const;
210 
216  void
217  log_statistics() const;
218 
223  std::size_t
225 
231 
238 
239 protected:
244 
245 private:
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 
302 template <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 
313 template <typename number>
314 inline 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 
331 template <typename number>
332 inline 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:
359  Assert(false, ExcNotImplemented());
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:
403  Assert(false, ExcNotImplemented());
404  }
405  }
406 }
407 
408 
409 template <typename number>
410 inline unsigned int
412 {
413  return n_diagonal_blocks;
414 }
415 
416 template <typename number>
417 template <typename number2>
418 inline 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:
440  Assert(false, ExcNotImplemented());
441  }
442 }
443 
444 
445 template <typename number>
446 template <typename number2>
447 inline 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:
469  Assert(false, ExcNotImplemented());
470  }
471 }
472 
473 
474 template <typename number>
475 inline 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 
486 template <typename number>
487 inline 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 
498 template <typename number>
499 inline 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 
510 template <typename number>
511 inline 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 
524 template <typename number>
525 inline 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 
538 template <typename number>
539 inline 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 
552 template <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 
566 template <typename number>
567 inline 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 
580 template <typename number>
581 inline bool
583 {
584  return var_same_diagonal;
585 }
586 
587 
588 template <typename number>
589 inline bool
591 {
592  return var_store_diagonals;
593 }
594 
595 
596 template <typename number>
597 inline void
599 {
600  var_inverses_ready = x;
601 }
602 
603 
604 template <typename number>
605 inline bool
607 {
608  return var_inverses_ready;
609 }
610 
611 
612 template <typename number>
613 inline 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  {
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  else if (inversion == gauss_jordan)
655  {}
656  else
657  {
658  Assert(false, ExcNotImplemented());
659  }
660 }
661 
662 
663 template <typename number>
664 inline std::size_t
666 {
667  std::size_t mem = sizeof(*this);
668  for (size_type i = 0; i < var_inverse_full.size(); ++i)
669  mem += MemoryConsumption::memory_consumption(var_inverse_full[i]);
670  for (size_type i = 0; i < var_diagonal.size(); ++i)
671  mem += MemoryConsumption::memory_consumption(var_diagonal[i]);
672  return mem;
673 }
674 
675 
677 
678 #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
types::global_dof_index size_type
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
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcDiagonalsNotStored()
#define DeclException0(Exception0)
Definition: exceptions.h:465
static ::ExceptionBase & ExcInverseNotAvailable()
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1855
LogStream deallog
Definition: logstream.cc:37
@ matrix
Contents is actually a matrix.
@ svd
Matrix contains singular value decomposition,.
@ inverse_svd
Matrix is the inverse of a singular value decomposition.
static const char U
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
std::string compress(const std::string &input)
Definition: utilities.cc:390
unsigned int global_dof_index
Definition: types.h:82