From bf33ccd56bb4fd1cd5a4422c0436c2ce3a1f7455 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 9 Apr 2001 09:38:34 +0000 Subject: [PATCH] Make the functions in MatrixCreator use multiple threads if supported. git-svn-id: https://svn.dealii.org/trunk@4401 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/numerics/matrices.h | 248 ++++++++++++---- deal.II/deal.II/source/numerics/matrices.cc | 301 +++++++++++++++++--- 2 files changed, 459 insertions(+), 90 deletions(-) diff --git a/deal.II/deal.II/include/numerics/matrices.h b/deal.II/deal.II/include/numerics/matrices.h index d0fdc993f7..9000db266a 100644 --- a/deal.II/deal.II/include/numerics/matrices.h +++ b/deal.II/deal.II/include/numerics/matrices.h @@ -11,6 +11,7 @@ #include +#include #include @@ -205,7 +206,11 @@ class MatrixCreator * Assemble the mass matrix. If no * coefficient is given, it is assumed * to be unity. - * + * + * If the library is configured + * to use multithreading, this + * function works in parallel. + * * See the general doc of this class * for more information. */ @@ -213,7 +218,7 @@ class MatrixCreator const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function *a = 0); + const Function * const a = 0); /** * Calls the @p{create_mass_matrix} @@ -223,16 +228,20 @@ class MatrixCreator static void create_mass_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function *a = 0); + const Function * const a = 0); /** - * Assemble the mass matrix and a right - * hand side vector. If no - * coefficient is given, it is assumed - * to be unity. + * Assemble the mass matrix and a + * right hand side vector. If no + * coefficient is given, it is + * assumed to be unity. * - * See the general doc of this class - * for more information. + * If the library is configured + * to use multithreading, this + * function works in parallel. + * + * See the general doc of this + * class for more information. */ static void create_mass_matrix (const Mapping &mapping, const DoFHandler &dof, @@ -240,7 +249,7 @@ class MatrixCreator SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function *a = 0); + const Function * const a = 0); /** * Calls the @p{create_mass_matrix} @@ -252,52 +261,65 @@ class MatrixCreator SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function *a = 0); + const Function * const a = 0); /** - * Assemble the mass matrix and a right - * hand side vector along the boundary. - * If no - * coefficient is given, it is assumed - * to be constant one. + * Assemble the mass matrix and a + * right hand side vector along + * the boundary. If no + * coefficient is given, it is + * assumed to be constant one. * - * The matrix is assumed to already be - * initialized with a suiting sparsity - * pattern (the @ref{DoFHandler} provides an + * The matrix is assumed to + * already be initialized with a + * suiting sparsity pattern (the + * @ref{DoFHandler} provides an * appropriate function). * - * See the general doc of this class - * for more information. + * If the library is configured + * to use multithreading, this + * function works in parallel. + * + * See the general doc of this + * class for more information. */ - static void create_boundary_mass_matrix (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const FunctionMap &boundary_functions, - Vector &rhs_vector, - std::vector&dof_to_boundary_mapping, - const Function *a = 0); + static + void create_boundary_mass_matrix (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const FunctionMap &boundary_functions, + Vector &rhs_vector, + std::vector&dof_to_boundary_mapping, + const Function *a = 0); /** - * Calls the @p{create_boundary_mass_matrix} + * Calls the + * @p{create_boundary_mass_matrix} * function, see above, with * @p{mapping=MappingQ1()}. */ - static void create_boundary_mass_matrix (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const FunctionMap &boundary_functions, - Vector &rhs_vector, - std::vector&dof_to_boundary_mapping, - const Function *a = 0); + static + void create_boundary_mass_matrix (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const FunctionMap &boundary_functions, + Vector &rhs_vector, + std::vector&dof_to_boundary_mapping, + const Function *a = 0); /** - * Assemble the Laplace matrix. If no - * coefficient is given, it is assumed - * to be constant one. + * Assemble the Laplace + * matrix. If no coefficient is + * given, it is assumed to be + * constant one. * - * See the general doc of this class - * for more information. + * If the library is configured + * to use multithreading, this + * function works in parallel. + * + * See the general doc of this + * class for more information. */ static void create_laplace_matrix (const Mapping &mapping, const DoFHandler &dof, @@ -316,10 +338,11 @@ class MatrixCreator const Function *a = 0); /** - * Generate Laplace matrix for a given level. + * Generate Laplace matrix for a + * given level. * - * See the general doc of this class - * for more information. + * See the general doc of this + * class for more information. */ static void create_level_laplace_matrix (unsigned int level, const MGDoFHandler& dof, @@ -328,13 +351,18 @@ class MatrixCreator const Function* a = 0); /** - * Assemble the Laplace matrix and a right - * hand side vector. If no - * coefficient is given, it is assumed - * to be constant one. + * Assemble the Laplace matrix + * and a right hand side + * vector. If no coefficient is + * given, it is assumed to be + * constant one. * - * See the general doc of this class - * for more information. + * If the library is configured + * to use multithreading, this + * function works in parallel. + * + * See the general doc of this + * class for more information. */ static void create_laplace_matrix (const Mapping &mapping, const DoFHandler &dof, @@ -342,7 +370,7 @@ class MatrixCreator SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function *a = 0); + const Function * const a = 0); /** * Calls the @p{create_laplace_matrix} @@ -360,6 +388,122 @@ class MatrixCreator * Exception */ DeclException0 (ExcComponentMismatch); + + private: + /** + * Convenience abbreviation for + * DoF handler cell iterators. + */ + typedef typename DoFHandler::active_cell_iterator active_cell_iterator; + + /** + * Pair of iterators denoting a + * half-open range. + */ + typedef std::pair IteratorRange; + + + /** + * Version of the same function + * (without suffix @p{_1}) with + * the same argument list that + * operates only on an interval + * of iterators. Used for + * parallelization. The mutex is + * used to synchronise access to + * the matrix. + */ + static + void create_mass_matrix_1 (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function * const a, + const IteratorRange &range, + Threads::ThreadMutex &mutex); + + /** + * Version of the same function + * (without suffix @p{_2}) with + * the same argument list that + * operates only on an interval + * of iterators. Used for + * parallelization. The mutex is + * used to synchronise access to + * the matrix. + */ + static + void create_mass_matrix_2 (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const a, + const IteratorRange &range, + Threads::ThreadMutex &mutex); + + /** + * Version of the same function + * (without suffix @p{_1}) with + * the same argument list that + * operates only on an interval + * of iterators. Used for + * parallelization. The mutex is + * used to synchronise access to + * the matrix. + */ + static + void create_laplace_matrix_1 (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function * const a, + const IteratorRange &range, + Threads::ThreadMutex &mutex); + + /** + * Version of the same function + * (without suffix @p{_2}) with + * the same argument list that + * operates only on an interval + * of iterators. Used for + * parallelization. The mutex is + * used to synchronise access to + * the matrix. + */ + static + void create_laplace_matrix_2 (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const a, + const IteratorRange &range, + Threads::ThreadMutex &mutex); + + /** + * Version of the same function + * (without suffix @p{_1}) with + * the same argument list that + * operates only on an interval + * of iterators. Used for + * parallelization. The mutex is + * used to synchronise access to + * the matrix. + */ + static + void create_boundary_mass_matrix_1 (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const FunctionMap &boundary_functions, + Vector &rhs_vector, + std::vector&dof_to_boundary_mapping, + const Function * const a, + const IteratorRange &range, + Threads::ThreadMutex &mutex); }; diff --git a/deal.II/deal.II/source/numerics/matrices.cc b/deal.II/deal.II/source/numerics/matrices.cc index 64bd7d2a6b..88334a0133 100644 --- a/deal.II/deal.II/source/numerics/matrices.cc +++ b/deal.II/deal.II/source/numerics/matrices.cc @@ -13,6 +13,7 @@ #include +#include #include #include #include @@ -40,14 +41,46 @@ using namespace std; -// TODO:[RH, WB] extend this function to use vector valued coefficient functions for system elements. -// TODO:[WB] implement multithreading for this function. template void MatrixCreator::create_mass_matrix (const Mapping &mapping, const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, const Function * const coefficient) +{ + const unsigned int n_threads = multithread_info.n_default_threads; + Threads::ThreadManager thread_manager; + + // define starting and end point + // for each thread + std::vector thread_ranges + = Threads::split_range (dof.begin_active(), + dof.end(), n_threads); + + // mutex to synchronise access to + // the matrix + Threads::ThreadMutex mutex; + + // then assemble in parallel + for (unsigned int thread=0; thread::create_mass_matrix_1) + .collect_args (mapping, dof, q, matrix, coefficient, + thread_ranges[thread], mutex)); + thread_manager.wait (); +}; + + + +// TODO:[RH] extend this function to use vector valued coefficient functions for system elements. +template +void MatrixCreator::create_mass_matrix_1 (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function * const coefficient, + const IteratorRange &range, + Threads::ThreadMutex &mutex) { UpdateFlags update_flags = UpdateFlags(update_values | update_JxW_values); if (coefficient != 0) @@ -65,8 +98,8 @@ void MatrixCreator::create_mass_matrix (const Mapping &mapping, std::vector dof_indices (dofs_per_cell); - typename DoFHandler::active_cell_iterator cell = dof.begin_active(); - for (; cell!=dof.end(); ++cell) + typename DoFHandler::active_cell_iterator cell = range.first; + for (; cell!=range.second; ++cell) { fe_values.reinit (cell); @@ -104,10 +137,15 @@ void MatrixCreator::create_mass_matrix (const Mapping &mapping, // transfer everything into the // global object + mutex.acquire (); for (unsigned int i=0; i::create_mass_matrix (const DoFHandler &dof, } -// TODO:[RH, WB] extend this function to use vector valued coefficient functions for system elements. -// TODO:[WB] implement multithreading for this function. template void MatrixCreator::create_mass_matrix (const Mapping &mapping, const DoFHandler &dof, @@ -133,6 +169,45 @@ void MatrixCreator::create_mass_matrix (const Mapping &mapping, const Function &rhs, Vector &rhs_vector, const Function * const coefficient) +{ + const unsigned int n_threads = multithread_info.n_default_threads; + Threads::ThreadManager thread_manager; + + // define starting and end point + // for each thread + std::vector thread_ranges + = Threads::split_range (dof.begin_active(), + dof.end(), n_threads); + + // mutex to synchronise access to + // the matrix + Threads::ThreadMutex mutex; + + // then assemble in parallel + for (unsigned int thread=0; thread:: + create_mass_matrix_2) + .collect_args (mapping, dof, q, matrix, rhs, + rhs_vector, coefficient, + thread_ranges[thread], mutex)); + thread_manager.wait (); +}; + + + +// TODO:[RH] extend this function to use vector valued coefficient functions for system elements. +template +void +MatrixCreator::create_mass_matrix_2 (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, + const IteratorRange &range, + Threads::ThreadMutex &mutex) { UpdateFlags update_flags = UpdateFlags(update_values | update_q_points | @@ -154,8 +229,8 @@ void MatrixCreator::create_mass_matrix (const Mapping &mapping, std::vector dof_indices (dofs_per_cell); - typename DoFHandler::active_cell_iterator cell = dof.begin_active(); - for (; cell!=dof.end(); ++cell) + typename DoFHandler::active_cell_iterator cell = range.first; + for (; cell!=range.second; ++cell) { fe_values.reinit (cell); @@ -205,12 +280,17 @@ void MatrixCreator::create_mass_matrix (const Mapping &mapping, // transfer everything into the // global object + mutex.acquire (); for (unsigned int i=0; i::create_boundary_mass_matrix (const Mapping<1> &, #endif -// TODO:[RH, WB] extend this function to use vector valued coefficient functions for system elements. -// TODO:[WB] implement multithreading for this function. +// TODO:[RH] extend this function to use vector valued coefficient functions for system elements. template -void MatrixCreator::create_boundary_mass_matrix (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const FunctionMap &boundary_functions, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function *a) +void +MatrixCreator::create_boundary_mass_matrix (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const FunctionMap &boundary_functions, + Vector &rhs_vector, + std::vector &dof_to_boundary_mapping, + const Function *a) +{ + const unsigned int n_threads = multithread_info.n_default_threads; + Threads::ThreadManager thread_manager; + + // define starting and end point + // for each thread + std::vector thread_ranges + = Threads::split_range (dof.begin_active(), + dof.end(), n_threads); + + // mutex to synchronise access to + // the matrix + Threads::ThreadMutex mutex; + + // then assemble in parallel + for (unsigned int thread=0; thread:: + create_boundary_mass_matrix_1) + .collect_args (mapping, dof, q, matrix, + boundary_functions, rhs_vector, + dof_to_boundary_mapping, a, + thread_ranges[thread], mutex)); + thread_manager.wait (); +}; + + +#if deal_II_dimension == 1 + +template <> +void +MatrixCreator<1>:: +create_boundary_mass_matrix_1 (const Mapping<1> &, + const DoFHandler<1> &, + const Quadrature<0> &, + SparseMatrix &, + const FunctionMap &, + Vector &, + std::vector &, + const Function<1> * const , + const IteratorRange &, + Threads::ThreadMutex &) +{ + Assert (false, ExcNotImplemented()); +}; + +#endif + + +template +void +MatrixCreator:: +create_boundary_mass_matrix_1 (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const FunctionMap &boundary_functions, + Vector &rhs_vector, + std::vector &dof_to_boundary_mapping, + const Function * const a, + const IteratorRange &range, + Threads::ThreadMutex &mutex) { const FiniteElement &fe = dof.get_fe(); const unsigned int n_components = fe.n_components(); @@ -325,9 +467,8 @@ void MatrixCreator::create_boundary_mass_matrix (const Mapping std::vector dofs_on_face_vector (dofs_per_face); std::set dofs_on_face; - DoFHandler::active_cell_iterator cell = dof.begin_active (), - endc = dof.end (); - for (; cell!=endc; ++cell) + DoFHandler::active_cell_iterator cell = range.first; + for (; cell!=range.second; ++cell) for (unsigned int face=0; face::faces_per_cell; ++face) // check if this face is on that part of // the boundary we are interested in @@ -516,7 +657,8 @@ void MatrixCreator::create_boundary_mass_matrix (const Mapping //in the innermost loop, we traverse a set twice in "find", but //this could be made much faster, e.g. by first building a vector //that already stores the result of find beforehand -// (then remove dofs_on_face altogether) +// (then remove dofs_on_face altogether) + mutex.acquire (); for (unsigned int i=0; i::create_boundary_mass_matrix (const Mapping cell_matrix(i,j)); else { +//TODO:[?] We assume here that shape functions that have support also on the boundary also have their support point on the boundary (by having their indices in dofs_on_face). This is not true, however, e.g. for discontinuous elements. // compare here for relative // smallness Assert (fabs(cell_matrix(i,j)) <= 1e-10 * max_diag_entry, @@ -542,6 +685,7 @@ void MatrixCreator::create_boundary_mass_matrix (const Mapping Assert (fabs(cell_vector(j)) <= 1e-10 * max_diag_entry, ExcInternalError()); }; + mutex.release (); }; }; @@ -561,14 +705,47 @@ void MatrixCreator::create_boundary_mass_matrix (const DoFHandler } -// TODO:[RH, WB] extend this function to use vector valued coefficient functions for system elements. -// TODO:[WB] implement multithreading for this function. + template void MatrixCreator::create_laplace_matrix (const Mapping &mapping, const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, const Function * const coefficient) +{ + const unsigned int n_threads = multithread_info.n_default_threads; + Threads::ThreadManager thread_manager; + + // define starting and end point + // for each thread + std::vector thread_ranges + = Threads::split_range (dof.begin_active(), + dof.end(), n_threads); + + // mutex to synchronise access to + // the matrix + Threads::ThreadMutex mutex; + + // then assemble in parallel + for (unsigned int thread=0; thread::create_laplace_matrix_1) + .collect_args (mapping, dof, q, matrix, coefficient, + thread_ranges[thread], mutex)); + thread_manager.wait (); +}; + + + +// TODO:[RH] extend this function to use vector valued coefficient functions for system elements. +template +void MatrixCreator::create_laplace_matrix_1 (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function * const coefficient, + const IteratorRange &range, + Threads::ThreadMutex &mutex) { UpdateFlags update_flags = UpdateFlags(update_JxW_values | update_gradients); @@ -587,8 +764,8 @@ void MatrixCreator::create_laplace_matrix (const Mapping &mappin std::vector dof_indices (dofs_per_cell); - typename DoFHandler::active_cell_iterator cell = dof.begin_active(); - for (; cell!=dof.end(); ++cell) + typename DoFHandler::active_cell_iterator cell = range.first; + for (; cell!=range.second; ++cell) { fe_values.reinit (cell); @@ -627,10 +804,15 @@ void MatrixCreator::create_laplace_matrix (const Mapping &mappin // transfer everything into the // global object + mutex.acquire (); for (unsigned int i=0; i::create_level_laplace_matrix (unsigned int level, -// TODO:[RH, WB] extend this function to use vector valued coefficient functions for system elements. -// TODO:[WB] implement multithreading for this function. template void MatrixCreator::create_laplace_matrix (const Mapping &mapping, const DoFHandler &dof, @@ -694,6 +874,45 @@ void MatrixCreator::create_laplace_matrix (const Mapping &mappin const Function &rhs, Vector &rhs_vector, const Function * const coefficient) +{ + const unsigned int n_threads = multithread_info.n_default_threads; + Threads::ThreadManager thread_manager; + + // define starting and end point + // for each thread + std::vector thread_ranges + = Threads::split_range (dof.begin_active(), + dof.end(), n_threads); + + // mutex to synchronise access to + // the matrix + Threads::ThreadMutex mutex; + + // then assemble in parallel + for (unsigned int thread=0; thread:: + create_laplace_matrix_2) + .collect_args (mapping, dof, q, matrix, rhs, + rhs_vector, coefficient, + thread_ranges[thread], mutex)); + thread_manager.wait (); +}; + + + +// TODO:[RH] extend this function to use vector valued coefficient functions for system elements. +template +void +MatrixCreator::create_laplace_matrix_2 (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs, + Vector &rhs_vector, + const Function * const coefficient, + const IteratorRange &range, + Threads::ThreadMutex &mutex) { UpdateFlags update_flags = UpdateFlags(update_values | update_gradients | @@ -716,8 +935,8 @@ void MatrixCreator::create_laplace_matrix (const Mapping &mappin std::vector dof_indices (dofs_per_cell); - typename DoFHandler::active_cell_iterator cell = dof.begin_active(); - for (; cell!=dof.end(); ++cell) + typename DoFHandler::active_cell_iterator cell = range.first; + for (; cell!=range.second; ++cell) { fe_values.reinit (cell); @@ -769,12 +988,17 @@ void MatrixCreator::create_laplace_matrix (const Mapping &mappin // transfer everything into the // global object + mutex.acquire (); for (unsigned int i=0; i::create_laplace_matrix (const DoFHandler &dof, } + template template void -- 2.39.5