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>
30 DEAL_II_NAMESPACE_OPEN
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;
96 if (permutation.size() == M.m())
97 invert_permuted_diagblocks(permutation, inverse_permutation);
104 template <
class MATRIX,
typename inverse_type>
107 const std::vector<size_type> &permutation,
108 const std::vector<size_type> &inverse_permutation,
111 set_permutation(permutation, inverse_permutation);
112 initialize(M, parameters);
115 template <
class MATRIX,
typename inverse_type>
117 const std::vector<size_type> &permutation,
118 const std::vector<size_type> &inverse_permutation)
124 Assert (this->inverses_ready()==0, ExcInverseMatricesAlreadyExist());
130 if (this->same_diagonal())
132 deallog <<
"PreconditionBlock uses only one diagonal block" << std::endl;
134 for (
size_type row_cell=0; row_cell<blocksize; ++row_cell)
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)
140 if (entry->column() < blocksize)
141 M_cell(row_cell, entry->column()) = entry->value();
145 if (this->store_diagonals())
146 this->diagonal(0) = M_cell;
147 switch (this->inversion)
150 this->inverse(0).invert(M_cell);
153 this->inverse_householder(0).initialize(M_cell);
156 this->inverse_svd(0) = M_cell;
157 this->inverse_svd(0).compute_inverse_svd(0.);
175 for (
unsigned int cell=0; cell<this->size(); ++cell)
177 const size_type cell_start = cell*blocksize;
178 for (
size_type row_cell=0; row_cell<blocksize; ++row_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();
200 if (this->store_diagonals())
201 this->diagonal(cell) = M_cell;
202 switch (this->inversion)
205 this->inverse(cell).invert(M_cell);
208 this->inverse_householder(cell).initialize(M_cell);
211 this->inverse_svd(cell) = M_cell;
212 this->inverse_svd(cell).compute_inverse_svd(0.);
220 this->inverses_computed(
true);
225 template <
class MATRIX,
typename inverse_type>
226 template <
typename number2>
231 const bool transpose_diagonal)
const
237 if (permutation.size() != 0)
238 Assert (permutation.size() == M.m() || permutation.size() == this->size(),
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());
242 const bool cell_permuted = (permutation.size() == this->size());
265 for (
unsigned int rawcell=0; rawcell < this->size(); ++rawcell)
267 const unsigned int cell = cell_permuted ? permutation[rawcell] : rawcell;
268 const size_type block_start = cell*this->blocksize;
269 const size_type permuted_block_start = permuted
270 ? permutation[block_start]
276 for (row = permuted_block_start, row_cell = 0;
277 row_cell < this->blocksize;
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
289 ? inverse_permutation[column]
291 b_cell_row -= entry->value() * prev(column);
293 if (!this->inverses_ready()
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;
306 if (this->inverses_ready())
308 if (transpose_diagonal)
309 this->inverse_Tvmult(cell, x_cell, b_cell);
311 this->inverse_vmult(cell, x_cell, b_cell);
320 for (row=permuted_block_start, row_cell=0;
321 row_cell<this->blocksize;
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
340 if (permutation.size() != 0)
341 Assert (permutation.size() == M.m() || permutation.size() == this->size(),
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());
345 const bool cell_permuted = (permutation.size() == this->size());
360 for (
unsigned int rawcell=this->size(); rawcell!=0 ;)
363 const unsigned int cell = cell_permuted ? permutation[rawcell] : rawcell;
364 const size_type block_start = cell*this->blocksize;
365 const size_type block_end = block_start + this->blocksize;
366 const size_type permuted_block_start = permuted
367 ? permutation[block_start]
369 for (row = permuted_block_start, row_cell = 0;
370 row_cell<this->blocksize;
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
381 ? inverse_permutation[column]
383 b_cell_row -= entry->value() * prev(column);
384 if (!this->inverses_ready()
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;
403 if (this->inverses_ready())
405 if (transpose_diagonal)
406 this->inverse_Tvmult(cell, x_cell, b_cell);
408 this->inverse_vmult(cell, x_cell, b_cell);
418 for (row=permuted_block_start, row_cell=0;
419 row_cell<this->blocksize;
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>
441 Assert (this->inverses_ready()==0, ExcInverseMatricesAlreadyExist());
445 if (this->same_diagonal())
447 deallog <<
"PreconditionBlock uses only one diagonal block" << std::endl;
448 for (
size_type row_cell=0; row_cell<blocksize; ++row_cell)
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)
454 if (entry->column() < blocksize)
455 M_cell(row_cell, entry->column()) = entry->value();
459 if (this->store_diagonals())
460 this->diagonal(0) = M_cell;
461 switch (this->inversion)
464 this->inverse(0).invert(M_cell);
467 this->inverse_householder(0).initialize(M_cell);
470 this->inverse_svd(0) = M_cell;
471 this->inverse_svd(0).compute_inverse_svd(1.e-12);
481 for (
unsigned int cell=0; cell<this->size(); ++cell)
483 const size_type cell_start = cell*blocksize;
484 for (
size_type row_cell=0; row_cell<blocksize; ++row_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();
502 if (this->store_diagonals())
503 this->diagonal(cell) = M_cell;
504 switch (this->inversion)
507 this->inverse(cell).invert(M_cell);
510 this->inverse_householder(cell).initialize(M_cell);
513 this->inverse_svd(cell) = M_cell;
514 this->inverse_svd(cell).compute_inverse_svd(1.e-12);
521 this->inverses_computed(
true);
526 template <
class MATRIX,
typename inverse_type>
528 const std::vector<size_type> &p,
529 const std::vector<size_type> &i)
533 if (this->inverses_ready())
538 permutation.resize(p.size());
539 inverse_permutation.resize(p.size());
540 for (
unsigned int k=0; k<p.size(); ++k)
542 permutation[k] = p[k];
543 inverse_permutation[k] = i[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;
585 if (!this->inverses_ready())
588 for (
unsigned int cell=0; cell < this->size(); ++cell)
590 for (row=cell*this->blocksize, row_cell=0;
591 row_cell<this->blocksize;
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;
603 row_cell<this->blocksize;
606 dst(row)+=x_cell(row_cell);
608 dst(row)=x_cell(row_cell);
610 begin_diag_block+=this->blocksize;
614 for (
unsigned int cell=0; cell < this->size(); ++cell)
616 for (row=cell*this->blocksize, row_cell=0;
617 row_cell<this->blocksize;
620 b_cell(row_cell)=src(row);
622 this->inverse_vmult(cell, x_cell, b_cell);
624 for (row=cell*this->blocksize, row_cell=0;
625 row_cell<this->blocksize;
628 dst(row)+=x_cell(row_cell);
630 dst(row)=x_cell(row_cell);
632 begin_diag_block+=this->blocksize;
634 dst *= this->relaxation;
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>
688 this->forward_step(*aux, dst, src,
false);
693 template <
class MATRIX,
typename inverse_type>
694 template <
typename number2>
703 this->backward_step(*aux, dst, src,
true);
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
760 ? this->permutation[block_start]
763 for (row = permuted_block_start, row_cell = 0;
764 row_cell < this->blocksize;
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
775 ? this->inverse_permutation[column]
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;
791 if (this->inverses_ready())
793 if (transpose_diagonal)
794 this->inverse_Tvmult(cell, x_cell, b_cell);
796 this->inverse_vmult(cell, x_cell, b_cell);
805 for (row=permuted_block_start, row_cell=0;
806 row_cell<this->blocksize;
808 dst(row)=this->relaxation*x_cell(row_cell);
810 block_start+=this->blocksize;
815 template <
class MATRIX,
typename inverse_type>
816 template <
typename number2>
820 const bool transpose_diagonal,
826 const bool permuted = (this->permutation.size() != 0);
842 size_type block_end=this->blocksize * this->size();
846 for (
unsigned int cell=this->size(); cell!=0 ;)
849 const size_type block_start = block_end - this->blocksize;
851 const size_type permuted_block_start = (this->permutation.size() != 0)
852 ? this->permutation[block_start]
854 for (row = permuted_block_start, row_cell = 0;
855 row_cell<this->blocksize;
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
866 ? this->inverse_permutation[column]
868 if (inverse_permuted_column >= block_end)
869 b_cell_row -= entry->value() * dst(column);
870 else if (!this->inverses_ready() && column >= block_start)
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;
887 if (this->inverses_ready())
889 if (transpose_diagonal)
890 this->inverse_Tvmult(cell, x_cell, b_cell);
892 this->inverse_vmult(cell, x_cell, b_cell);
902 for (row=permuted_block_start, row_cell=0;
903 row_cell<this->blocksize;
905 dst(row)=this->relaxation*x_cell(row_cell);
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>
938 backward(dst, src,
true,
false);
942 template <
class MATRIX,
typename inverse_type>
943 template <
typename number2>
948 backward(dst, src,
true,
true);
953 template <
class MATRIX,
typename inverse_type>
954 template <
typename number2>
959 this->forward_step(dst, dst, src,
false);
963 template <
class MATRIX,
typename inverse_type>
964 template <
typename number2>
969 this->backward_step(dst, dst, src,
true);
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);
997 const double scaling = (2.-this->relaxation)/this->relaxation;
1000 for (
unsigned int cell=0; cell < this->size(); ++cell)
1004 for (
size_type row_cell=0; row_cell<this->blocksize; ++row_cell)
1005 cell_src(row_cell)=help(row+row_cell);
1007 this->diagonal(cell).vmult(cell_dst, cell_src);
1009 for (
size_type row_cell=0; row_cell<this->blocksize; ++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);
1028 const double scaling = (2.-this->relaxation)/this->relaxation;
1031 for (
unsigned int cell=0; cell < this->size(); ++cell)
1035 for (
size_type row_cell=0; row_cell<this->blocksize; ++row_cell)
1036 cell_src(row_cell)=help(row+row_cell);
1038 this->diagonal(cell).Tvmult(cell_dst, cell_src);
1040 for (
size_type row_cell=0; row_cell<this->blocksize; ++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>
1054 this->forward_step(dst, dst, src,
false);
1055 this->backward_step(dst, dst, src,
false);
1059 template <
class MATRIX,
typename inverse_type>
1060 template <
typename number2>
1065 this->backward_step(dst, dst, src,
true);
1066 this->forward_step(dst, dst, src,
true);
1071 DEAL_II_NAMESPACE_CLOSE
void Tvmult(Vector< number2 > &, const Vector< number2 > &) const
void Tvmult(Vector< number2 > &, const Vector< number2 > &) 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
::ExceptionBase & ExcMessage(std::string arg1)
void backward(Vector< number2 > &, const Vector< number2 > &, const bool transpose_diagonal, const bool adding) const
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 step(Vector< number2 > &dst, const Vector< number2 > &rhs) const
void initialize(const MATRIX &A, const AdditionalData parameters)
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
void step(Vector< number2 > &dst, const Vector< number2 > &rhs) const
#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
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)
virtual void reinit(const size_type N, const bool fast=false)
void set_permutation(const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation)
PreconditionBlockBase< inverse_type >::Inversion inversion
void vmult(Vector< number2 > &, const Vector< number2 > &) const
void Tstep(Vector< number2 > &dst, const Vector< number2 > &rhs) const
PreconditionBlock(bool store_diagonals=false)
::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 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
void Tstep(Vector< number2 > &dst, const Vector< number2 > &rhs) const
void Tvmult_add(Vector< number2 > &, const Vector< number2 > &) const
types::global_dof_index size_type
size_type block_size() const