From ee9c9cde498accf97f66ad357270ca4e288778cf Mon Sep 17 00:00:00 2001 From: David Wells Date: Sun, 22 Jan 2017 12:04:16 -0500 Subject: [PATCH] Move the matrix-free project copy operation up the call stack. The matrix-free projection solves the linear system with a specific vector type and then copies the result into the given one. This commit moves that copy operation into a separate function so that the actual matrix free function no longer depends on the type of the input vector. Put another way: we no longer have to instantiate the matrix free projections for all possible vector types. Before moving the copy operation: File name Memory, MB Time, s source/numerics/vector_tools_project.cc 2774 76 source/numerics/vector_tools_project_inst2.cc 2829 81 source/numerics/vector_tools_project_inst3.cc 2823 86 After: File name Memory, MB Time, s source/numerics/vector_tools_project.cc 1855 47 source/numerics/vector_tools_project_inst2.cc 1995 51 source/numerics/vector_tools_project_inst3.cc 2183 55 --- .../deal.II/numerics/vector_tools.templates.h | 166 +++++++++++------- 1 file changed, 98 insertions(+), 68 deletions(-) diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index 1d93df16a3..2548266fea 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -899,16 +899,17 @@ namespace VectorTools * MatrixFree implementation of project() for an arbitrary number of * components and arbitrary degree of the FiniteElement. */ - template - void project_matrix_free (const Mapping &mapping, - const DoFHandler &dof, - const ConstraintMatrix &constraints, - const Quadrature &quadrature, - const Function &function, - VectorType &vec_result, - const bool enforce_zero_boundary, - const Quadrature &q_boundary, - const bool project_to_boundary_first) + template + void project_matrix_free + (const Mapping &mapping, + const DoFHandler &dof, + const ConstraintMatrix &constraints, + const Quadrature &quadrature, + const Function::value_type> &function, + LinearAlgebra::distributed::Vector &work_result, + const bool enforce_zero_boundary, + const Quadrature &q_boundary, + const bool project_to_boundary_first) { Assert (project_to_boundary_first == false, ExcNotImplemented()); Assert (enforce_zero_boundary == false, ExcNotImplemented()); @@ -916,34 +917,30 @@ namespace VectorTools (void) project_to_boundary_first; (void) q_boundary; - typedef typename VectorType::value_type number; Assert (dof.get_fe().n_components() == function.n_components, ExcDimensionMismatch(dof.get_fe().n_components(), function.n_components)); - Assert (vec_result.size() == dof.n_dofs(), - ExcDimensionMismatch (vec_result.size(), dof.n_dofs())); Assert (dof.get_fe().degree == fe_degree, ExcDimensionMismatch(fe_degree, dof.get_fe().degree)); Assert (dof.get_fe().n_components() == components, ExcDimensionMismatch(components, dof.get_fe().n_components())); // set up mass matrix and right hand side - typename MatrixFree::AdditionalData additional_data; + typename MatrixFree::AdditionalData additional_data; additional_data.tasks_parallel_scheme = - MatrixFree::AdditionalData::partition_color; + MatrixFree::AdditionalData::partition_color; additional_data.mapping_update_flags = (update_values | update_JxW_values); - std_cxx11::shared_ptr > matrix_free( - new MatrixFree ()); + std_cxx11::shared_ptr > matrix_free( + new MatrixFree ()); matrix_free->reinit (mapping, dof, constraints, QGauss<1>(fe_degree+2), additional_data); - typedef MatrixFreeOperators::MassOperator MatrixType; + typedef MatrixFreeOperators::MassOperator MatrixType; MatrixType mass_matrix; mass_matrix.initialize(matrix_free); mass_matrix.compute_diagonal(); - typedef LinearAlgebra::distributed::Vector LocalVectorType; - LocalVectorType vec, rhs, inhomogeneities; - matrix_free->initialize_dof_vector(vec); + LinearAlgebra::distributed::Vector rhs, inhomogeneities; + matrix_free->initialize_dof_vector(work_result); matrix_free->initialize_dof_vector(rhs); matrix_free->initialize_dof_vector(inhomogeneities); constraints.distribute(inhomogeneities); @@ -958,17 +955,13 @@ namespace VectorTools // badly conditioned matrices. This behavior can be observed, e.g. for // FE_Q_Hierarchical for degree higher than three. ReductionControl control(5.*rhs.size(), 0., 1e-12, false, false); - SolverCG > cg(control); + SolverCG > cg(control); PreconditionJacobi preconditioner; preconditioner.initialize(mass_matrix, 1.); - cg.solve (mass_matrix, vec, rhs, preconditioner); - vec+=inhomogeneities; + cg.solve (mass_matrix, work_result, rhs, preconditioner); + work_result+=inhomogeneities; - constraints.distribute (vec); - IndexSet::ElementIterator it = locally_owned_dofs.begin(); - for (; it!=locally_owned_dofs.end(); ++it) - vec_result(*it) = vec(*it); - vec_result.compress(VectorOperation::insert); + constraints.distribute (work_result); } @@ -979,64 +972,65 @@ namespace VectorTools * FiniteElement and call project_matrix_free with the appropriate * template arguments. */ - template - void project_matrix_free_degree (const Mapping &mapping, - const DoFHandler &dof, - const ConstraintMatrix &constraints, - const Quadrature &quadrature, - const Function &function, - VectorType &vec_result, - const bool enforce_zero_boundary, - const Quadrature &q_boundary, - const bool project_to_boundary_first) + template + void project_matrix_free_degree + (const Mapping &mapping, + const DoFHandler &dof, + const ConstraintMatrix &constraints, + const Quadrature &quadrature, + const Function::value_type> &function, + LinearAlgebra::distributed::Vector &work_result, + const bool enforce_zero_boundary, + const Quadrature &q_boundary, + const bool project_to_boundary_first) { switch (dof.get_fe().degree) { case 1: - VectorTools::project_matrix_free - (mapping, dof, constraints, quadrature, function, vec_result, + project_matrix_free + (mapping, dof, constraints, quadrature, function, work_result, enforce_zero_boundary, q_boundary, project_to_boundary_first); break; case 2: - VectorTools::project_matrix_free - (mapping, dof, constraints, quadrature, function, vec_result, + project_matrix_free + (mapping, dof, constraints, quadrature, function, work_result, enforce_zero_boundary, q_boundary, project_to_boundary_first); break; case 3: - VectorTools::project_matrix_free - (mapping, dof, constraints, quadrature, function, vec_result, + project_matrix_free + (mapping, dof, constraints, quadrature, function, work_result, enforce_zero_boundary, q_boundary, project_to_boundary_first); break; case 4: project_matrix_free - (mapping, dof, constraints, quadrature, function, vec_result, + (mapping, dof, constraints, quadrature, function, work_result, enforce_zero_boundary, q_boundary, project_to_boundary_first); break; case 5: project_matrix_free - (mapping, dof, constraints, quadrature, function, vec_result, + (mapping, dof, constraints, quadrature, function, work_result, enforce_zero_boundary, q_boundary, project_to_boundary_first); break; case 6: project_matrix_free - (mapping, dof, constraints, quadrature, function, vec_result, + (mapping, dof, constraints, quadrature, function, work_result, enforce_zero_boundary, q_boundary, project_to_boundary_first); break; case 7: project_matrix_free - (mapping, dof, constraints, quadrature, function, vec_result, + (mapping, dof, constraints, quadrature, function, work_result, enforce_zero_boundary, q_boundary, project_to_boundary_first); break; case 8: project_matrix_free - (mapping, dof, constraints, quadrature, function, vec_result, + (mapping, dof, constraints, quadrature, function, work_result, enforce_zero_boundary, q_boundary, project_to_boundary_first); break; @@ -1049,40 +1043,41 @@ namespace VectorTools // Helper interface for the matrix-free implementation of project(). // Used to determine the number of components. - template - void project_matrix_free_component (const Mapping &mapping, - const DoFHandler &dof, - const ConstraintMatrix &constraints, - const Quadrature &quadrature, - const Function &function, - VectorType &vec_result, - const bool enforce_zero_boundary, - const Quadrature &q_boundary, - const bool project_to_boundary_first) + template + void project_matrix_free_component + (const Mapping &mapping, + const DoFHandler &dof, + const ConstraintMatrix &constraints, + const Quadrature &quadrature, + const Function::value_type> &function, + LinearAlgebra::distributed::Vector &work_result, + const bool enforce_zero_boundary, + const Quadrature &q_boundary, + const bool project_to_boundary_first) { switch (dof.get_fe().n_components()) { case 1: project_matrix_free_degree<1> - (mapping, dof, constraints, quadrature, function, vec_result, + (mapping, dof, constraints, quadrature, function, work_result, enforce_zero_boundary, q_boundary, project_to_boundary_first); break; case 2: project_matrix_free_degree<2> - (mapping, dof, constraints, quadrature, function, vec_result, + (mapping, dof, constraints, quadrature, function, work_result, enforce_zero_boundary, q_boundary, project_to_boundary_first); break; case 3: project_matrix_free_degree<3> - (mapping, dof, constraints, quadrature, function, vec_result, + (mapping, dof, constraints, quadrature, function, work_result, enforce_zero_boundary, q_boundary, project_to_boundary_first); break; case 4: project_matrix_free_degree<4> - (mapping, dof, constraints, quadrature, function, vec_result, + (mapping, dof, constraints, quadrature, function, work_result, enforce_zero_boundary, q_boundary, project_to_boundary_first); break; @@ -1093,6 +1088,41 @@ namespace VectorTools + /** + * Helper interface for the matrix-free implementation of project(): avoid + * instantiating the other helper functions for more than one VectorType + * by copying from a LinearAlgebra::distributed::Vector. + */ + template + void project_matrix_free_copy_vector + (const Mapping &mapping, + const DoFHandler &dof, + const ConstraintMatrix &constraints, + const Quadrature &quadrature, + const Function &function, + VectorType &vec_result, + const bool enforce_zero_boundary, + const Quadrature &q_boundary, + const bool project_to_boundary_first) + { + Assert(vec_result.size() == dof.n_dofs(), + ExcDimensionMismatch(vec_result.size(), dof.n_dofs())); + + LinearAlgebra::distributed::Vector work_result; + project_matrix_free_component (mapping, dof, constraints, quadrature, + function, work_result, + enforce_zero_boundary, q_boundary, + project_to_boundary_first); + + const IndexSet locally_owned_dofs = dof.locally_owned_dofs(); + IndexSet::ElementIterator it = locally_owned_dofs.begin(); + for (; it!=locally_owned_dofs.end(); ++it) + vec_result(*it) = work_result(*it); + vec_result.compress(VectorOperation::insert); + } + + + /** * Specialization of project() for the case dim==spacedim. * Check if we can use the MatrixFree implementation or need @@ -1121,10 +1151,10 @@ namespace VectorTools use_matrix_free = false; if (use_matrix_free) - project_matrix_free_component (mapping, dof, constraints, quadrature, - function, vec_result, - enforce_zero_boundary, q_boundary, - project_to_boundary_first); + project_matrix_free_copy_vector (mapping, dof, constraints, quadrature, + function, vec_result, + enforce_zero_boundary, q_boundary, + project_to_boundary_first); else { Assert((dynamic_cast* > (&(dof.get_triangulation()))==0), -- 2.39.5