Reference documentation for deal.II version 8.1.0
precondition_block_base.h
1 // ---------------------------------------------------------------------
2 // @f$Id: precondition_block_base.h 30036 2013-07-18 16:55:32Z maier @f$
3 //
4 // Copyright (C) 1999 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__precondition_block_base_h
18 #define __deal2__precondition_block_base_h
19 
20 
21 #include <deal.II/base/config.h>
23 #include <deal.II/base/subscriptor.h>
24 #include <deal.II/base/smartpointer.h>
25 #include <deal.II/base/memory_consumption.h>
26 #include <deal.II/lac/householder.h>
27 #include <deal.II/lac/lapack_full_matrix.h>
28 
29 #include <vector>
30 
32 
33 template <typename number> class FullMatrix;
34 template <typename number> class Vector;
35 
55 template <typename number>
57 {
58 public:
63 
69  enum Inversion
70  {
87  };
88 
94  Inversion method = gauss_jordan);
95 
100 
108  void clear();
109 
118  void reinit(unsigned int nblocks, size_type blocksize, bool compress,
119  Inversion method = gauss_jordan);
120 
125  void inverses_computed(bool are_they);
126 
140  void set_same_diagonal ();
141 
146  bool same_diagonal () const;
147 
153  bool store_diagonals() const;
154 
159  bool inverses_ready () const;
160 
164  bool empty () const;
165 
169  unsigned int size() const;
170 
177  number el(size_type i, size_type j) const;
178 
183  template <typename number2>
184  void inverse_vmult(size_type i, Vector<number2> &dst, const Vector<number2> &src) const;
185 
190  template <typename number2>
191  void inverse_Tvmult(size_type i, Vector<number2> &dst, const Vector<number2> &src) const;
192 
197  FullMatrix<number> &inverse (size_type i);
198 
204 
209  LAPACKFullMatrix<number> &inverse_svd (size_type i);
210 
215  const FullMatrix<number> &inverse (size_type i) const;
216 
221  const Householder<number> &inverse_householder (size_type i) const;
222 
227  const LAPACKFullMatrix<number> &inverse_svd (size_type i) const;
228 
233  FullMatrix<number> &diagonal (size_type i);
234 
239  const FullMatrix<number> &diagonal (size_type i) const;
240 
250  void log_statistics () const;
251 
257  std::size_t memory_consumption () const;
258 
265  DeclException0 (ExcDiagonalsNotStored);
266 
274  DeclException0 (ExcInverseNotAvailable);
275 
276 protected:
281 
282 private:
288  unsigned int n_diagonal_blocks;
289 
302  std::vector<FullMatrix<number> > var_inverse_full;
303 
316  std::vector<Householder<number> > var_inverse_householder;
317 
330  std::vector<LAPACKFullMatrix<number> > var_inverse_svd;
331 
337  std::vector<FullMatrix<number> > var_diagonal;
338 
339 
345 
351 
358 };
359 
360 //----------------------------------------------------------------------//
361 
362 template <typename number>
363 inline
365  bool store, Inversion method)
366  :
367  inversion(method),
369  var_store_diagonals(store),
370  var_same_diagonal(false),
371  var_inverses_ready(false)
372 {}
373 
374 
375 template <typename number>
376 inline
377 void
379 {
380  if (var_inverse_full.size()!=0)
381  var_inverse_full.erase(var_inverse_full.begin(), var_inverse_full.end());
382  if (var_inverse_householder.size()!=0)
384  if (var_inverse_svd.size()!=0)
385  var_inverse_svd.erase(var_inverse_svd.begin(), var_inverse_svd.end());
386  if (var_diagonal.size()!=0)
387  var_diagonal.erase(var_diagonal.begin(), var_diagonal.end());
388  var_same_diagonal = false;
389  var_inverses_ready = false;
390  n_diagonal_blocks = 0;
391 }
392 
393 template <typename number>
394 inline
395 void
396 PreconditionBlockBase<number>::reinit(unsigned int n, size_type b, bool compress,
397  Inversion method)
398 {
399  inversion = method;
400  var_same_diagonal = compress;
401  var_inverses_ready = false;
402  n_diagonal_blocks = n;
403 
404  if (compress)
405  {
406  switch (inversion)
407  {
408  case gauss_jordan:
409  var_inverse_full.resize(1);
410  var_inverse_full[0].reinit(b,b);
411  break;
412  case householder:
413  var_inverse_householder.resize(1);
414  break;
415  case svd:
416  var_inverse_svd.resize(1);
417  var_inverse_svd[0].reinit(b,b);
418  break;
419  default:
420  Assert(false, ExcNotImplemented());
421  }
422 
423  if (store_diagonals())
424  {
425  var_diagonal.resize(1);
426  var_diagonal[0].reinit(b,b);
427  }
428  }
429  else
430  {
431  // set the arrays to the right
432  // size. we could do it like this:
433  // var_inverse = vector<>(nblocks,FullMatrix<>())
434  // but this would involve copying many
435  // FullMatrix objects.
436  //
437  // the following is a neat trick which
438  // avoids copying
439  if (store_diagonals())
440  {
441  std::vector<FullMatrix<number> >
442  tmp(n, FullMatrix<number>(b));
443  var_diagonal.swap (tmp);
444  }
445 
446  switch (inversion)
447  {
448  case gauss_jordan:
449  {
450  std::vector<FullMatrix<number> >
451  tmp(n, FullMatrix<number>(b));
452  var_inverse_full.swap (tmp);
453  break;
454  }
455  case householder:
456  var_inverse_householder.resize(n);
457  break;
458  case svd:
459  {
460  std::vector<LAPACKFullMatrix<number> >
461  tmp(n, LAPACKFullMatrix<number>(b));
462  var_inverse_svd.swap (tmp);
463  break;
464  }
465  default:
466  Assert(false, ExcNotImplemented());
467  }
468  }
469 }
470 
471 
472 template <typename number>
473 inline
474 unsigned int
476 {
477  return n_diagonal_blocks;
478 }
479 
480 template <typename number>
481 template <typename number2>
482 inline
483 void
485  size_type i, Vector<number2> &dst, const Vector<number2> &src) const
486 {
487  const size_type ii = same_diagonal() ? 0U : i;
488 
489  switch (inversion)
490  {
491  case gauss_jordan:
492  AssertIndexRange (ii, var_inverse_full.size());
493  var_inverse_full[ii].vmult(dst, src);
494  break;
495  case householder:
497  var_inverse_householder[ii].vmult(dst, src);
498  break;
499  case svd:
500  AssertIndexRange (ii, var_inverse_svd.size());
501  var_inverse_svd[ii].vmult(dst, src);
502  break;
503  default:
504  Assert(false, ExcNotImplemented());
505  }
506 }
507 
508 
509 template <typename number>
510 template <typename number2>
511 inline
512 void
514  size_type i, Vector<number2> &dst, const Vector<number2> &src) const
515 {
516  const size_type ii = same_diagonal() ? 0U : i;
517 
518  switch (inversion)
519  {
520  case gauss_jordan:
521  AssertIndexRange (ii, var_inverse_full.size());
522  var_inverse_full[ii].Tvmult(dst, src);
523  break;
524  case householder:
526  var_inverse_householder[ii].Tvmult(dst, src);
527  break;
528  case svd:
529  AssertIndexRange (ii, var_inverse_svd.size());
530  var_inverse_svd[ii].Tvmult(dst, src);
531  break;
532  default:
533  Assert(false, ExcNotImplemented());
534  }
535 }
536 
537 
538 template <typename number>
539 inline
540 const FullMatrix<number> &
542 {
543  if (same_diagonal())
544  return var_inverse_full[0];
545 
546  Assert (i < var_inverse_full.size(), ExcIndexRange(i,0,var_inverse_full.size()));
547  return var_inverse_full[i];
548 }
549 
550 
551 template <typename number>
552 inline
553 const Householder<number> &
555 {
556  if (same_diagonal())
557  return var_inverse_householder[0];
558 
560  return var_inverse_householder[i];
561 }
562 
563 
564 template <typename number>
565 inline
568 {
569  if (same_diagonal())
570  return var_inverse_svd[0];
571 
572  AssertIndexRange (i, var_inverse_svd.size());
573  return var_inverse_svd[i];
574 }
575 
576 
577 template <typename number>
578 inline
579 const FullMatrix<number> &
581 {
582  Assert(store_diagonals(), ExcDiagonalsNotStored());
583 
584  if (same_diagonal())
585  return var_diagonal[0];
586 
587  Assert (i < var_diagonal.size(), ExcIndexRange(i,0,var_diagonal.size()));
588  return var_diagonal[i];
589 }
590 
591 
592 template <typename number>
593 inline
596 {
597  Assert(var_inverse_full.size() != 0, ExcInverseNotAvailable());
598 
599  if (same_diagonal())
600  return var_inverse_full[0];
601 
602  Assert (i < var_inverse_full.size(), ExcIndexRange(i,0,var_inverse_full.size()));
603  return var_inverse_full[i];
604 }
605 
606 
607 template <typename number>
608 inline
611 {
612  Assert(var_inverse_householder.size() != 0, ExcInverseNotAvailable());
613 
614  if (same_diagonal())
615  return var_inverse_householder[0];
616 
618  return var_inverse_householder[i];
619 }
620 
621 
622 template <typename number>
623 inline
626 {
627  Assert(var_inverse_svd.size() != 0, ExcInverseNotAvailable());
628 
629  if (same_diagonal())
630  return var_inverse_svd[0];
631 
632  AssertIndexRange (i, var_inverse_svd.size());
633  return var_inverse_svd[i];
634 }
635 
636 
637 template <typename number>
638 inline
641 {
642  Assert(store_diagonals(), ExcDiagonalsNotStored());
643 
644  if (same_diagonal())
645  return var_diagonal[0];
646 
647  Assert (i < var_diagonal.size(), ExcIndexRange(i,0,var_diagonal.size()));
648  return var_diagonal[i];
649 }
650 
651 
652 template <typename number>
653 inline bool
655 {
656  return var_same_diagonal;
657 }
658 
659 
660 template <typename number>
661 inline bool
663 {
664  return var_store_diagonals;
665 }
666 
667 
668 template <typename number>
669 inline void
671 {
672  var_inverses_ready = x;
673 }
674 
675 
676 template <typename number>
677 inline bool
679 {
680  return var_inverses_ready;
681 }
682 
683 
684 template <typename number>
685 inline void
687 {
688  deallog << "PreconditionBlockBase: " << size() << " blocks; ";
689 
690  if (inversion == svd)
691  {
692  unsigned int kermin = 100000000, kermax = 0;
693  double sigmin = 1.e300, sigmax= -1.e300;
694  double kappamin = 1.e300, kappamax= -1.e300;
695 
696  for (size_type b=0; b<size(); ++b)
697  {
698  const LAPACKFullMatrix<number> &matrix = inverse_svd(b);
699  size_type k=1;
700  while (k <= matrix.n_cols() && matrix.singular_value(matrix.n_cols()-k) == 0)
701  ++k;
702  const double s0 = matrix.singular_value(0);
703  const double sm = matrix.singular_value(matrix.n_cols()-k);
704  const double co = sm/s0;
705 
706  if (kermin > k) kermin = k-1;
707  if (kermax < k) kermax = k-1;
708  if (s0 < sigmin) sigmin = s0;
709  if (sm > sigmax) sigmax = sm;
710  if (co < kappamin) kappamin = co;
711  if (co > kappamax) kappamax = co;
712  }
713  deallog << "dim ker [" << kermin << ':' << kermax
714  << "] sigma [" << sigmin << ':' << sigmax
715  << "] kappa [" << kappamin << ':' << kappamax << ']' << std::endl;
716 
717  }
718  else if (inversion == householder)
719  {
720  }
721  else if (inversion == gauss_jordan)
722  {
723  }
724  else
725  {
726  Assert(false, ExcNotImplemented());
727  }
728 }
729 
730 
731 template <typename number>
732 inline
733 std::size_t
735 {
736  std::size_t mem = sizeof(*this);
737  for (size_type i=0; i<var_inverse_full.size(); ++i)
739  for (size_type i=0; i<var_diagonal.size(); ++i)
741  return mem;
742 }
743 
744 
745 DEAL_II_NAMESPACE_CLOSE
746 
747 #endif
void inverse_Tvmult(size_type i, Vector< number2 > &dst, const Vector< number2 > &src) const
number singular_value(const size_type i) const
types::global_dof_index size_type
DeclException0(ExcDiagonalsNotStored)
#define AssertIndexRange(index, range)
Definition: exceptions.h:888
unsigned int size() const
void inverses_computed(bool are_they)
std::vector< FullMatrix< number > > var_diagonal
FullMatrix< number > & inverse(size_type i)
PreconditionBlockBase(bool store_diagonals=false, Inversion method=gauss_jordan)
unsigned int global_dof_index
Definition: types.h:100
#define Assert(cond, exc)
Definition: exceptions.h:299
number el(size_type i, size_type j) const
std::size_t memory_consumption(const T &t)
bool empty() const
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
std::vector< FullMatrix< number > > var_inverse_full
void reinit(unsigned int nblocks, size_type blocksize, bool compress, Inversion method=gauss_jordan)
std::size_t memory_consumption() const
std::vector< LAPACKFullMatrix< number > > var_inverse_svd
LAPACKFullMatrix< number > & inverse_svd(size_type i)
Householder< number > & inverse_householder(size_type i)
::ExceptionBase & ExcNotImplemented()
std::vector< Householder< number > > var_inverse_householder
FullMatrix< number > & diagonal(size_type i)
unsigned int n_cols() const
void inverse_vmult(size_type i, Vector< number2 > &dst, const Vector< number2 > &src) const