From 89226aa618ac7fee7980f4d4600efbfd1605892c Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Wed, 2 Nov 2016 08:34:37 +0100 Subject: [PATCH] use DiagonalMatrix in MatrixFreeOperators::Base To that end add a few missing function to DiagonalMatrix and also switch to Jacobi preconditioner in VectorTools::project() from quadrature points. --- include/deal.II/lac/diagonal_matrix.h | 28 ++++++++++ include/deal.II/matrix_free/operators.h | 51 +++++++++---------- .../deal.II/numerics/vector_tools.templates.h | 10 ++-- tests/matrix_free/mass_operator_02.cc | 4 +- 4 files changed, 58 insertions(+), 35 deletions(-) diff --git a/include/deal.II/lac/diagonal_matrix.h b/include/deal.II/lac/diagonal_matrix.h index 3238b340e4..ff9e27b8fd 100644 --- a/include/deal.II/lac/diagonal_matrix.h +++ b/include/deal.II/lac/diagonal_matrix.h @@ -53,6 +53,11 @@ public: */ VectorType &get_vector(); + /** + * Clear content of this object and reset to the state of default constructor. + */ + void clear(); + /** * Returns a read-only reference to the underlying vector. */ @@ -143,6 +148,11 @@ public: void Tvmult_add (VectorType &dst, const VectorType &src) const; + /** + * Return the memory consumption of this object. + */ + std::size_t memory_consumption () const; + private: /** * The stored vector. @@ -162,6 +172,24 @@ DiagonalMatrix::DiagonalMatrix() +template +void +DiagonalMatrix::clear() +{ + diagonal.reinit(0); +} + + + +template +std::size_t +DiagonalMatrix::memory_consumption () const +{ + return diagonal.memory_consumption(); +} + + + template void DiagonalMatrix::reinit(const VectorType &vec) diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h index 4e5a78c132..3b17805338 100644 --- a/include/deal.II/matrix_free/operators.h +++ b/include/deal.II/matrix_free/operators.h @@ -22,6 +22,7 @@ #include #include +#include #include #include #include @@ -38,7 +39,7 @@ namespace MatrixFreeOperators * the finest mesh or at a certain level in geometric multigrid. * * A derived class has to implement apply_add() method as well as - * compute_diagonal() to fill the protected member inverse_diagonal_entries. + * compute_diagonal() to initialize the protected member inverse_diagonal_entries. * In case of a non-symmetric operator, Tapply_add() should be additionally * implemented. * @@ -160,7 +161,7 @@ namespace MatrixFreeOperators /** * Get read access to the inverse diagonal of this operator. */ - const LinearAlgebra::distributed::Vector &get_matrix_diagonal_inverse() const; + const DiagonalMatrix > &get_matrix_diagonal_inverse() const; /** * Apply the Jacobi preconditioner, which multiplies every element of the @@ -201,7 +202,7 @@ namespace MatrixFreeOperators /** * A vector to store inverse of diagonal elements. */ - LinearAlgebra::distributed::Vector inverse_diagonal_entries; + DiagonalMatrix > inverse_diagonal_entries; private: @@ -611,7 +612,7 @@ namespace MatrixFreeOperators Base::clear () { data = NULL; - inverse_diagonal_entries.reinit(0); + inverse_diagonal_entries.clear(); } @@ -623,8 +624,8 @@ namespace MatrixFreeOperators { (void) col; Assert (row == col, ExcNotImplemented()); - Assert (inverse_diagonal_entries.size() > 0, ExcNotInitialized()); - return 1.0/inverse_diagonal_entries(row); + Assert (inverse_diagonal_entries.m() > 0, ExcNotInitialized()); + return 1.0/inverse_diagonal_entries(row,row); } @@ -866,10 +867,10 @@ namespace MatrixFreeOperators template - const LinearAlgebra::distributed::Vector & + const DiagonalMatrix > & Base::get_matrix_diagonal_inverse() const { - Assert(inverse_diagonal_entries.size() > 0, ExcNotInitialized()); + Assert(inverse_diagonal_entries.m() > 0, ExcNotInitialized()); return inverse_diagonal_entries; } @@ -891,10 +892,8 @@ namespace MatrixFreeOperators const LinearAlgebra::distributed::Vector &src, const Number omega) const { - Assert(inverse_diagonal_entries.size() > 0, ExcNotInitialized()); - - dst = src; - dst.scale(inverse_diagonal_entries); + Assert(inverse_diagonal_entries.m() > 0, ExcNotInitialized()); + inverse_diagonal_entries.vmult(dst,src); dst*= omega; } @@ -919,24 +918,22 @@ namespace MatrixFreeOperators Assert((Base::data != NULL), ExcNotInitialized()); LinearAlgebra::distributed::Vector ones; - Base::initialize_dof_vector(Base::inverse_diagonal_entries); - Base::initialize_dof_vector(ones); + LinearAlgebra::distributed::Vector &inverse_diagonal_entries = Base::inverse_diagonal_entries.get_vector(); + this->initialize_dof_vector(ones); + this->initialize_dof_vector(inverse_diagonal_entries); ones = 1.; ones.update_ghost_values(); - apply_add(Base::inverse_diagonal_entries, ones); + apply_add(inverse_diagonal_entries, ones); - const std::vector &constrained_dofs - = Base::data->get_constrained_dofs(); - for (unsigned int i=0; i< constrained_dofs.size(); ++i) - Base::inverse_diagonal_entries.local_element(constrained_dofs[i]) = 1.; + this->set_constrained_entries_to_one(inverse_diagonal_entries); - const unsigned int local_size = Base::inverse_diagonal_entries.local_size(); + const unsigned int local_size = inverse_diagonal_entries.local_size(); for (unsigned int i=0; i::inverse_diagonal_entries.local_element(i) - =1./Base::inverse_diagonal_entries.local_element(i); + inverse_diagonal_entries.local_element(i) + =1./inverse_diagonal_entries.local_element(i); - Base::inverse_diagonal_entries.compress(VectorOperation::insert); - Base::inverse_diagonal_entries.update_ghost_values(); + inverse_diagonal_entries.compress(VectorOperation::insert); + inverse_diagonal_entries.update_ghost_values(); } @@ -1028,7 +1025,7 @@ namespace MatrixFreeOperators Assert((Base::data != NULL), ExcNotInitialized()); unsigned int dummy = 0; - LinearAlgebra::distributed::Vector &inverse_diagonal_entries = Base::inverse_diagonal_entries; + LinearAlgebra::distributed::Vector &inverse_diagonal_entries = Base::inverse_diagonal_entries.get_vector(); this->initialize_dof_vector(inverse_diagonal_entries); Base:: data->cell_loop (&LaplaceOperator::local_diagonal_cell, @@ -1042,8 +1039,8 @@ namespace MatrixFreeOperators else inverse_diagonal_entries.local_element(i) = 1.; - Base::inverse_diagonal_entries.compress(VectorOperation::insert); - Base::inverse_diagonal_entries.update_ghost_values(); + inverse_diagonal_entries.compress(VectorOperation::insert); + inverse_diagonal_entries.update_ghost_values(); } diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index 453e4aa351..1decd4bbfd 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -970,9 +970,8 @@ namespace VectorTools //now invert the matrix ReductionControl control(rhs.size(), 0., 1e-12, false, false); SolverCG > cg(control); - typename PreconditionChebyshev::AdditionalData data; - data.matrix_diagonal_inverse = mass_matrix.get_matrix_diagonal_inverse(); - PreconditionChebyshev preconditioner; + typename PreconditionJacobi::AdditionalData data(0.8); + PreconditionJacobi preconditioner; preconditioner.initialize(mass_matrix, data); cg.solve (mass_matrix, vec, rhs, preconditioner); vec+=inhomogeneities; @@ -1041,9 +1040,8 @@ namespace VectorTools //now invert the matrix ReductionControl control(rhs.size(), 0., 1e-12, false, false); SolverCG > cg(control); - typename PreconditionChebyshev::AdditionalData data; - data.matrix_diagonal_inverse = mass_matrix.get_matrix_diagonal_inverse(); - PreconditionChebyshev preconditioner; + typename PreconditionJacobi::AdditionalData data(0.8); + PreconditionJacobi preconditioner; preconditioner.initialize(mass_matrix, data); cg.solve (mass_matrix, vec, rhs, preconditioner); vec+=inhomogeneities; diff --git a/tests/matrix_free/mass_operator_02.cc b/tests/matrix_free/mass_operator_02.cc index 1b83bb0c1b..4bcc628f4e 100644 --- a/tests/matrix_free/mass_operator_02.cc +++ b/tests/matrix_free/mass_operator_02.cc @@ -83,8 +83,8 @@ void test () MatrixFreeOperators::MassOperator mf; mf.initialize(mf_data); mf.compute_diagonal(); - const LinearAlgebra::distributed::Vector diagonal - = mf.get_matrix_diagonal_inverse(); + const LinearAlgebra::distributed::Vector &diagonal + = mf.get_matrix_diagonal_inverse().get_vector(); LinearAlgebra::distributed::Vector in, out, ref; mf_data.initialize_dof_vector (in); -- 2.39.5