From: Martin Kronbichler Date: Sat, 20 Apr 2013 16:29:05 +0000 (+0000) Subject: Support BlockVector in matrix free class: access individual blocks of these vectors... X-Git-Tag: v8.0.0~649 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e5efae06177bcaaf03cd26a0c17b7cfc6eb7f035;p=dealii.git Support BlockVector in matrix free class: access individual blocks of these vectors (it was possible also before to directly access a BlockVector with operator(), but that is very inefficient). Add/update tests. git-svn-id: https://svn.dealii.org/trunk@29350 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/matrix_free/fe_evaluation.h b/deal.II/include/deal.II/matrix_free/fe_evaluation.h index 3a5b84c092..b3387178f9 100644 --- a/deal.II/include/deal.II/matrix_free/fe_evaluation.h +++ b/deal.II/include/deal.II/matrix_free/fe_evaluation.h @@ -1781,6 +1781,36 @@ namespace internal } }; + // allows to select between block vectors and non-block vectors, which + // allows to use a unified interface for extracting blocks on block vectors + // and doing nothing on usual vectors + template + struct BlockVectorSelector {}; + + template + struct BlockVectorSelector + { + typedef typename VectorType::BlockType BaseVectorType; + + static BaseVectorType* get_vector_component (VectorType &vec, + const unsigned int component) + { + AssertIndexRange (component, vec.n_blocks()); + return &vec.block(component); + } + }; + + template + struct BlockVectorSelector + { + typedef VectorType BaseVectorType; + + static BaseVectorType* get_vector_component (VectorType &vec, + const unsigned int) + { + return &vec; + } + }; } @@ -2110,13 +2140,12 @@ void FEEvaluationBase ::read_dof_values (const VectorType &src) { - AssertDimension (n_components_, n_fe_components); - - // only need one component, but to silent compiler warnings, use - // n_components copies here (but these will not be used) - VectorType *src_data[n_components]; + // select between block vectors and non-block vectors. Note that the number + // of components is checked in the internal data + typename internal::BlockVectorSelector::value>::BaseVectorType *src_data[n_components]; for (unsigned int d=0; d(&src); + src_data[d] = internal::BlockVectorSelector::value>::get_vector_component(const_cast(src), d); internal::VectorReader reader; read_write_operation (reader, src_data); @@ -2194,13 +2223,13 @@ void FEEvaluationBase ::read_dof_values_plain (const VectorType &src) { - AssertDimension (n_components_, n_fe_components); - - // only need one component, but to avoid compiler warnings, use n_components - // copies here (but these will not be used) - const VectorType *src_data[n_components]; + // select between block vectors and non-block vectors. Note that the number + // of components is checked in the internal data + typename internal::BlockVectorSelector::value>::BaseVectorType *src_data[n_components]; for (unsigned int d=0; d::value>::get_vector_component(src, d); + read_dof_values_plain (src_data); } @@ -2258,15 +2287,15 @@ void FEEvaluationBase ::distribute_local_to_global (VectorType &dst) const { - AssertDimension (n_components_, n_fe_components); Assert (dof_values_initialized==true, internal::ExcAccessToUninitializedField()); - // only need one component, but to avoid compiler warnings, use n_components - // copies here (but these will not be used) - VectorType *dst_data [n_components]; + // select between block vectors and non-block vectors. Note that the number + // of components is checked in the internal data + typename internal::BlockVectorSelector::value>::BaseVectorType *dst_data[n_components]; for (unsigned int d=0; d::value>::get_vector_component(dst, d); internal::VectorDistributorLocalToGlobal distributor; read_write_operation (distributor, dst_data); @@ -2336,15 +2365,15 @@ void FEEvaluationBase ::set_dof_values (VectorType &dst) const { - AssertDimension (n_components_, n_fe_components); Assert (dof_values_initialized==true, internal::ExcAccessToUninitializedField()); - // only need one component, but to avoid compiler warnings, use n_components - // copies here (but these will not be used) - VectorType *dst_data [n_components]; + // select between block vectors and non-block vectors. Note that the number + // of components is checked in the internal data + typename internal::BlockVectorSelector::value>::BaseVectorType *dst_data[n_components]; for (unsigned int d=0; d::value>::get_vector_component(dst, d); internal::VectorSetter setter; read_write_operation (setter, dst_data); diff --git a/deal.II/include/deal.II/matrix_free/matrix_free.h b/deal.II/include/deal.II/matrix_free/matrix_free.h index 2ec5a1ac7f..fa52dc196f 100644 --- a/deal.II/include/deal.II/matrix_free/matrix_free.h +++ b/deal.II/include/deal.II/matrix_free/matrix_free.h @@ -18,11 +18,13 @@ #include #include #include +#include #include #include #include #include #include +#include #include #include #include @@ -616,7 +618,7 @@ public: * different MPI processors need to be accessed at a certain time when * accessing remote data and overlapping communication with computation. */ - typename DoFHandler::active_cell_iterator + typename DoFHandler::cell_iterator get_cell_iterator (const unsigned int macro_cell_number, const unsigned int vector_number, const unsigned int fe_component = 0) const; @@ -1191,7 +1193,7 @@ MatrixFree::get_dof_handler (const unsigned int dof_index) const template inline -typename DoFHandler::active_cell_iterator +typename DoFHandler::cell_iterator MatrixFree::get_cell_iterator(const unsigned int macro_cell_number, const unsigned int vector_number, const unsigned int dof_index) const @@ -1222,7 +1224,7 @@ MatrixFree::get_cell_iterator(const unsigned int macro_cell_number, std::pair index = cell_level_index[macro_cell_number*vectorization_length+vector_number]; - return typename DoFHandler::active_cell_iterator + return typename DoFHandler::cell_iterator (&dofh->get_tria(), index.first, index.second, dofh); } @@ -1638,160 +1640,253 @@ reinit(const Mapping &mapping, // internal helper functions that define how to call MPI data exchange // functions: for generic vectors, do nothing at all. For distributed vectors, // can call update_ghost_values_start function and so on. If we have -// collections of vectors, just do the individual functions of the components +// collections of vectors, just do the individual functions of the +// components. this is a bit messy for block vectors, which use some +// additional helper functions to select the blocks namespace internal { + template + void update_ghost_values_start_block (const VectorStruct &vec, + const unsigned int channel, + internal::bool2type); + template + void update_ghost_values_finish_block (const VectorStruct &vec, + internal::bool2type); + template + void compress_start_block (const VectorStruct &vec, + const unsigned int channel, + internal::bool2type); + template + void compress_finish_block (const VectorStruct &vec, + internal::bool2type); + + template + void update_ghost_values_start_block (const VectorStruct &, + const unsigned int, + internal::bool2type) + {} + template + void update_ghost_values_finish_block (const VectorStruct &, + internal::bool2type) + {} + template + void compress_start_block (const VectorStruct &, + const unsigned int, + internal::bool2type) + {} + template + void compress_finish_block (const VectorStruct &, + internal::bool2type) + {} + + + template inline - void update_ghost_values_start (const VectorStruct &, + void update_ghost_values_start (const VectorStruct &vec, const unsigned int channel = 0) { - (void)channel; + update_ghost_values_start_block(vec, channel, + internal::bool2type::value>()); } template inline - void update_ghost_values_start (const parallel::distributed::Vector &src, + void update_ghost_values_start (const parallel::distributed::Vector &vec, const unsigned int channel = 0) { - src.update_ghost_values_start(channel); + vec.update_ghost_values_start(channel); } template inline - void update_ghost_values_start (const std::vector &src) + void update_ghost_values_start (const std::vector &vec) { - for (unsigned int comp=0; comp inline - void update_ghost_values_start (const std::vector &src) + void update_ghost_values_start (const std::vector &vec) { - for (unsigned int comp=0; comp + inline + void update_ghost_values_start_block (const VectorStruct &vec, + const unsigned int channel, + internal::bool2type) + { + for (unsigned int i=0; i inline - void update_ghost_values_finish (const VectorStruct &) - {} + void update_ghost_values_finish (const VectorStruct &vec) + { + update_ghost_values_finish_block(vec, + internal::bool2type::value>()); + } template inline - void update_ghost_values_finish (const parallel::distributed::Vector &src) + void update_ghost_values_finish (const parallel::distributed::Vector &vec) { - src.update_ghost_values_finish(); + vec.update_ghost_values_finish(); } template inline - void update_ghost_values_finish (const std::vector &src) + void update_ghost_values_finish (const std::vector &vec) { - for (unsigned int comp=0; comp inline - void update_ghost_values_finish (const std::vector &src) + void update_ghost_values_finish (const std::vector &vec) { - for (unsigned int comp=0; comp inline - void compress_start (VectorStruct &, + void update_ghost_values_finish_block (const VectorStruct &vec, + internal::bool2type) + { + for (unsigned int i=0; i + inline + void compress_start (VectorStruct &vec, const unsigned int channel = 0) { - (void)channel; + compress_start_block (vec, channel, + internal::bool2type::value>()); } template inline - void compress_start (parallel::distributed::Vector &dst, + void compress_start (parallel::distributed::Vector &vec, const unsigned int channel = 0) { - dst.compress_start(channel); + vec.compress_start(channel); } template inline - void compress_start (std::vector &dst) + void compress_start (std::vector &vec) { - for (unsigned int comp=0; comp inline - void compress_start (std::vector &dst) + void compress_start (std::vector &vec) { - for (unsigned int comp=0; comp inline - void compress_finish (VectorStruct &) - {} + void compress_start_block (VectorStruct &vec, + const unsigned int channel, + internal::bool2type) + { + for (unsigned int i=0; i + inline + void compress_finish (VectorStruct &vec) + { + compress_finish_block(vec, + internal::bool2type::value>()); + } template inline - void compress_finish (parallel::distributed::Vector &dst) + void compress_finish (parallel::distributed::Vector &vec) + { + vec.compress_finish(); + } + + + + template + inline + void compress_finish (std::vector &vec) { - dst.compress_finish(); + for (unsigned int comp=0; comp inline - void compress_finish (std::vector &dst) + void compress_finish (std::vector &vec) { - for (unsigned int comp=0; comp inline - void compress_finish (std::vector &dst) + void compress_finish_block (VectorStruct &vec, + internal::bool2type) { - for (unsigned int comp=0; comp::local_apply, this, dst, src); }; @@ -128,9 +126,8 @@ void test () BlockSparsityPattern sparsity_pattern; BlockSparseMatrix system_matrix; - BlockVector solution; + BlockVector solution, mf_solution; BlockVector system_rhs; - std::vector > vec1, vec2; dof_handler.distribute_dofs (fe); dof_handler_sca.distribute_dofs (fe_sca); @@ -160,16 +157,7 @@ void test () solution.collect_sizes (); system_rhs.reinit (solution); - - vec1.resize (dim); - vec2.resize (dim); - vec1[0].reinit (dofs_per_block); - vec2[0].reinit (vec1[0]); - for (unsigned int i=1; i > VectorType; - MatrixFreeTest mf (mf_data); - mf.vmult (vec2, vec1); + MatrixFreeTest > mf (mf_data); + mf.vmult (mf_solution, system_rhs); // Verification - double error = 0.; - for (unsigned int i=0; i::local_apply, this, dst, src); }; @@ -144,7 +142,7 @@ void test () BlockVector solution; BlockVector system_rhs; - std::vector > vec1, vec2; + BlockVector mf_solution; dof_handler.distribute_dofs (fe); dof_handler_u.distribute_dofs (fe_u); @@ -268,14 +266,7 @@ void test () solution.collect_sizes (); system_rhs.reinit (solution); - - vec1.resize (2); - vec2.resize (2); - for (unsigned int d=0; d<2; ++d) - { - vec1[d].reinit (dofs_per_block[d]); - vec2[d].reinit (vec1[d]); - } + mf_solution.reinit (solution); // fill system_rhs with random numbers for (unsigned int j=0; j > VectorType; - MatrixFreeTest mf (mf_data); - mf.vmult (vec2, vec1); + MatrixFreeTest > mf (mf_data); + mf.vmult (mf_solution, system_rhs); // Verification - double error = 0.; - for (unsigned int i=0; i<2; ++i) - for (unsigned int j=0; j &data, { fe_eval.reinit (cell); - // compare values with the ones the FEValues - // gives us. Those are seen as reference fe_eval.read_dof_values (src); fe_eval.evaluate (true, true, false); for (unsigned int q=0; q &data, { fe_eval.reinit (cell); - // compare values with the ones the FEValues - // gives us. Those are seen as reference fe_eval.read_dof_values (src); fe_eval.evaluate (true, true, false); for (unsigned int q=0; q &data, { fe_eval.reinit (cell); - // compare values with the ones the FEValues - // gives us. Those are seen as reference fe_eval.read_dof_values (src); fe_eval.evaluate (true, true, false); for (unsigned int q=0; q +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +template +void +helmholtz_operator (const MatrixFree &data, + parallel::distributed::BlockVector &dst, + const parallel::distributed::BlockVector &src, + const std::pair &cell_range) +{ + FEEvaluation fe_eval (data); + const unsigned int n_q_points = fe_eval.n_q_points; + + for(unsigned int cell=cell_range.first;cell +class MatrixFreeTest +{ + public: + typedef VectorizedArray vector_t; + static const std::size_t n_vectors = VectorizedArray::n_array_elements; + + MatrixFreeTest(const MatrixFree &data_in): + data (data_in) + {}; + + void vmult (parallel::distributed::BlockVector &dst, + const parallel::distributed::BlockVector &src) const + { + dst = 0; + const std_cxx1x::function &, + parallel::distributed::BlockVector &, + const parallel::distributed::BlockVector &, + const std::pair&)> + wrap = helmholtz_operator; + data.cell_loop (wrap, dst, src); + }; + +private: + const MatrixFree &data; +}; + + + + + +template +void test () +{ + typedef double number; + + parallel::distributed::Triangulation tria (MPI_COMM_WORLD); + GridGenerator::hyper_cube (tria); + tria.refine_global(1); + typename Triangulation::active_cell_iterator + cell = tria.begin_active (), + endc = tria.end(); + cell = tria.begin_active (); + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) + if (cell->center().norm()<0.2) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + if (dim < 3 && fe_degree < 2) + tria.refine_global(2); + else + tria.refine_global(1); + if (tria.begin(tria.n_levels()-1)->is_locally_owned()) + tria.begin(tria.n_levels()-1)->set_refine_flag(); + if (tria.last()->is_locally_owned()) + tria.last()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + cell = tria.begin_active (); + for (unsigned int i=0; i<10-3*dim; ++i) + { + cell = tria.begin_active (); + unsigned int counter = 0; + for (; cell!=endc; ++cell, ++counter) + if (cell->is_locally_owned()) + if (counter % (7-i) == 0) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + + FE_Q fe (fe_degree); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + + IndexSet owned_set = dof.locally_owned_dofs(); + IndexSet relevant_set; + DoFTools::extract_locally_relevant_dofs (dof, relevant_set); + + ConstraintMatrix constraints (relevant_set); + DoFTools::make_hanging_node_constraints(dof, constraints); + VectorTools::interpolate_boundary_values (dof, 0, ZeroFunction(), + constraints); + constraints.close(); + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + //std::cout << "Number of cells: " << tria.n_global_active_cells() << std::endl; + //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl; + //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl; + + MatrixFree mf_data; + { + const QGauss<1> quad (fe_degree+1); + typename MatrixFree::AdditionalData data; + data.mpi_communicator = MPI_COMM_WORLD; + data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + data.tasks_block_size = 7; + mf_data.reinit (dof, constraints, quad, data); + } + + MatrixFreeTest mf (mf_data); + parallel::distributed::Vector ref; + parallel::distributed::BlockVector in(2), out(2); + for (unsigned int i=0; i<2; ++i) + { + mf_data.initialize_dof_vector (in.block(i)); + mf_data.initialize_dof_vector (out.block(i)); + } + in.collect_sizes(); + out.collect_sizes(); + + mf_data.initialize_dof_vector (ref); + + for (unsigned int i=0; i quadrature_formula(fe_degree+1); + + FEValues fe_values (dof.get_fe(), quadrature_formula, + update_values | update_gradients | + update_JxW_values); + + const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) + { + cell_matrix = 0; + fe_values.reinit (cell); + + for (unsigned int q_point=0; q_pointget_dof_indices(local_dof_indices); + constraints.distribute_local_to_global (cell_matrix, + local_dof_indices, + sparse_matrix); + } + } + sparse_matrix.compress(VectorOperation::add); + + deallog << "Norm of difference (component 1/2): "; + for (unsigned int i=0; i<2; ++i) + { + sparse_matrix.vmult (ref, in.block(i)); + out.block(i) -= ref; + const double diff_norm = out.block(i).linfty_norm(); + deallog << diff_norm << " "; + } + deallog << std::endl << std::endl; +} + + +int main (int argc, char ** argv) +{ + Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv); + + unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + deallog.push(Utilities::int_to_string(myid)); + + if (myid == 0) + { + std::ofstream logfile(output_file_for_mpi("matrix_free_04").c_str()); + deallog.attach(logfile); + deallog << std::setprecision(4); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + deallog.push("2d"); + test<2,1>(); + test<2,2>(); + deallog.pop(); + + deallog.push("3d"); + test<3,1>(); + test<3,2>(); + deallog.pop(); + } + else + { + deallog.depth_console(0); + test<2,1>(); + test<2,2>(); + test<3,1>(); + test<3,2>(); + } +} + diff --git a/tests/mpi/matrix_free_04/ncpu_1/cmp/generic b/tests/mpi/matrix_free_04/ncpu_1/cmp/generic new file mode 100644 index 0000000000..c81a56551e --- /dev/null +++ b/tests/mpi/matrix_free_04/ncpu_1/cmp/generic @@ -0,0 +1,13 @@ + +DEAL:0:2d::Testing FE_Q<2>(1) +DEAL:0:2d::Norm of difference (component 1/2): 0 0 +DEAL:0:2d:: +DEAL:0:2d::Testing FE_Q<2>(2) +DEAL:0:2d::Norm of difference (component 1/2): 0 0 +DEAL:0:2d:: +DEAL:0:3d::Testing FE_Q<3>(1) +DEAL:0:3d::Norm of difference (component 1/2): 0 0 +DEAL:0:3d:: +DEAL:0:3d::Testing FE_Q<3>(2) +DEAL:0:3d::Norm of difference (component 1/2): 0 0 +DEAL:0:3d:: diff --git a/tests/mpi/matrix_free_04/ncpu_2/cmp/generic b/tests/mpi/matrix_free_04/ncpu_2/cmp/generic new file mode 100644 index 0000000000..c81a56551e --- /dev/null +++ b/tests/mpi/matrix_free_04/ncpu_2/cmp/generic @@ -0,0 +1,13 @@ + +DEAL:0:2d::Testing FE_Q<2>(1) +DEAL:0:2d::Norm of difference (component 1/2): 0 0 +DEAL:0:2d:: +DEAL:0:2d::Testing FE_Q<2>(2) +DEAL:0:2d::Norm of difference (component 1/2): 0 0 +DEAL:0:2d:: +DEAL:0:3d::Testing FE_Q<3>(1) +DEAL:0:3d::Norm of difference (component 1/2): 0 0 +DEAL:0:3d:: +DEAL:0:3d::Testing FE_Q<3>(2) +DEAL:0:3d::Norm of difference (component 1/2): 0 0 +DEAL:0:3d::