17 #ifndef __deal2__precondition_block_base_h
18 #define __deal2__precondition_block_base_h
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>
34 template <
typename number>
class Vector;
55 template <
typename number>
118 void reinit(
unsigned int nblocks, size_type blocksize,
bool compress,
169 unsigned int size()
const;
177 number
el(size_type i, size_type j)
const;
183 template <
typename number2>
190 template <
typename number2>
362 template <
typename number>
368 n_diagonal_blocks(0),
369 var_store_diagonals(store),
370 var_same_diagonal(false),
371 var_inverses_ready(false)
375 template <
typename number>
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)
383 var_inverse_householder.erase(var_inverse_householder.begin(), var_inverse_householder.end());
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;
393 template <
typename number>
400 var_same_diagonal = compress;
401 var_inverses_ready =
false;
402 n_diagonal_blocks = n;
409 var_inverse_full.resize(1);
410 var_inverse_full[0].reinit(b,b);
413 var_inverse_householder.resize(1);
416 var_inverse_svd.resize(1);
417 var_inverse_svd[0].reinit(b,b);
423 if (store_diagonals())
425 var_diagonal.resize(1);
426 var_diagonal[0].reinit(b,b);
439 if (store_diagonals())
441 std::vector<FullMatrix<number> >
443 var_diagonal.swap (tmp);
450 std::vector<FullMatrix<number> >
452 var_inverse_full.swap (tmp);
456 var_inverse_householder.resize(n);
460 std::vector<LAPACKFullMatrix<number> >
462 var_inverse_svd.swap (tmp);
472 template <
typename number>
477 return n_diagonal_blocks;
480 template <
typename number>
481 template <
typename number2>
487 const size_type ii = same_diagonal() ? 0U : i;
493 var_inverse_full[ii].vmult(dst, src);
497 var_inverse_householder[ii].vmult(dst, src);
501 var_inverse_svd[ii].vmult(dst, src);
509 template <
typename number>
510 template <
typename number2>
516 const size_type ii = same_diagonal() ? 0U : i;
522 var_inverse_full[ii].Tvmult(dst, src);
526 var_inverse_householder[ii].Tvmult(dst, src);
530 var_inverse_svd[ii].Tvmult(dst, src);
538 template <
typename number>
544 return var_inverse_full[0];
547 return var_inverse_full[i];
551 template <
typename number>
557 return var_inverse_householder[0];
560 return var_inverse_householder[i];
564 template <
typename number>
570 return var_inverse_svd[0];
573 return var_inverse_svd[i];
577 template <
typename number>
582 Assert(store_diagonals(), ExcDiagonalsNotStored());
585 return var_diagonal[0];
588 return var_diagonal[i];
592 template <
typename number>
597 Assert(var_inverse_full.size() != 0, ExcInverseNotAvailable());
600 return var_inverse_full[0];
603 return var_inverse_full[i];
607 template <
typename number>
612 Assert(var_inverse_householder.size() != 0, ExcInverseNotAvailable());
615 return var_inverse_householder[0];
618 return var_inverse_householder[i];
622 template <
typename number>
627 Assert(var_inverse_svd.size() != 0, ExcInverseNotAvailable());
630 return var_inverse_svd[0];
633 return var_inverse_svd[i];
637 template <
typename number>
642 Assert(store_diagonals(), ExcDiagonalsNotStored());
645 return var_diagonal[0];
648 return var_diagonal[i];
652 template <
typename number>
656 return var_same_diagonal;
660 template <
typename number>
664 return var_store_diagonals;
668 template <
typename number>
672 var_inverses_ready = x;
676 template <
typename number>
680 return var_inverses_ready;
684 template <
typename number>
688 deallog <<
"PreconditionBlockBase: " << size() <<
" blocks; ";
690 if (inversion == svd)
692 unsigned int kermin = 100000000, kermax = 0;
693 double sigmin = 1.e300, sigmax= -1.e300;
694 double kappamin = 1.e300, kappamax= -1.e300;
704 const double co = sm/s0;
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;
713 deallog <<
"dim ker [" << kermin <<
':' << kermax
714 <<
"] sigma [" << sigmin <<
':' << sigmax
715 <<
"] kappa [" << kappamin <<
':' << kappamax <<
']' << std::endl;
718 else if (inversion == householder)
721 else if (inversion == gauss_jordan)
731 template <
typename number>
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)
745 DEAL_II_NAMESPACE_CLOSE
void inverse_Tvmult(size_type i, Vector< number2 > &dst, const Vector< number2 > &src) const
bool store_diagonals() const
types::global_dof_index size_type
DeclException0(ExcDiagonalsNotStored)
#define AssertIndexRange(index, range)
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
#define Assert(cond, exc)
number el(size_type i, size_type j) const
std::size_t memory_consumption(const T &t)
bool same_diagonal() 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)
void log_statistics() const
unsigned int n_diagonal_blocks
std::size_t memory_consumption() const
number singular_value(const size_type i) 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
bool inverses_ready() const
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