From 5918079cdf64d9285fa6581861add0872c12ce34 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Thu, 1 Oct 2009 12:00:20 +0000 Subject: [PATCH] Change the parallelization scheme: Now we use WorkStream instead of parallel_for with locks. Now this tutorial program should be finished. Need to read through the comments once again, though. git-svn-id: https://svn.dealii.org/trunk@19635 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-37/doc/intro.dox | 8 +- deal.II/examples/step-37/doc/results.dox | 25 +- deal.II/examples/step-37/step-37.cc | 460 ++++++++++++++--------- 3 files changed, 295 insertions(+), 198 deletions(-) diff --git a/deal.II/examples/step-37/doc/intro.dox b/deal.II/examples/step-37/doc/intro.dox index f8e8367b47..456a02c1a9 100644 --- a/deal.II/examples/step-37/doc/intro.dox +++ b/deal.II/examples/step-37/doc/intro.dox @@ -431,8 +431,12 @@ matrix-vector product implementation efficient on a GPU. For our program, we choose to follow a simple strategy: We let several processors work together by splitting the cells into several chunks. The -threading building blocks implemenation of a %parallel for loop -implements this concept using the parallel::apply_to_subranges() function. +threading building blocks implemenation of a %parallel pipeline implements +this concept using the WorkStream::run() function. What the pipeline does is +pretty much the same as a for loop. However, it can be instructed to do some +part of the loop body by just one process at a time and in natural order. We +need to use this for writing the local contributions into the global vector, +in order to avoid a race condition.

Combination with multigrid

diff --git a/deal.II/examples/step-37/doc/results.dox b/deal.II/examples/step-37/doc/results.dox index a444ddf67a..e8bce13790 100644 --- a/deal.II/examples/step-37/doc/results.dox +++ b/deal.II/examples/step-37/doc/results.dox @@ -9,20 +9,20 @@ output: @code Cycle 0 Number of degrees of freedom: 337 -System matrix memory consumption: 0.0257 MBytes. -Multigrid objects memory consumption: 0.05071 MBytes. +System matrix memory consumption: 0.02573 MBytes. +Multigrid objects memory consumption: 0.05083 MBytes. Convergence in 10 CG iterations. Cycle 1 Number of degrees of freedom: 1313 -System matrix memory consumption: 0.0925 MBytes. -Multigrid objects memory consumption: 0.1792 MBytes. +System matrix memory consumption: 0.09254 MBytes. +Multigrid objects memory consumption: 0.1794 MBytes. Convergence in 10 CG iterations. Cycle 2 Number of degrees of freedom: 5185 -System matrix memory consumption: 0.3551 MBytes. -Multigrid objects memory consumption: 0.6775 MBytes. +System matrix memory consumption: 0.3552 MBytes. +Multigrid objects memory consumption: 0.6777 MBytes. Convergence in 10 CG iterations. Cycle 3 @@ -40,7 +40,7 @@ Convergence in 10 CG iterations. Cycle 5 Number of degrees of freedom: 328193 System matrix memory consumption: 22.1 MBytes. -Multigrid objects memory consumption: 41.63 MBytes. +Multigrid objects memory consumption: 41.64 MBytes. Convergence in 10 CG iterations. @endcode @@ -51,14 +51,14 @@ program in three dimensions: @code Cycle 0 Number of degrees of freedom: 517 -System matrix memory consumption: 0.1 MBytes. -Multigrid objects memory consumption: 0.1462 MBytes. +System matrix memory consumption: 0.1001 MBytes. +Multigrid objects memory consumption: 0.1463 MBytes. Convergence in 9 CG iterations. Cycle 1 Number of degrees of freedom: 3817 -System matrix memory consumption: 0.6612 MBytes. -Multigrid objects memory consumption: 0.8894 MBytes. +System matrix memory consumption: 0.6613 MBytes. +Multigrid objects memory consumption: 0.8896 MBytes. Convergence in 10 CG iterations. Cycle 2 @@ -75,9 +75,8 @@ Convergence in 11 CG iterations. Cycle 4 Number of degrees of freedom: 1847617 -System matrix memory consumption: 321.9 MBytes. +System matrix memory consumption: 322 MBytes. Multigrid objects memory consumption: 415.1 MBytes. -Convergence in 11 CG iterations. @endcode A comparison of the memory consumption of the MatrixFree class with a diff --git a/deal.II/examples/step-37/step-37.cc b/deal.II/examples/step-37/step-37.cc index 1d3a3edd84..c6e7d51f82 100644 --- a/deal.II/examples/step-37/step-37.cc +++ b/deal.II/examples/step-37/step-37.cc @@ -16,6 +16,7 @@ #include #include #include +#include #include #include @@ -109,6 +110,57 @@ void Coefficient::value_list (const std::vector > &points, // @sect3{Matrix-free implementation.} + // First com a few variables that we use for + // defining the %parallel layout of the vector + // multiplication function with the WorkStream + // concept in the Matrix-free class. +namespace WorkStreamData +{ + template + struct ScratchData + { + ScratchData (); + ScratchData (const ScratchData &scratch); + FullMatrix solutions; + }; + + template + ScratchData::ScratchData () + : + solutions () + {} + + template + ScratchData::ScratchData (const ScratchData &scratch) + : + solutions () + {} + + template + struct CopyData : public ScratchData + { + CopyData (); + CopyData (const CopyData &scratch); + unsigned int first_dof; + unsigned int n_dofs; + }; + + template + CopyData::CopyData () + : + ScratchData () + {} + + template + CopyData::CopyData (const CopyData &scratch) + : + ScratchData () + {} + +} + + + // Next comes the implemenation of the // matrix-free class. It provides some // standard information we expect for @@ -196,11 +248,18 @@ class MatrixFree : public Subscriptor // as well as a few other variables that // store matrix properties. private: + typedef std::vector >::const_iterator + CellChunkIterator; template - void vmult_on_subrange (const unsigned int first_cell, - const unsigned int last_cell, - Vector &dst, - const Vector &src) const; + void local_vmult (CellChunkIterator cell_range, + WorkStreamData::ScratchData &scratch, + WorkStreamData::CopyData ©, + const Vector &src) const; + + template + void + copy_local_to_global (const WorkStreamData::CopyData ©, + Vector &dst) const; FullMatrix small_matrix; Table<2,unsigned int> indices_local_to_global; @@ -216,6 +275,7 @@ class MatrixFree : public Subscriptor unsigned int n_dofs, n_cells; unsigned int m, n; unsigned int n_points, n_comp; + std::vector > chunks; } matrix_sizes; }; @@ -232,81 +292,6 @@ MatrixFree::MatrixFree () - // This function initializes the structures - // of the matrix. It writes the number of - // total degrees of freedom in the problem - // as well as the number of cells to the - // MatrixSizes struct and copies the small - // matrix that transforms the solution from - // support points to quadrature points. It - // uses the small matrix for determining - // the number of degrees of freedom per - // cell (number of rows in - // small_matrix). The number - // of quadrature points needs to be passed - // through the last variable - // n_points_per_cell, since - // the number of columns in the small - // matrix is - // dim*n_points_per_cell for - // the Laplace problem (the Laplacian is a - // tensor and has dim - // components). In this function, we also - // give the fields containing the - // derivative information and the local dof - // indices the correct sizes. They will be - // filled by calling the respective set - // function. -template -void MatrixFree:: -reinit (const unsigned int n_dofs_in, - const unsigned int n_cells_in, - const FullMatrix &small_matrix_in, - const unsigned int n_points_per_cell) -{ - small_matrix = small_matrix_in; - - derivatives.reinit (n_cells_in, n_points_per_cell); - indices_local_to_global.reinit (n_cells_in, small_matrix.m()); - - diagonal_is_calculated = false; - - matrix_sizes.n_dofs = n_dofs_in; - matrix_sizes.n_cells = n_cells_in; - matrix_sizes.m = small_matrix.m(); - matrix_sizes.n = small_matrix.n(); - matrix_sizes.n_points = n_points_per_cell; - matrix_sizes.n_comp = small_matrix.n()/matrix_sizes.n_points; - Assert(matrix_sizes.n_comp * n_points_per_cell == small_matrix.n(), - ExcInternalError()); -} - - - - // This function we need if we want to - // delete the content of the matrix, - // e.g. when we are finished with one grid - // level and continue to the next one. Just - // put all the field sizes to 0. -template -void -MatrixFree::clear () -{ - small_matrix.reinit(0,0); - derivatives.reinit (0,0); - indices_local_to_global.reinit(0,0); - - constraints.clear(); - - diagonal_values.reinit (0); - diagonal_is_calculated = false; - - matrix_sizes.n_dofs = 0; - matrix_sizes.n_cells = 0; -} - - - // This function returns the number of rows // of the global matrix, and the next one // the number of columns (which is the @@ -403,13 +388,14 @@ set_derivative_data (const unsigned int cell_no, // This is the central function of the // matrix-free class, implementing the // multiplication of the matrix with a - // vector. This function actually not work - // on all the cells, but only a subset of - // cells. Since this function operates - // similarly irrespective on which cell - // chunk we are sitting, we can parallelize - // it and get very regular operation - // patterns. + // vector. This function does actually not + // work on all the cells, but only a subset + // of cells, specified by the first argument + // cell_range. Since this + // function operates similarly irrespective + // on which cell chunk we are sitting, we can + // parallelize it and get very regular + // operation patterns. // // Following the discussion in the // introduction, we try to work on multiple @@ -445,65 +431,17 @@ set_derivative_data (const unsigned int cell_no, // cell for the first and the number of // quadrature points times the number of // components per point for the latter. - // - // One more thing to make this work - // efficiently is to decide how many cells - // should be included in the matrix that - // contains the solution values at local - // dofs for several cells. If we choose too - // few cells, then the gains from using the - // matrix-matrix product will not be fully - // utilized (dgemm tends to provide more - // efficiency the larger the matrix - // dimensions get). If we choose too many, - // we will firstly degrade parallelization - // (which is based on some these chunks), - // and secondly introduce an inefficiency - // that comes from the computer - // architecture: Right after the first - // matrix-matrix multiplication, we - // transform the solution on quadrature - // points by using derivatives. Obviously, - // we want to have fast access to that - // data, so it should still be present in - // L2 cache and not to be fetched from main - // memory. The total memory usage of the - // data on quadrature points should be not - // more than about half the cache size of - // the processor in order to be on the safe - // side. Since most today's processors - // provide 512 kBytes or more cache memory - // per core, we choose about 250 kB as a - // size. Clearly, this is an - // architecture-dependent value and the - // interested user can squeeze out some - // extra performance by hand-tuning this - // parameter. Once we have chosen the - // number of cells we collect in one chunk, - // we determine how many chunks we have on - // the given cell range and recalculate the - // actual chunk size in order to evenly - // distribute the chunks. template template void MatrixFree:: -vmult_on_subrange (const unsigned int first_cell, - const unsigned int last_cell, - Vector &dst, - const Vector &src) const + local_vmult (CellChunkIterator cell_range, + WorkStreamData::ScratchData &scratch, + WorkStreamData::CopyData ©, + const Vector &src) const { - FullMatrix solution_cells, solution_points; - - const unsigned int divisor = 250000/(matrix_sizes.n*sizeof(number)); - const unsigned int n_chunks = (last_cell-first_cell)/divisor + 1; - const unsigned int chunk_size = - (last_cell-first_cell)/n_chunks + ((last_cell-first_cell)%n_chunks>0); - - for (unsigned int k=first_cell; klast_cell ? last_cell-k : chunk_size; + const unsigned int first_cell = cell_range->first, + chunk_size = cell_range->second - cell_range->first; // OK, now we are sitting in the loop that // goes over our chunks of cells. What we @@ -559,33 +497,49 @@ vmult_on_subrange (const unsigned int first_cell, // are related to each other. Since we // simultaneously apply the constraints, we // hand this task off to the ConstraintMatrix - // object. Most often, itis used to work on - // one cell at a time, but since we work on a - // whole chunk of dofs, we can do that just - // as easily for all the cells at once. - solution_cells.reinit (current_chunk_size,matrix_sizes.m, true); - solution_points.reinit (current_chunk_size,matrix_sizes.n, true); - - const unsigned int n_cell_entries = current_chunk_size*matrix_sizes.m; - constraints.get_dof_values(src, &indices_local_to_global(k,0), - &solution_cells(0,0), - &solution_cells(0,0)+n_cell_entries); - - solution_cells.mmult (solution_points, small_matrix); - - for (unsigned int i=0; icopy_local_to_global + // function). Most often, the + // ConstraintMatrix function is used to be + // applied to data from one cell at a time, + // but since we work on a whole chunk of + // dofs, we can do that just as easily for + // all the cells at once. + copy.solutions.reinit (chunk_size,matrix_sizes.m, true); + copy.first_dof = first_cell; + copy.n_dofs = chunk_size*matrix_sizes.m; + scratch.solutions.reinit (chunk_size,matrix_sizes.n, true); + + constraints.get_dof_values(src, &indices_local_to_global(copy.first_dof,0), + ©.solutions(0,0), + ©.solutions(0,0)+copy.n_dofs); + + copy.solutions.mmult (scratch.solutions, small_matrix); + + for (unsigned int i=0, k = first_cell; i +template +void +MatrixFree:: + copy_local_to_global (const WorkStreamData::CopyData ©, + Vector &dst) const +{ + constraints.distribute_local_to_global (©.solutions(0,0), + ©.solutions(0,0)+copy.n_dofs, + &indices_local_to_global(copy.first_dof,0), + dst); } @@ -640,27 +594,29 @@ MatrixFree::Tvmult_add (Vector &dst, // vector src and adds the // result to vector dst. We call // a %parallel function that applies the - // multiplication on a subrange of cells - // (cf. the @ref threads module). - // - // TODO: Use WorkStream for parallelization - // instead of apply_to_subranges, once we - // have realized the best way for doing - // that. + // multiplication on a chunk of cells at once + // using the WorkStream module (cf. also the + // @ref threads module). The subdivision into + // chunks will be performed in the reinit + // function and is stored in the field + // matrix_sizes.chunks. template template void MatrixFree::vmult_add (Vector &dst, const Vector &src) const { - parallel::apply_to_subranges (0, matrix_sizes.n_cells, - std_cxx1x::bind(&MatrixFree:: - template vmult_on_subrange, - this, - _1,_2, - boost::ref(dst), - boost::cref(src)), - 200); + + WorkStream::run (matrix_sizes.chunks.begin(), matrix_sizes.chunks.end(), + std_cxx1x::bind(&MatrixFree:: + template local_vmult, + this, _1, _2, _3, boost::cref(src)), + std_cxx1x::bind(&MatrixFree:: + template copy_local_to_global, + this, _1, boost::ref(dst)), + WorkStreamData::ScratchData(), + WorkStreamData::CopyData(), + 4,1); // One thing to be cautious about: The // deal.II classes expect that the matrix @@ -683,6 +639,142 @@ MatrixFree::vmult_add (Vector &dst, + // This function initializes the structures + // of the matrix. It writes the number of + // total degrees of freedom in the problem + // as well as the number of cells to the + // MatrixSizes struct and copies the small + // matrix that transforms the solution from + // support points to quadrature points. It + // uses the small matrix for determining + // the number of degrees of freedom per + // cell (number of rows in + // small_matrix). The number + // of quadrature points needs to be passed + // through the last variable + // n_points_per_cell, since + // the number of columns in the small + // matrix is + // dim*n_points_per_cell for + // the Laplace problem (the Laplacian is a + // tensor and has dim + // components). In this function, we also + // give the fields containing the + // derivative information and the local dof + // indices the correct sizes. They will be + // filled by calling the respective set + // function. +template +void MatrixFree:: +reinit (const unsigned int n_dofs_in, + const unsigned int n_cells_in, + const FullMatrix &small_matrix_in, + const unsigned int n_points_per_cell) +{ + small_matrix = small_matrix_in; + + derivatives.reinit (n_cells_in, n_points_per_cell); + indices_local_to_global.reinit (n_cells_in, small_matrix.m()); + + diagonal_is_calculated = false; + + matrix_sizes.n_dofs = n_dofs_in; + matrix_sizes.n_cells = n_cells_in; + matrix_sizes.m = small_matrix.m(); + matrix_sizes.n = small_matrix.n(); + matrix_sizes.n_points = n_points_per_cell; + matrix_sizes.n_comp = small_matrix.n()/matrix_sizes.n_points; + Assert(matrix_sizes.n_comp * n_points_per_cell == small_matrix.n(), + ExcInternalError()); + + // One thing to make the matrix-vector + // product with this class efficient is to + // decide how many cells should be summarized + // to one chunk, which will determine the + // size of the full matrix that we work + // on. If we choose too few cells, then the + // gains from using the matrix-matrix product + // will not be fully utilized (dgemm tends to + // provide more efficiency the larger the + // matrix dimensions get). If we choose too + // many, we will firstly degrade + // parallelization (which is based on some + // these chunks), and secondly introduce an + // inefficiency that comes from the computer + // architecture: In the actual working + // function above, right after the first + // matrix-matrix multiplication, we transform + // the solution on quadrature points by using + // derivatives. Obviously, we want to have + // fast access to that data, so it should + // still be present in L2 cache and not + // needed to be fetched from main memory. The + // total memory usage of the data on + // quadrature points should be not more than + // about half the cache size of the processor + // in order to be on the safe side. Since + // most today's processors provide 512 kBytes + // or more cache memory per core, we choose + // about 50 kB as a size to be on the safe + // side (other things need to be stored in + // the CPU as well). Clearly, this is an + // architecture-dependent value and the + // interested user can squeeze out some extra + // performance by hand-tuning this + // parameter. Once we have chosen the number + // of cells we collect in one chunk, we + // determine how many chunks we have on the + // given cell range and recalculate the + // actual chunk size in order to evenly + // distribute the chunks. + const unsigned int divisor = 50000/(matrix_sizes.m*sizeof(double)); + unsigned int n_chunks = matrix_sizes.n_cells/divisor + 1; + if (n_chunks<2*multithread_info.n_default_threads) + n_chunks = 2*multithread_info.n_default_threads; + + const unsigned int chunk_size = (matrix_sizes.n_cells/n_chunks + + (matrix_sizes.n_cells%n_chunks>0)); + matrix_sizes.chunks.resize (n_chunks); + for (unsigned int i=0; i matrix_sizes.n_cells) + { + matrix_sizes.chunks[i].second = matrix_sizes.n_cells; + break; + } + else + matrix_sizes.chunks[i].second = (i+1)*chunk_size; + } +} + + + + // This function we need if we want to + // delete the content of the matrix, + // e.g. when we are finished with one grid + // level and continue to the next one. Just + // put all the field sizes to 0. +template +void +MatrixFree::clear () +{ + small_matrix.reinit(0,0); + derivatives.reinit (0,0); + indices_local_to_global.reinit(0,0); + + constraints.clear(); + + diagonal_values.reinit (0); + diagonal_is_calculated = false; + + matrix_sizes.n_dofs = 0; + matrix_sizes.n_cells = 0; + matrix_sizes.chunks.clear(); +} + + + // This function returns the entries of the // matrix. Since this class is intended not // to store the matrix entries, it would @@ -778,7 +870,9 @@ std::size_t MatrixFree::memory_consumption () const indices_local_to_global.memory_consumption() + constraints.memory_consumption() + small_matrix.memory_consumption() + - diagonal_values.memory_consumption() + sizeof(*this); + diagonal_values.memory_consumption() + + matrix_sizes.chunks.size()*2*sizeof(unsigned int) + + sizeof(*this); return glob_size; } -- 2.39.5