17 #ifndef __deal2__precondition_block_templates_h 18 #define __deal2__precondition_block_templates_h 21 #include <deal.II/base/config.h> 23 #include <deal.II/base/logstream.h> 24 #include <deal.II/base/memory_consumption.h> 25 #include <deal.II/lac/householder.h> 26 #include <deal.II/lac/precondition_block.h> 27 #include <deal.II/lac/vector.h> 28 #include <deal.II/lac/full_matrix.h> 33 template<
class MATRIX,
typename inverse_type>
36 const double relaxation,
37 const bool invert_diagonal,
38 const bool same_diagonal)
40 relaxation (relaxation),
41 block_size(block_size),
42 invert_diagonal(invert_diagonal),
43 same_diagonal(same_diagonal),
49 template <
typename number>
54 template <
class MATRIX,
typename inverse_type>
58 A(0, typeid(*this).name())
62 template <
class MATRIX,
typename inverse_type>
67 template <
class MATRIX,
typename inverse_type>
76 template <
class MATRIX,
typename inverse_type>
84 Assert (M.m() == M.n(), ExcNotQuadratic());
87 Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m()));
90 const unsigned int nblocks = A->m()/bsize;
104 template <
class MATRIX,
typename inverse_type>
115 template <
class MATRIX,
typename inverse_type>
132 deallog <<
"PreconditionBlock uses only one diagonal block" << std::endl;
136 typename MATRIX::const_iterator entry = M.begin(row_cell);
137 const typename MATRIX::const_iterator row_end = M.end(row_cell);
138 while (entry != row_end)
141 M_cell(row_cell, entry->column()) = entry->value();
175 for (
unsigned int cell=0; cell<this->
size(); ++cell)
180 const size_type urow = row_cell + cell_start;
184 typename MATRIX::const_iterator entry = M.begin(row);
185 const typename MATRIX::const_iterator row_end = M.end(row);
187 for (; entry != row_end; ++entry)
190 if (inverse_permutation[entry->column()]<cell_start)
193 const size_type column_cell = inverse_permutation[entry->column()]-cell_start;
194 if (column_cell >= blocksize)
196 M_cell(row_cell, column_cell) = entry->value();
225 template <
class MATRIX,
typename inverse_type>
226 template <
typename number2>
231 const bool transpose_diagonal)
const 239 ExcMessage(
"Permutation vector size must be equal to either the number of blocks or the dimension of the system"));
241 const bool permuted = (
permutation.size() == M.m());
265 for (
unsigned int rawcell=0; rawcell < this->
size(); ++rawcell)
267 const unsigned int cell = cell_permuted ?
permutation[rawcell] : rawcell;
269 const size_type permuted_block_start = permuted
276 for (row = permuted_block_start, row_cell = 0;
281 const typename MATRIX::const_iterator row_end = M.end(row);
282 typename MATRIX::const_iterator entry = M.begin(row);
285 for (; entry != row_end; ++entry)
287 const size_type column = entry->column();
288 const size_type inverse_permuted_column = permuted
291 b_cell_row -= entry->value() * prev(column);
294 && inverse_permuted_column >= block_start
295 && inverse_permuted_column < block_start + this->blocksize)
297 const size_type column_cell = column - block_start;
298 if (transpose_diagonal)
299 M_cell(column_cell, row_cell) = entry->value();
301 M_cell(row_cell, column_cell) = entry->value();
304 b_cell(row_cell)=b_cell_row;
308 if (transpose_diagonal)
320 for (row=permuted_block_start, row_cell=0;
323 dst(row) = prev(row) + this->
relaxation*x_cell(row_cell);
328 template <
class MATRIX,
typename inverse_type>
329 template <
typename number2>
334 const bool transpose_diagonal)
const 342 ExcMessage(
"Permutation vector size must be equal to either the number of blocks or the dimension of the system"));
344 const bool permuted = (
permutation.size() == M.m());
360 for (
unsigned int rawcell=this->
size(); rawcell!=0 ;)
363 const unsigned int cell = cell_permuted ?
permutation[rawcell] : rawcell;
366 const size_type permuted_block_start = permuted
369 for (row = permuted_block_start, row_cell = 0;
373 const typename MATRIX::const_iterator row_end = M.end(row);
374 typename MATRIX::const_iterator entry = M.begin(row);
377 for (; entry != row_end; ++entry)
379 const size_type column = entry->column();
380 const size_type inverse_permuted_column = permuted
383 b_cell_row -= entry->value() * prev(column);
385 && inverse_permuted_column < block_end
386 && column >= block_start)
388 const size_type column_cell = column - block_start;
395 if (transpose_diagonal)
396 M_cell(column_cell, row_cell) = entry->value();
398 M_cell(row_cell, column_cell) = entry->value();
401 b_cell(row_cell)=b_cell_row;
405 if (transpose_diagonal)
418 for (row=permuted_block_start, row_cell=0;
421 dst(row) = prev(row) + this->
relaxation*x_cell(row_cell);
426 template <
class MATRIX,
typename inverse_type>
434 template <
class MATRIX,
typename inverse_type>
447 deallog <<
"PreconditionBlock uses only one diagonal block" << std::endl;
450 typename MATRIX::const_iterator entry = M.begin(row_cell);
451 const typename MATRIX::const_iterator row_end = M.end(row_cell);
452 while (entry != row_end)
455 M_cell(row_cell, entry->column()) = entry->value();
481 for (
unsigned int cell=0; cell<this->
size(); ++cell)
486 const size_type row = row_cell + cell_start;
487 typename MATRIX::const_iterator entry = M.begin(row);
488 const typename MATRIX::const_iterator row_end = M.end(row);
490 for (; entry != row_end; ++entry)
492 if (entry->column()<cell_start)
495 const size_type column_cell = entry->column()-cell_start;
496 if (column_cell >= blocksize)
498 M_cell(row_cell, column_cell) = entry->value();
526 template <
class MATRIX,
typename inverse_type>
528 const std::vector<size_type> &p,
529 const std::vector<size_type> &i)
540 for (
unsigned int k=0; k<p.size(); ++k)
548 template <
class MATRIX,
typename inverse_type>
552 return (
sizeof(*
this)
563 template <
class MATRIX,
typename inverse_type>
564 template <
typename number2>
583 size_type row, row_cell, begin_diag_block=0;
588 for (
unsigned int cell=0; cell < this->
size(); ++cell)
590 for (row=cell*this->
blocksize, row_cell=0;
594 b_cell(row_cell)=src(row);
595 for (
size_type column_cell=0, column=cell*this->blocksize;
596 column_cell<this->
blocksize; ++column_cell, ++column)
597 M_cell(row_cell,column_cell)=M(row,column);
602 for (row=cell*this->blocksize, row_cell=0;
606 dst(row)+=x_cell(row_cell);
608 dst(row)=x_cell(row_cell);
614 for (
unsigned int cell=0; cell < this->
size(); ++cell)
616 for (row=cell*this->
blocksize, row_cell=0;
620 b_cell(row_cell)=src(row);
624 for (row=cell*this->blocksize, row_cell=0;
628 dst(row)+=x_cell(row_cell);
630 dst(row)=x_cell(row_cell);
638 template <
class MATRIX,
typename inverse_type>
639 template <
typename number2>
644 do_vmult(dst, src,
false);
648 template <
class MATRIX,
typename inverse_type>
649 template <
typename number2>
654 do_vmult(dst, src,
false);
658 template <
class MATRIX,
typename inverse_type>
659 template <
typename number2>
664 do_vmult(dst, src,
true);
668 template <
class MATRIX,
typename inverse_type>
669 template <
typename number2>
674 do_vmult(dst, src,
true);
678 template <
class MATRIX,
typename inverse_type>
679 template <
typename number2>
693 template <
class MATRIX,
typename inverse_type>
694 template <
typename number2>
713 template <
class MATRIX,
typename inverse_type>
719 template <
class MATRIX,
typename inverse_type>
725 template <
class MATRIX,
typename inverse_type>
726 template <
typename number2>
730 const bool transpose_diagonal,
736 const bool permuted = (this->
permutation.size() != 0);
757 for (
unsigned int cell=0; cell < this->
size(); ++cell)
759 const size_type permuted_block_start = permuted
763 for (row = permuted_block_start, row_cell = 0;
767 const typename MATRIX::const_iterator row_end = M.end(row);
768 typename MATRIX::const_iterator entry = M.begin(row);
771 for (; entry != row_end; ++entry)
773 const size_type column = entry->column();
774 const size_type inverse_permuted_column = permuted
778 if (inverse_permuted_column < block_start)
779 b_cell_row -= entry->value() * dst(column);
780 else if (!this->
inverses_ready() && column < block_start + this->blocksize)
782 const size_type column_cell = column - block_start;
783 if (transpose_diagonal)
784 M_cell(column_cell, row_cell) = entry->value();
786 M_cell(row_cell, column_cell) = entry->value();
789 b_cell(row_cell)=b_cell_row;
793 if (transpose_diagonal)
805 for (row=permuted_block_start, row_cell=0;
815 template <
class MATRIX,
typename inverse_type>
816 template <
typename number2>
820 const bool transpose_diagonal,
826 const bool permuted = (this->
permutation.size() != 0);
846 for (
unsigned int cell=this->
size(); cell!=0 ;)
854 for (row = permuted_block_start, row_cell = 0;
858 const typename MATRIX::const_iterator row_end = M.end(row);
859 typename MATRIX::const_iterator entry = M.begin(row);
862 for (; entry != row_end; ++entry)
864 const size_type column = entry->column();
865 const size_type inverse_permuted_column = permuted
868 if (inverse_permuted_column >= block_end)
869 b_cell_row -= entry->value() * dst(column);
872 const size_type column_cell = column - block_start;
879 if (transpose_diagonal)
880 M_cell(column_cell, row_cell) = entry->value();
882 M_cell(row_cell, column_cell) = entry->value();
885 b_cell(row_cell)=b_cell_row;
889 if (transpose_diagonal)
902 for (row=permuted_block_start, row_cell=0;
906 block_end = block_start;
912 template <
class MATRIX,
typename inverse_type>
913 template <
typename number2>
918 forward(dst, src,
false,
false);
922 template <
class MATRIX,
typename inverse_type>
923 template <
typename number2>
928 forward(dst, src,
false,
true);
932 template <
class MATRIX,
typename inverse_type>
933 template <
typename number2>
942 template <
class MATRIX,
typename inverse_type>
943 template <
typename number2>
953 template <
class MATRIX,
typename inverse_type>
954 template <
typename number2>
963 template <
class MATRIX,
typename inverse_type>
964 template <
typename number2>
978 template <
class MATRIX,
typename inverse_type>
985 template <
class MATRIX,
typename inverse_type>
986 template <
typename number2>
993 this->
forward(help, src,
false,
false);
1000 for (
unsigned int cell=0; cell < this->
size(); ++cell)
1005 cell_src(row_cell)=help(row+row_cell);
1010 help(row+row_cell) = scaling * cell_dst(row_cell);
1013 this->
backward(dst, help,
false,
false);
1016 template <
class MATRIX,
typename inverse_type>
1017 template <
typename number2>
1024 this->
backward(help, src,
true,
false);
1031 for (
unsigned int cell=0; cell < this->
size(); ++cell)
1036 cell_src(row_cell)=help(row+row_cell);
1041 help(row+row_cell) = scaling * cell_dst(row_cell);
1044 this->
forward(dst, help,
true,
false);
1048 template <
class MATRIX,
typename inverse_type>
1049 template <
typename number2>
1059 template <
class MATRIX,
typename inverse_type>
1060 template <
typename number2>
1071 DEAL_II_NAMESPACE_CLOSE
void Tvmult(Vector< number2 > &, const Vector< number2 > &) const
void Tvmult(Vector< number2 > &, const Vector< number2 > &) const
void inverse_Tvmult(size_type i, Vector< number2 > &dst, const Vector< number2 > &src) const
#define AssertDimension(dim1, dim2)
void backward_step(Vector< number2 > &dst, const Vector< number2 > &prev, const Vector< number2 > &src, const bool transpose_diagonal) const
void step(Vector< number2 > &dst, const Vector< number2 > &rhs) const
bool store_diagonals() const
::ExceptionBase & ExcMessage(std::string arg1)
unsigned int size() const
void backward(Vector< number2 > &, const Vector< number2 > &, const bool transpose_diagonal, const bool adding) const
void invert(const FullMatrix< number2 > &M)
AdditionalData(const size_type block_size, const double relaxation=1., const bool invert_diagonal=true, const bool same_diagonal=false)
void vmult(Vector< number2 > &, const Vector< number2 > &) const
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
std::vector< size_type > inverse_permutation
void step(Vector< number2 > &dst, const Vector< number2 > &rhs) const
void initialize(const MATRIX &A, const AdditionalData parameters)
void inverses_computed(bool are_they)
void do_vmult(Vector< number2 > &, const Vector< number2 > &, bool adding) const
void Tstep(Vector< number2 > &dst, const Vector< number2 > &rhs) const
void forward_step(Vector< number2 > &dst, const Vector< number2 > &prev, const Vector< number2 > &src, const bool transpose_diagonal) const
std::size_t memory_consumption() const
FullMatrix< inverse_type > & inverse(size_type i)
void step(Vector< number2 > &dst, const Vector< number2 > &rhs) const
virtual void reinit(const size_type N, const bool fast=false)
#define Assert(cond, exc)
void forward(Vector< number2 > &, const Vector< number2 > &, const bool transpose_diagonal, const bool adding) const
void Tvmult_add(Vector< number2 > &, const Vector< number2 > &) const
bool same_diagonal() const
void invert_permuted_diagblocks(const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation)
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
void Tvmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
void reinit(unsigned int nblocks, size_type blocksize, bool compress, Inversion method=gauss_jordan)
void set_permutation(const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation)
void compute_inverse_svd(const double threshold=0.)
PreconditionBlockBase< inverse_type >::Inversion inversion
void vmult(Vector< number2 > &, const Vector< number2 > &) const
void Tstep(Vector< number2 > &dst, const Vector< number2 > &rhs) const
std::vector< size_type > permutation
PreconditionBlock(bool store_diagonals=false)
LAPACKFullMatrix< inverse_type > & inverse_svd(size_type i)
Householder< inverse_type > & inverse_householder(size_type i)
::ExceptionBase & ExcNotImplemented()
::ExceptionBase & ExcNotInitialized()
void Tvmult(Vector< number2 > &, const Vector< number2 > &) const
void vmult_add(Vector< number2 > &, const Vector< number2 > &) const
void vmult_add(Vector< number2 > &, const Vector< number2 > &) const
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void initialize(const FullMatrix< number2 > &)
void vmult(Vector< number2 > &, const Vector< number2 > &) const
SmartPointer< const MATRIX, PreconditionBlock< MATRIX, inverse_type > > A
double least_squares(Vector< number2 > &dst, const Vector< number2 > &src) const
bool inverses_ready() const
FullMatrix< inverse_type > & diagonal(size_type i)
void Tstep(Vector< number2 > &dst, const Vector< number2 > &rhs) const
void Tvmult_add(Vector< number2 > &, const Vector< number2 > &) const
void inverse_vmult(size_type i, Vector< number2 > &dst, const Vector< number2 > &src) const
types::global_dof_index size_type
size_type block_size() const