Reference documentation for deal.II version 9.2.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\}}\)
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 
60 template <typename number>
62 {
63 public:
68 
73  enum Inversion
74  {
88  };
89 
94  Inversion method = gauss_jordan);
95 
99  ~PreconditionBlockBase() = default;
100 
105  void
106  clear();
107 
112  void
113  reinit(unsigned int nblocks,
114  size_type blocksize,
115  bool compress,
116  Inversion method = gauss_jordan);
117 
121  void
122  inverses_computed(bool are_they);
123 
127  bool
128  same_diagonal() const;
129 
133  bool
134  store_diagonals() const;
135 
139  bool
140  inverses_ready() const;
141 
145  unsigned int
146  size() const;
147 
151  template <typename number2>
152  void
154  Vector<number2> & dst,
155  const Vector<number2> &src) const;
156 
160  template <typename number2>
161  void
163  Vector<number2> & dst,
164  const Vector<number2> &src) const;
165 
170  inverse(size_type i);
171 
177 
183 
187  const FullMatrix<number> &
188  inverse(size_type i) const;
189 
193  const Householder<number> &
195 
200  inverse_svd(size_type i) const;
201 
206  diagonal(size_type i);
207 
211  const FullMatrix<number> &
212  diagonal(size_type i) const;
213 
219  void
220  log_statistics() const;
221 
226  std::size_t
227  memory_consumption() const;
228 
234 
241 
242 protected:
247 
248 private:
252  unsigned int n_diagonal_blocks;
253 
260  std::vector<FullMatrix<number>> var_inverse_full;
261 
268  std::vector<Householder<number>> var_inverse_householder;
269 
276  std::vector<LAPACKFullMatrix<number>> var_inverse_svd;
277 
283  std::vector<FullMatrix<number>> var_diagonal;
284 
285 
290 
295 
301 };
302 
303 //----------------------------------------------------------------------//
304 
305 template <typename number>
307  Inversion method)
308  : inversion(method)
309  , n_diagonal_blocks(0)
310  , var_store_diagonals(store)
311  , var_same_diagonal(false)
312  , var_inverses_ready(false)
313 {}
314 
315 
316 template <typename number>
317 inline void
319 {
320  if (var_inverse_full.size() != 0)
321  var_inverse_full.erase(var_inverse_full.begin(), var_inverse_full.end());
322  if (var_inverse_householder.size() != 0)
325  if (var_inverse_svd.size() != 0)
326  var_inverse_svd.erase(var_inverse_svd.begin(), var_inverse_svd.end());
327  if (var_diagonal.size() != 0)
328  var_diagonal.erase(var_diagonal.begin(), var_diagonal.end());
329  var_same_diagonal = false;
330  var_inverses_ready = false;
331  n_diagonal_blocks = 0;
332 }
333 
334 template <typename number>
335 inline void
337  size_type b,
338  bool compress,
339  Inversion method)
340 {
341  inversion = method;
343  var_inverses_ready = false;
344  n_diagonal_blocks = n;
345 
346  if (compress)
347  {
348  switch (inversion)
349  {
350  case gauss_jordan:
351  var_inverse_full.resize(1);
352  var_inverse_full[0].reinit(b, b);
353  break;
354  case householder:
355  var_inverse_householder.resize(1);
356  break;
357  case svd:
358  var_inverse_svd.resize(1);
359  var_inverse_svd[0].reinit(b, b);
360  break;
361  default:
362  Assert(false, ExcNotImplemented());
363  }
364 
365  if (store_diagonals())
366  {
367  var_diagonal.resize(1);
368  var_diagonal[0].reinit(b, b);
369  }
370  }
371  else
372  {
373  // set the arrays to the right
374  // size. we could do it like this:
375  // var_inverse = vector<>(nblocks,FullMatrix<>())
376  // but this would involve copying many
377  // FullMatrix objects.
378  //
379  // the following is a neat trick which
380  // avoids copying
381  if (store_diagonals())
382  {
383  std::vector<FullMatrix<number>> tmp(n, FullMatrix<number>(b));
384  var_diagonal.swap(tmp);
385  }
386 
387  switch (inversion)
388  {
389  case gauss_jordan:
390  {
391  std::vector<FullMatrix<number>> tmp(n, FullMatrix<number>(b));
392  var_inverse_full.swap(tmp);
393  break;
394  }
395  case householder:
396  var_inverse_householder.resize(n);
397  break;
398  case svd:
399  {
400  std::vector<LAPACKFullMatrix<number>> tmp(
402  var_inverse_svd.swap(tmp);
403  break;
404  }
405  default:
406  Assert(false, ExcNotImplemented());
407  }
408  }
409 }
410 
411 
412 template <typename number>
413 inline unsigned int
415 {
416  return n_diagonal_blocks;
417 }
418 
419 template <typename number>
420 template <typename number2>
421 inline void
423  Vector<number2> & dst,
424  const Vector<number2> &src) const
425 {
426  const size_type ii = same_diagonal() ? 0U : i;
427 
428  switch (inversion)
429  {
430  case gauss_jordan:
432  var_inverse_full[ii].vmult(dst, src);
433  break;
434  case householder:
436  var_inverse_householder[ii].vmult(dst, src);
437  break;
438  case svd:
439  AssertIndexRange(ii, var_inverse_svd.size());
440  var_inverse_svd[ii].vmult(dst, src);
441  break;
442  default:
443  Assert(false, ExcNotImplemented());
444  }
445 }
446 
447 
448 template <typename number>
449 template <typename number2>
450 inline void
452  Vector<number2> & dst,
453  const Vector<number2> &src) const
454 {
455  const size_type ii = same_diagonal() ? 0U : i;
456 
457  switch (inversion)
458  {
459  case gauss_jordan:
461  var_inverse_full[ii].Tvmult(dst, src);
462  break;
463  case householder:
465  var_inverse_householder[ii].Tvmult(dst, src);
466  break;
467  case svd:
468  AssertIndexRange(ii, var_inverse_svd.size());
469  var_inverse_svd[ii].Tvmult(dst, src);
470  break;
471  default:
472  Assert(false, ExcNotImplemented());
473  }
474 }
475 
476 
477 template <typename number>
478 inline const FullMatrix<number> &
480 {
481  if (same_diagonal())
482  return var_inverse_full[0];
483 
485  return var_inverse_full[i];
486 }
487 
488 
489 template <typename number>
490 inline const Householder<number> &
492 {
493  if (same_diagonal())
494  return var_inverse_householder[0];
495 
497  return var_inverse_householder[i];
498 }
499 
500 
501 template <typename number>
502 inline const LAPACKFullMatrix<number> &
504 {
505  if (same_diagonal())
506  return var_inverse_svd[0];
507 
509  return var_inverse_svd[i];
510 }
511 
512 
513 template <typename number>
514 inline const FullMatrix<number> &
516 {
518 
519  if (same_diagonal())
520  return var_diagonal[0];
521 
522  AssertIndexRange(i, var_diagonal.size());
523  return var_diagonal[i];
524 }
525 
526 
527 template <typename number>
528 inline FullMatrix<number> &
530 {
532 
533  if (same_diagonal())
534  return var_inverse_full[0];
535 
537  return var_inverse_full[i];
538 }
539 
540 
541 template <typename number>
542 inline Householder<number> &
544 {
546 
547  if (same_diagonal())
548  return var_inverse_householder[0];
549 
551  return var_inverse_householder[i];
552 }
553 
554 
555 template <typename number>
558 {
560 
561  if (same_diagonal())
562  return var_inverse_svd[0];
563 
565  return var_inverse_svd[i];
566 }
567 
568 
569 template <typename number>
570 inline FullMatrix<number> &
572 {
574 
575  if (same_diagonal())
576  return var_diagonal[0];
577 
578  AssertIndexRange(i, var_diagonal.size());
579  return var_diagonal[i];
580 }
581 
582 
583 template <typename number>
584 inline bool
586 {
587  return var_same_diagonal;
588 }
589 
590 
591 template <typename number>
592 inline bool
594 {
595  return var_store_diagonals;
596 }
597 
598 
599 template <typename number>
600 inline void
602 {
603  var_inverses_ready = x;
604 }
605 
606 
607 template <typename number>
608 inline bool
610 {
611  return var_inverses_ready;
612 }
613 
614 
615 template <typename number>
616 inline void
618 {
619  deallog << "PreconditionBlockBase: " << size() << " blocks; ";
620 
621  if (inversion == svd)
622  {
623  unsigned int kermin = 100000000, kermax = 0;
624  double sigmin = 1.e300, sigmax = -1.e300;
625  double kappamin = 1.e300, kappamax = -1.e300;
626 
627  for (size_type b = 0; b < size(); ++b)
628  {
630  size_type k = 1;
631  while (k <= matrix.n_cols() &&
632  matrix.singular_value(matrix.n_cols() - k) == 0)
633  ++k;
634  const double s0 = matrix.singular_value(0);
635  const double sm = matrix.singular_value(matrix.n_cols() - k);
636  const double co = sm / s0;
637 
638  if (kermin > k)
639  kermin = k - 1;
640  if (kermax < k)
641  kermax = k - 1;
642  if (s0 < sigmin)
643  sigmin = s0;
644  if (sm > sigmax)
645  sigmax = sm;
646  if (co < kappamin)
647  kappamin = co;
648  if (co > kappamax)
649  kappamax = co;
650  }
651  deallog << "dim ker [" << kermin << ':' << kermax << "] sigma [" << sigmin
652  << ':' << sigmax << "] kappa [" << kappamin << ':' << kappamax
653  << ']' << std::endl;
654  }
655  else if (inversion == householder)
656  {}
657  else if (inversion == gauss_jordan)
658  {}
659  else
660  {
661  Assert(false, ExcNotImplemented());
662  }
663 }
664 
665 
666 template <typename number>
667 inline std::size_t
669 {
670  std::size_t mem = sizeof(*this);
671  for (size_type i = 0; i < var_inverse_full.size(); ++i)
673  for (size_type i = 0; i < var_diagonal.size(); ++i)
675  return mem;
676 }
677 
678 
680 
681 #endif
LogStream deallog
Definition: logstream.cc:37
Contents is actually a matrix.
size_type n_cols() const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
static const char U
~PreconditionBlockBase()=default
void inverses_computed(bool are_they)
void inverse_vmult(size_type i, Vector< number2 > &dst, const Vector< number2 > &src) const
std::vector< FullMatrix< number > > var_diagonal
FullMatrix< number > & inverse(size_type i)
std::string compress(const std::string &input)
Definition: utilities.cc:392
PreconditionBlockBase(bool store_diagonals=false, Inversion method=gauss_jordan)
#define Assert(cond, exc)
Definition: exceptions.h:1419
#define DeclException0(Exception0)
Definition: exceptions.h:473
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
std::vector< FullMatrix< number > > var_inverse_full
void reinit(unsigned int nblocks, size_type blocksize, bool compress, Inversion method=gauss_jordan)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
number singular_value(const size_type i) const
void inverse_Tvmult(size_type i, Vector< number2 > &dst, const Vector< number2 > &src) const
unsigned int global_dof_index
Definition: types.h:76
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
static ::ExceptionBase & ExcDiagonalsNotStored()
std::vector< LAPACKFullMatrix< number > > var_inverse_svd
LAPACKFullMatrix< number > & inverse_svd(size_type i)
Householder< number > & inverse_householder(size_type i)
static ::ExceptionBase & ExcNotImplemented()
std::size_t memory_consumption() const
std::vector< Householder< number > > var_inverse_householder
FullMatrix< number > & diagonal(size_type i)
static ::ExceptionBase & ExcInverseNotAvailable()
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)