17 #ifndef __deal2__precondition_h 18 #define __deal2__precondition_h 22 #include <deal.II/base/config.h> 23 #include <deal.II/base/smartpointer.h> 24 #include <deal.II/base/utilities.h> 25 #include <deal.II/base/parallel.h> 26 #include <deal.II/base/template_constraints.h> 27 #include <deal.II/lac/tridiagonal_matrix.h> 28 #include <deal.II/lac/solver_cg.h> 29 #include <deal.II/lac/vector_memory.h> 35 template <
typename number>
class Vector;
41 template <
typename number>
class Vector;
98 template <
class MATRIX>
99 void initialize (
const MATRIX &matrix,
105 template<
class VECTOR>
106 void vmult (VECTOR &,
const VECTOR &)
const;
115 template<
class VECTOR>
116 void Tvmult (VECTOR &,
const VECTOR &)
const;
121 template<
class VECTOR>
122 void vmult_add (VECTOR &,
const VECTOR &)
const;
131 template<
class VECTOR>
132 void Tvmult_add (VECTOR &,
const VECTOR &)
const;
203 template <
class MATRIX>
204 void initialize (
const MATRIX &,
210 template<
class VECTOR>
211 void vmult (VECTOR &,
const VECTOR &)
const;
220 template<
class VECTOR>
221 void Tvmult (VECTOR &,
const VECTOR &)
const;
225 template<
class VECTOR>
226 void vmult_add (VECTOR &,
const VECTOR &)
const;
235 template<
class VECTOR>
236 void Tvmult_add (VECTOR &,
const VECTOR &)
const;
299 template<
class MATRIX = SparseMatrix<
double>,
class VECTOR = Vector<
double> >
307 typedef void (
MATRIX::* function_ptr)(VECTOR &,
const VECTOR &)
const;
319 const function_ptr method);
327 void vmult (VECTOR &dst,
328 const VECTOR &src)
const;
352 template<
class MATRIX = SparseMatrix<
double> >
383 void initialize (
const MATRIX &A,
431 template <
class MATRIX = SparseMatrix<
double> >
438 template<
class VECTOR>
439 void vmult (VECTOR &,
const VECTOR &)
const;
448 template<
class VECTOR>
449 void Tvmult (VECTOR &,
const VECTOR &)
const;
456 template<
class VECTOR>
457 void step (VECTOR &x,
const VECTOR &rhs)
const;
464 template<
class VECTOR>
465 void Tstep (VECTOR &x,
const VECTOR &rhs)
const;
515 template <
class MATRIX = SparseMatrix<
double> >
522 template<
class VECTOR>
523 void vmult (VECTOR &,
const VECTOR &)
const;
529 template<
class VECTOR>
530 void Tvmult (VECTOR &,
const VECTOR &)
const;
537 template<
class VECTOR>
538 void step (VECTOR &x,
const VECTOR &rhs)
const;
545 template<
class VECTOR>
546 void Tstep (VECTOR &x,
const VECTOR &rhs)
const;
576 template <
class MATRIX = SparseMatrix<
double> >
601 void initialize (
const MATRIX &A,
607 template<
class VECTOR>
608 void vmult (VECTOR &,
const VECTOR &)
const;
617 template<
class VECTOR>
618 void Tvmult (VECTOR &,
const VECTOR &)
const;
626 template<
class VECTOR>
627 void step (VECTOR &x,
const VECTOR &rhs)
const;
634 template<
class VECTOR>
635 void Tstep (VECTOR &x,
const VECTOR &rhs)
const;
678 template <
class MATRIX = SparseMatrix<
double> >
706 void initialize (
const MATRIX &A,
707 const std::vector<size_type> &permutation,
708 const std::vector<size_type> &inverse_permutation,
715 template<
class VECTOR>
716 void vmult (VECTOR &,
const VECTOR &)
const;
722 template<
class VECTOR>
723 void Tvmult (VECTOR &,
const VECTOR &)
const;
789 template<
class SOLVER,
class MATRIX = SparseMatrix<
double>,
class PRECONDITION = PreconditionIdentity>
805 void initialize (
SOLVER &,
807 const PRECONDITION &);
812 template<
class VECTOR>
813 void vmult (VECTOR &,
const VECTOR &)
const;
847 template<
class MATRIX,
class PRECOND,
class VECTOR>
865 void vmult (VECTOR &dst,
const VECTOR &src)
const;
871 void Tvmult (VECTOR &dst,
const VECTOR &src)
const;
876 double residual (VECTOR &dst,
const VECTOR &src,
const VECTOR &rhs)
const;
909 template <
class MATRIX=SparseMatrix<
double>,
class VECTOR=Vector<
double> >
928 const double smoothing_range = 0.,
929 const bool nonzero_starting =
false,
930 const unsigned int eig_cg_n_iterations = 8,
931 const double eig_cg_residual = 1e-2,
932 const double max_eigenvalue = 1);
1005 void initialize (
const MATRIX &matrix,
1012 void vmult (VECTOR &dst,
1013 const VECTOR &src)
const;
1019 void Tvmult (VECTOR &dst,
1020 const VECTOR &src)
const;
1073 template <
class MATRIX>
1081 template<
class VECTOR>
1090 template<
class VECTOR>
1097 template<
class VECTOR>
1106 template<
class VECTOR>
1117 const double relaxation)
1119 relaxation(relaxation)
1129 relaxation=add_data.relaxation;
1143 template <
class MATRIX>
1154 template<
class VECTOR>
1158 dst.equ(relaxation,src);
1163 template<
class VECTOR>
1167 dst.equ(relaxation,src);
1170 template<
class VECTOR>
1174 dst.add(relaxation,src);
1179 template<
class VECTOR>
1183 dst.add(relaxation,src);
1188 template <
class MATRIX>
1194 relaxation = parameters.relaxation;
1198 template <
class MATRIX>
1208 template <
class MATRIX>
1209 template<
class VECTOR>
1214 this->A->precondition_Jacobi (dst, src, this->relaxation);
1219 template <
class MATRIX>
1220 template<
class VECTOR>
1225 this->A->precondition_Jacobi (dst, src, this->relaxation);
1230 template <
class MATRIX>
1231 template<
class VECTOR>
1236 this->A->Jacobi_step (dst, src, this->relaxation);
1241 template <
class MATRIX>
1242 template<
class VECTOR>
1253 template <
class MATRIX>
1254 template<
class VECTOR>
1259 this->A->precondition_SOR (dst, src, this->relaxation);
1264 template <
class MATRIX>
1265 template<
class VECTOR>
1270 this->A->precondition_TSOR (dst, src, this->relaxation);
1275 template <
class MATRIX>
1276 template<
class VECTOR>
1281 this->A->SOR_step (dst, src, this->relaxation);
1286 template <
class MATRIX>
1287 template<
class VECTOR>
1292 this->A->TSOR_step (dst, src, this->relaxation);
1299 template <
class MATRIX>
1302 const typename BaseClass::AdditionalData ¶meters)
1314 const size_type n = this->A->
n();
1315 pos_right_of_diagonal.resize(n, static_cast<std::size_t>(-1));
1316 for (size_type row=0; row<n; ++row)
1323 it = mat->
begin(row)+1;
1324 for ( ; it < mat->
end(row); ++it)
1325 if (it->column() > row)
1327 pos_right_of_diagonal[row] = it - mat->
begin();
1333 template <
class MATRIX>
1334 template<
class VECTOR>
1339 this->A->precondition_SSOR (dst, src, this->relaxation, pos_right_of_diagonal);
1344 template <
class MATRIX>
1345 template<
class VECTOR>
1350 this->A->precondition_SSOR (dst, src, this->relaxation, pos_right_of_diagonal);
1355 template <
class MATRIX>
1356 template<
class VECTOR>
1361 this->A->SSOR_step (dst, src, this->relaxation);
1366 template <
class MATRIX>
1367 template<
class VECTOR>
1378 template <
class MATRIX>
1382 const std::vector<size_type> &p,
1383 const std::vector<size_type> &ip,
1387 inverse_permutation = &ip;
1392 template <
class MATRIX>
1393 template<
class VECTOR>
1399 this->A->PSOR (dst, *permutation, *inverse_permutation, this->relaxation);
1404 template <
class MATRIX>
1405 template<
class VECTOR>
1411 this->A->TPSOR (dst, *permutation, *inverse_permutation, this->relaxation);
1418 template<
class MATRIX,
class VECTOR>
1420 const function_ptr method)
1422 matrix(M), precondition(method)
1427 template<
class MATRIX,
class VECTOR>
1430 const VECTOR &src)
const 1432 (matrix.*precondition)(dst, src);
1437 template<
class MATRIX>
1442 relaxation (relaxation)
1449 template<
class SOLVER,
class MATRIX,
class PRECONDITION>
1453 solver(0), matrix(0), precondition(0)
1457 template<
class SOLVER,
class MATRIX,
class PRECONDITION>
1462 const PRECONDITION &p)
1470 template<
class SOLVER,
class MATRIX,
class PRECONDITION>
1471 template<
class VECTOR>
1474 const VECTOR &src)
const 1476 Assert (solver !=0 && matrix != 0 && precondition != 0,
1479 solver->solve(*matrix, dst, src, *precondition);
1485 template<
class MATRIX,
class PRECOND,
class VECTOR>
1491 A(A), P(P), mem(mem)
1495 template<
class MATRIX,
class PRECOND,
class VECTOR>
1499 const VECTOR &src)
const 1501 VECTOR *h = mem.
alloc();
1510 template<
class MATRIX,
class PRECOND,
class VECTOR>
1514 const VECTOR &src)
const 1516 VECTOR *h = mem.
alloc();
1525 template<
class MATRIX,
class PRECOND,
class VECTOR>
1530 const VECTOR &rhs)
const 1532 VECTOR *h = mem.
alloc();
1537 dst.sadd(-1.,1.,rhs);
1538 return dst.l2_norm ();
1555 template <
typename VECTOR>
1558 vector_updates (
const VECTOR &src,
1559 const VECTOR &matrix_diagonal_inverse,
1560 const bool start_zero,
1561 const double factor1,
1562 const double factor2,
1569 dst.equ (factor2, src);
1570 dst.scale (matrix_diagonal_inverse);
1571 update1.equ(-1.,dst);
1576 update2.scale (matrix_diagonal_inverse);
1578 update1.equ(factor2, update2);
1580 update1.sadd(factor1, factor2, update2);
1586 template <
typename Number>
1591 VectorUpdatesRange (
const size_t size,
1593 const Number *matrix_diagonal_inverse,
1594 const bool start_zero,
1595 const Number factor1,
1596 const Number factor2,
1602 matrix_diagonal_inverse (matrix_diagonal_inverse),
1603 start_zero (start_zero),
1610 if (size < internal::Vector::minimum_parallel_grain_size)
1611 apply_to_subrange (0, size);
1613 apply_parallel (0, size,
1614 internal::Vector::minimum_parallel_grain_size);
1617 ~VectorUpdatesRange()
1621 apply_to_subrange (
const size_t begin,
1622 const size_t end)
const 1624 if (factor1 == Number())
1627 for (size_type i=begin; i<end; ++i)
1629 dst[i] = factor2 * src[i] * matrix_diagonal_inverse[i];
1630 update1[i] = -dst[i];
1633 for (size_type i=begin; i<end; ++i)
1635 update1[i] = ((update2[i]-src[i]) *
1636 factor2*matrix_diagonal_inverse[i]);
1637 dst[i] -= update1[i];
1641 for (size_type i=begin; i<end; ++i)
1643 const Number update2i = ((update2[i] - src[i]) *
1644 matrix_diagonal_inverse[i]);
1645 update1[i] = factor1 * update1[i] + factor2 * update2i;
1646 dst[i] -= update1[i];
1651 const Number *matrix_diagonal_inverse;
1652 const bool start_zero;
1653 const Number factor1;
1654 const Number factor2;
1655 mutable Number *update1;
1656 mutable Number *update2;
1657 mutable Number *dst;
1661 template <
typename Number>
1664 vector_updates (const ::Vector<Number> &src,
1665 const ::Vector<Number> &matrix_diagonal_inverse,
1666 const bool start_zero,
1667 const double factor1,
1668 const double factor2,
1673 VectorUpdatesRange<Number>(src.size(), src.begin(),
1674 matrix_diagonal_inverse.begin(),
1675 start_zero, factor1, factor2,
1676 update1.begin(), update2.begin(), dst.begin());
1680 template <
typename Number>
1685 const bool start_zero,
1686 const double factor1,
1687 const double factor2,
1693 matrix_diagonal_inverse.
begin(),
1694 start_zero, factor1, factor2,
1698 template <
typename VECTOR>
1699 struct DiagonalPreconditioner
1701 DiagonalPreconditioner (
const VECTOR &vector)
1703 diagonal_vector(vector)
1706 void vmult (VECTOR &dst,
1707 const VECTOR &src)
const 1710 dst.scale(diagonal_vector);
1713 const VECTOR &diagonal_vector;
1720 template <
class MATRIX,
class VECTOR>
1724 const double smoothing_range,
1725 const bool nonzero_starting,
1726 const unsigned int eig_cg_n_iterations,
1727 const double eig_cg_residual,
1728 const double max_eigenvalue)
1731 smoothing_range (smoothing_range),
1732 nonzero_starting (nonzero_starting),
1733 eig_cg_n_iterations (eig_cg_n_iterations),
1734 eig_cg_residual (eig_cg_residual),
1735 max_eigenvalue (max_eigenvalue)
1740 template <
class MATRIX,
class VECTOR>
1744 is_initialized (
false)
1749 template <
class MATRIX,
class VECTOR>
1755 matrix_ptr = &matrix;
1756 data = additional_data;
1757 if (data.matrix_diagonal_inverse.size() != matrix.m())
1759 Assert(data.matrix_diagonal_inverse.size() == 0,
1760 ExcMessage(
"Matrix diagonal vector set but not sized correctly"));
1761 data.matrix_diagonal_inverse.reinit(matrix.m());
1762 for (
unsigned int i=0; i<matrix.m(); ++i)
1763 data.matrix_diagonal_inverse(i) = 1./matrix.el(i,i);
1770 double max_eigenvalue, min_eigenvalue;
1771 if (data.eig_cg_n_iterations > 0)
1774 ExcMessage (
"Need to set at least two iterations to find eigenvalues."));
1778 static_cast<std::ostream *
>(0);
1782 std::ostringstream log_msg;
1788 VECTOR *rhs = memory.
alloc();
1789 VECTOR *dummy = memory.
alloc();
1790 rhs->reinit(data.matrix_diagonal_inverse,
true);
1791 dummy->reinit(data.matrix_diagonal_inverse);
1792 *rhs = 1./std::sqrt(static_cast<double>(matrix.m()));
1797 internal::PreconditionChebyshev::DiagonalPreconditioner<VECTOR>
1798 preconditioner(data.matrix_diagonal_inverse);
1801 solver.
solve(matrix, *dummy, *rhs, preconditioner);
1813 std::string cg_message = log_msg.str();
1814 const std::size_t pos = cg_message.find(
"cg:: ");
1815 if (pos != std::string::npos)
1817 cg_message.erase(0, pos+5);
1818 std::string first = cg_message;
1820 if (cg_message.find_first_of(
" ") != std::string::npos)
1821 first.erase(cg_message.find_first_of(
" "), std::string::npos);
1822 std::istringstream(first) >> min_eigenvalue;
1824 if (cg_message.find_last_of(
" ") != std::string::npos)
1826 cg_message.erase(0, cg_message.find_last_of(
" ")+1);
1827 std::istringstream(cg_message) >> max_eigenvalue;
1829 else max_eigenvalue = min_eigenvalue;
1832 min_eigenvalue = max_eigenvalue = 1;
1837 deallog.
attach(*old_stream,
false);
1841 max_eigenvalue *= 1.2;
1845 max_eigenvalue = data.max_eigenvalue;
1846 min_eigenvalue = data.max_eigenvalue/data.smoothing_range;
1849 const double alpha = (data.smoothing_range > 1. ?
1850 max_eigenvalue / data.smoothing_range :
1851 std::min(0.9*max_eigenvalue, min_eigenvalue));
1852 delta = (max_eigenvalue-alpha)*0.5;
1853 theta = (max_eigenvalue+alpha)*0.5;
1855 update1.reinit (data.matrix_diagonal_inverse,
true);
1856 update2.reinit (data.matrix_diagonal_inverse,
true);
1858 is_initialized =
true;
1863 template <
class MATRIX,
class VECTOR>
1867 const VECTOR &src)
const 1870 double rhok = delta / theta, sigma = theta / delta;
1871 if (data.nonzero_starting && !dst.all_zero())
1873 matrix_ptr->vmult (update2, dst);
1874 internal::PreconditionChebyshev::vector_updates
1875 (src, data.matrix_diagonal_inverse,
false, 0., 1./theta, update1,
1879 internal::PreconditionChebyshev::vector_updates
1880 (src, data.matrix_diagonal_inverse,
true, 0., 1./theta, update1,
1883 for (
unsigned int k=0; k<data.degree; ++k)
1885 matrix_ptr->vmult (update2, dst);
1886 const double rhokp = 1./(2.*sigma-rhok);
1887 const double factor1 = rhokp * rhok, factor2 = 2.*rhokp/delta;
1889 internal::PreconditionChebyshev::vector_updates
1890 (src, data.matrix_diagonal_inverse,
false, factor1, factor2, update1,
1897 template <
class MATRIX,
class VECTOR>
1901 const VECTOR &src)
const 1904 double rhok = delta / theta, sigma = theta / delta;
1905 if (data.nonzero_starting && !dst.all_zero())
1907 matrix_ptr->Tvmult (update2, dst);
1908 internal::PreconditionChebyshev::vector_updates
1909 (src, data.matrix_diagonal_inverse,
false, 0., 1./theta, update1,
1913 internal::PreconditionChebyshev::vector_updates
1914 (src, data.matrix_diagonal_inverse,
true, 0., 1./theta, update1,
1917 for (
unsigned int k=0; k<data.degree; ++k)
1919 matrix_ptr->Tvmult (update2, dst);
1920 const double rhokp = 1./(2.*sigma-rhok);
1921 const double factor1 = rhokp * rhok, factor2 = 2.*rhokp/delta;
1923 internal::PreconditionChebyshev::vector_updates
1924 (src, data.matrix_diagonal_inverse,
false, factor1, factor2, update1,
1931 template <
class MATRIX,
class VECTOR>
1935 is_initialized =
false;
1937 data.matrix_diagonal_inverse.reinit(0);
1946 DEAL_II_NAMESPACE_CLOSE
void initialize(const MATRIX &matrix, const AdditionalData &additional_data=AdditionalData())
void Tvmult(VECTOR &u, const VECTOR &v) const
const std::vector< size_type > * inverse_permutation
void Tvmult(VECTOR &, const VECTOR &) const
void initialize(const MATRIX &A, const AdditionalData ¶meters=AdditionalData())
void Tvmult_add(VECTOR &, const VECTOR &) const
void Tstep(VECTOR &x, const VECTOR &rhs) const
void Tvmult(VECTOR &, const VECTOR &) const
void vmult(VECTOR &u, const VECTOR &v) const
const function_ptr precondition
void vmult(VECTOR &dst, const VECTOR &src) const
void vmult(VECTOR &dst, const VECTOR &src) const
void Tvmult(VECTOR &dst, const VECTOR &src) const
void step(VECTOR &x, const VECTOR &rhs) const
SmartPointer< const MATRIX, PreconditionChebyshev< MATRIX, VECTOR > > matrix_ptr
::ExceptionBase & ExcMessage(std::string arg1)
void Tvmult(VECTOR &dst, const VECTOR &src) const
types::global_dof_index size_type
void Tvmult(VECTOR &, const VECTOR &) const
void vmult(VECTOR &dst, const VECTOR &src) const
virtual VECTOR * alloc()=0
void Tvmult(VECTOR &, const VECTOR &) const
double residual(VECTOR &dst, const VECTOR &src, const VECTOR &rhs) const
AdditionalData(const double relaxation=1.)
VectorMemory< VECTOR > & mem
void vmult_add(VECTOR &, const VECTOR &) const
std::ostream & get_file_stream()
const_iterator begin() const
VECTOR matrix_diagonal_inverse
SmartPointer< SOLVER, PreconditionLACSolver< SOLVER, MATRIX, PRECONDITION > > solver
void initialize(const AdditionalData ¶meters)
void vmult(VECTOR &, const VECTOR &) const
unsigned int global_dof_index
void Tvmult(VECTOR &, const VECTOR &) const
#define Assert(cond, exc)
void Tvmult(VECTOR &, const VECTOR &) const
void Tvmult_add(VECTOR &, const VECTOR &) const
PreconditionedMatrix(const MATRIX &A, const PRECOND &P, VectorMemory< VECTOR > &mem)
AdditionalData(const unsigned int degree=0, const double smoothing_range=0., const bool nonzero_starting=false, const unsigned int eig_cg_n_iterations=8, const double eig_cg_residual=1e-2, const double max_eigenvalue=1)
void step(VECTOR &x, const VECTOR &rhs) const
void vmult(VECTOR &, const VECTOR &) const
SmartPointer< const PRECONDITION, PreconditionLACSolver< SOLVER, MATRIX, PRECONDITION > > precondition
void initialize(const MATRIX &A, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const typename PreconditionRelaxation< MATRIX >::AdditionalData ¶meters=typename PreconditionRelaxation< MATRIX >::AdditionalData())
types::global_dof_index size_type
SmartPointer< const MATRIX, PreconditionRelaxation< MATRIX > > A
void initialize(SOLVER &, const MATRIX &, const PRECONDITION &)
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
const_iterator end() const
unsigned int last_step() const
void attach(std::ostream &o, const bool print_job_id=true)
types::global_dof_index size_type
virtual void free(const VECTOR *const)=0
void step(VECTOR &x, const VECTOR &rhs) const
void vmult(VECTOR &, const VECTOR &) const
unsigned int eig_cg_n_iterations
std::vector< std::size_t > pos_right_of_diagonal
void solve(const MATRIX &A, VECTOR &x, const VECTOR &b, const PRECONDITIONER &precondition)
const std::vector< size_type > * permutation
PreconditionUseMatrix(const MATRIX &M, const function_ptr method)
void initialize(const MATRIX &A, const typename BaseClass::AdditionalData ¶meters=typename BaseClass::AdditionalData())
void Tstep(VECTOR &x, const VECTOR &rhs) const
void vmult_add(VECTOR &, const VECTOR &) const
::ExceptionBase & ExcNotInitialized()
void vmult(VECTOR &, const VECTOR &) const
AdditionalData(const double relaxation=1.)
void vmult(VECTOR &, const VECTOR &) const
void Tstep(VECTOR &x, const VECTOR &rhs) const
SmartPointer< const MATRIX, PreconditionLACSolver< SOLVER, MATRIX, PRECONDITION > > matrix
virtual void free(const VECTOR *const)
size_type local_size() const
void vmult(VECTOR &, const VECTOR &) const
void vmult(VECTOR &, const VECTOR &) const
void initialize(const MATRIX &matrix, const AdditionalData &additional_data=AdditionalData())
PreconditionRelaxation< MATRIX > BaseClass