Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
precondition_block_base.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2018 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 
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/memory_consumption.h>
24 #include <deal.II/base/smartpointer.h>
25 #include <deal.II/base/subscriptor.h>
26 
27 #include <deal.II/lac/householder.h>
28 #include <deal.II/lac/lapack_full_matrix.h>
29 
30 #include <vector>
31 
32 DEAL_II_NAMESPACE_OPEN
33 
34 template <typename number>
35 class FullMatrix;
36 template <typename number>
37 class Vector;
38 
57 template <typename number>
59 {
60 public:
65 
70  enum Inversion
71  {
85  };
86 
91  Inversion method = gauss_jordan);
92 
96  ~PreconditionBlockBase() = default;
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
131  store_diagonals() const;
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 
167  inverse(size_type i);
168 
174 
180 
184  const FullMatrix<number> &
185  inverse(size_type i) const;
186 
190  const Householder<number> &
192 
197  inverse_svd(size_type i) const;
198 
203  diagonal(size_type i);
204 
208  const FullMatrix<number> &
209  diagonal(size_type i) const;
210 
216  void
217  log_statistics() const;
218 
223  std::size_t
224  memory_consumption() const;
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  Assert(i < var_inverse_full.size(),
482  ExcIndexRange(i, 0, var_inverse_full.size()));
483  return var_inverse_full[i];
484 }
485 
486 
487 template <typename number>
488 inline const Householder<number> &
490 {
491  if (same_diagonal())
492  return var_inverse_householder[0];
493 
494  AssertIndexRange(i, var_inverse_householder.size());
495  return var_inverse_householder[i];
496 }
497 
498 
499 template <typename number>
500 inline const LAPACKFullMatrix<number> &
502 {
503  if (same_diagonal())
504  return var_inverse_svd[0];
505 
506  AssertIndexRange(i, var_inverse_svd.size());
507  return var_inverse_svd[i];
508 }
509 
510 
511 template <typename number>
512 inline const FullMatrix<number> &
514 {
515  Assert(store_diagonals(), ExcDiagonalsNotStored());
516 
517  if (same_diagonal())
518  return var_diagonal[0];
519 
520  Assert(i < var_diagonal.size(), ExcIndexRange(i, 0, var_diagonal.size()));
521  return var_diagonal[i];
522 }
523 
524 
525 template <typename number>
526 inline FullMatrix<number> &
528 {
529  Assert(var_inverse_full.size() != 0, ExcInverseNotAvailable());
530 
531  if (same_diagonal())
532  return var_inverse_full[0];
533 
534  Assert(i < var_inverse_full.size(),
535  ExcIndexRange(i, 0, var_inverse_full.size()));
536  return var_inverse_full[i];
537 }
538 
539 
540 template <typename number>
541 inline Householder<number> &
543 {
544  Assert(var_inverse_householder.size() != 0, ExcInverseNotAvailable());
545 
546  if (same_diagonal())
547  return var_inverse_householder[0];
548 
549  AssertIndexRange(i, var_inverse_householder.size());
550  return var_inverse_householder[i];
551 }
552 
553 
554 template <typename number>
557 {
558  Assert(var_inverse_svd.size() != 0, ExcInverseNotAvailable());
559 
560  if (same_diagonal())
561  return var_inverse_svd[0];
562 
563  AssertIndexRange(i, var_inverse_svd.size());
564  return var_inverse_svd[i];
565 }
566 
567 
568 template <typename number>
569 inline FullMatrix<number> &
571 {
572  Assert(store_diagonals(), ExcDiagonalsNotStored());
573 
574  if (same_diagonal())
575  return var_diagonal[0];
576 
577  Assert(i < var_diagonal.size(), ExcIndexRange(i, 0, var_diagonal.size()));
578  return var_diagonal[i];
579 }
580 
581 
582 template <typename number>
583 inline bool
585 {
586  return var_same_diagonal;
587 }
588 
589 
590 template <typename number>
591 inline bool
593 {
594  return var_store_diagonals;
595 }
596 
597 
598 template <typename number>
599 inline void
601 {
602  var_inverses_ready = x;
603 }
604 
605 
606 template <typename number>
607 inline bool
609 {
610  return var_inverses_ready;
611 }
612 
613 
614 template <typename number>
615 inline void
617 {
618  deallog << "PreconditionBlockBase: " << size() << " blocks; ";
619 
620  if (inversion == svd)
621  {
622  unsigned int kermin = 100000000, kermax = 0;
623  double sigmin = 1.e300, sigmax = -1.e300;
624  double kappamin = 1.e300, kappamax = -1.e300;
625 
626  for (size_type b = 0; b < size(); ++b)
627  {
628  const LAPACKFullMatrix<number> &matrix = inverse_svd(b);
629  size_type k = 1;
630  while (k <= matrix.n_cols() &&
631  matrix.singular_value(matrix.n_cols() - k) == 0)
632  ++k;
633  const double s0 = matrix.singular_value(0);
634  const double sm = matrix.singular_value(matrix.n_cols() - k);
635  const double co = sm / s0;
636 
637  if (kermin > k)
638  kermin = k - 1;
639  if (kermax < k)
640  kermax = k - 1;
641  if (s0 < sigmin)
642  sigmin = s0;
643  if (sm > sigmax)
644  sigmax = sm;
645  if (co < kappamin)
646  kappamin = co;
647  if (co > kappamax)
648  kappamax = co;
649  }
650  deallog << "dim ker [" << kermin << ':' << kermax << "] sigma [" << sigmin
651  << ':' << sigmax << "] kappa [" << kappamin << ':' << kappamax
652  << ']' << std::endl;
653  }
654  else if (inversion == householder)
655  {}
656  else if (inversion == gauss_jordan)
657  {}
658  else
659  {
660  Assert(false, ExcNotImplemented());
661  }
662 }
663 
664 
665 template <typename number>
666 inline 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 
678 DEAL_II_NAMESPACE_CLOSE
679 
680 #endif
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
~PreconditionBlockBase()=default
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
LinearAlgebra::distributed::Vector< Number > Vector
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)
PreconditionBlockBase(bool store_diagonals=false, Inversion method=gauss_jordan)
#define Assert(cond, exc)
Definition: exceptions.h:1407
#define DeclException0(Exception0)
Definition: exceptions.h:473
std::vector< FullMatrix< number > > var_inverse_full
void reinit(unsigned int nblocks, size_type blocksize, bool compress, Inversion method=gauss_jordan)
void inverse_Tvmult(size_type i, Vector< number2 > &dst, const Vector< number2 > &src) const
unsigned int global_dof_index
Definition: types.h:89
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)