From: Peter Munch Date: Wed, 9 Dec 2020 12:49:21 +0000 (+0100) Subject: Implement MatrixFreeTools::compute_matrix() X-Git-Tag: v9.3.0-rc1~770^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ac9ae5b09fcd1451a063c015352fd8838828bded;p=dealii.git Implement MatrixFreeTools::compute_matrix() --- diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h index 6ff02f3950..4d46c108c7 100644 --- a/include/deal.II/matrix_free/tools.h +++ b/include/deal.II/matrix_free/tools.h @@ -72,6 +72,64 @@ namespace MatrixFreeTools const unsigned int first_selected_component = 0); + /** + * Compute the matrix representation of a linear operator (@p matrix), given + * @p matrix_free and the local cell integral operation @p local_vmult. The + * vector is initialized to the right size in the function. + * + * The parameters @p dof_no, @p quad_no, and @p first_selected_component are + * passed to the constructor of the FEEvaluation that is internally set up. + */ + template + void + compute_matrix( + const MatrixFree & matrix_free, + const AffineConstraints & constraints, + MatrixType & matrix, + const std::function &)> &local_vmult, + const unsigned int dof_no = 0, + const unsigned int quad_no = 0, + const unsigned int first_selected_component = 0); + + /** + * Same as above but with a class and a function pointer. + */ + template + void + compute_matrix( + const MatrixFree &matrix_free, + const AffineConstraints & constraints, + MatrixType & matrix, + void (CLASS::*cell_operation)(FEEvaluation &) const, + const CLASS * owning_class, + const unsigned int dof_no = 0, + const unsigned int quad_no = 0, + const unsigned int first_selected_component = 0); + + // implementations @@ -509,6 +567,156 @@ namespace MatrixFreeTools false); } + template + void + compute_matrix( + const MatrixFree & matrix_free, + const AffineConstraints & constraints, + MatrixType & matrix, + const std::function &)> &local_vmult, + const unsigned int dof_no, + const unsigned int quad_no, + const unsigned int first_selected_component) + { + matrix_free.template cell_loop( + [&](const auto &, auto &dst, const auto &, const auto range) { + FEEvaluation + integrator( + matrix_free, range, dof_no, quad_no, first_selected_component); + + unsigned int const dofs_per_cell = integrator.dofs_per_cell; + + std::vector dof_indices(dofs_per_cell); + std::vector dof_indices_mf(dofs_per_cell); + + std::array, + VectorizedArrayType::size()> + matrices; + + std::fill_n(matrices.begin(), + VectorizedArrayType::size(), + FullMatrix(dofs_per_cell, + dofs_per_cell)); + + const auto lexicographic_numbering = + matrix_free + .get_shape_info(dof_no, + quad_no, + first_selected_component, + integrator.get_active_fe_index(), + integrator.get_active_quadrature_index()) + .lexicographic_numbering; + + for (auto cell = range.first; cell < range.second; ++cell) + { + integrator.reinit(cell); + + unsigned int const n_filled_lanes = + matrix_free.n_active_entries_per_cell_batch(cell); + + for (unsigned int v = 0; v < n_filled_lanes; ++v) + matrices[v] = 0.0; + + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + integrator.begin_dof_values()[i] = + static_cast(i == j); + + local_vmult(integrator); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int v = 0; v < n_filled_lanes; ++v) + matrices[v](i, j) = integrator.begin_dof_values()[i][v]; + } + + for (unsigned int v = 0; v < n_filled_lanes; ++v) + { + const auto cell_v = matrix_free.get_cell_iterator(cell, v); + + if (matrix_free.get_mg_level() != numbers::invalid_unsigned_int) + cell_v->get_mg_dof_indices(dof_indices); + else + cell_v->get_dof_indices(dof_indices); + + for (unsigned int j = 0; j < dof_indices.size(); ++j) + dof_indices_mf[j] = dof_indices[lexicographic_numbering[j]]; + + constraints.distribute_local_to_global(matrices[v], + dof_indices_mf, + dof_indices_mf, + dst); + } + } + }, + matrix, + matrix); + + matrix.compress(VectorOperation::add); + + for (unsigned int i = 0; i < matrix.m(); ++i) + if (matrix.get_row_length(i) > 0 && matrix(i, i) == 0.0 && + constraints.is_constrained(i)) + matrix.add(i, i, 1); + } + + template + void + compute_matrix( + const MatrixFree &matrix_free, + const AffineConstraints & constraints, + MatrixType & matrix, + void (CLASS::*cell_operation)(FEEvaluation &) const, + const CLASS * owning_class, + const unsigned int dof_no, + const unsigned int quad_no, + const unsigned int first_selected_component) + { + compute_matrix(matrix_free, + constraints, + matrix, + [&](auto &feeval) { + (owning_class->*cell_operation)(feeval); + }, + dof_no, + quad_no, + first_selected_component); + } + } // namespace MatrixFreeTools DEAL_II_NAMESPACE_CLOSE diff --git a/tests/matrix_free/compute_diagonal_01.cc b/tests/matrix_free/compute_diagonal_01.cc index 5013d77e72..1a9d46ba9c 100644 --- a/tests/matrix_free/compute_diagonal_01.cc +++ b/tests/matrix_free/compute_diagonal_01.cc @@ -65,6 +65,7 @@ test() Test test(matrix_free, + constraint, [](FEEvaluation test(matrix_free, + constraints, [](FEEvaluation test(matrix_free, + constraint, [](FEEvaluation #include +#include #include #include @@ -30,6 +31,9 @@ #include #include +#include +#include + #include #include #include @@ -55,7 +59,8 @@ class Test using VectorType = LinearAlgebra::distributed::Vector; public: - Test(const MatrixFree matrix_free, + Test(const MatrixFree &matrix_free, + const AffineConstraints & constraints, const std::function &)> &cell_operation) : matrix_free(matrix_free) + , constraints(constraints) , cell_operation(cell_operation) {} @@ -71,8 +77,27 @@ public: do_test() { // compute diagonal with the new function + VectorType diagonal_global, diagonal_global_reference; + + SparseMatrix A1, A2, A_ref; + SparsityPattern sparsity_pattern; + + const bool test_matrix = Utilities::MPI::n_mpi_processes( + matrix_free.get_task_info().communicator) == 1; + + if (test_matrix) + { + DynamicSparsityPattern dsp(matrix_free.get_dof_handler().n_dofs()); + DoFTools::make_sparsity_pattern(matrix_free.get_dof_handler(), + dsp, + constraints); + sparsity_pattern.copy_from(dsp); + A1.reinit(sparsity_pattern); + A2.reinit(sparsity_pattern); + A_ref.reinit(sparsity_pattern); + } + { - VectorType diagonal_global; MatrixFreeTools::compute_diagonal>( + matrix_free, constraints, A1, [&](auto &phi) { + this->cell_operation(phi); + }); + } + + if (test_matrix) + { + MatrixFreeTools::compute_matrix( + matrix_free, constraints, A2, &Test::cell_function, this); + } + // compute diagonal globally { - VectorType src, dst, temp; + VectorType src, temp; matrix_free.initialize_dof_vector(src); - matrix_free.initialize_dof_vector(dst); + matrix_free.initialize_dof_vector(diagonal_global_reference); matrix_free.initialize_dof_vector(temp); for (unsigned int i = 0; i < src.size(); i++) @@ -104,15 +149,40 @@ public: if (src.get_partitioner()->in_local_range(i)) { - dst[i] = temp[i]; - src[i] = 0.0; + diagonal_global_reference[i] = temp[i]; + src[i] = 0.0; } + + if (test_matrix) + { + for (unsigned int j = 0; j < src.size(); j++) + if (temp[j] != 0.0) + A_ref(j, i) = temp[j]; + else if (i == j) + A_ref(j, i) = 1.0; + } + temp = 0.0; } - dst.print(deallog.get_file_stream()); - deallog << dst.l2_norm() << std::endl; + diagonal_global_reference.print(deallog.get_file_stream()); + deallog << diagonal_global_reference.l2_norm() << std::endl; } + + if (test_matrix) + { + auto a1 = A1.begin(); + auto a2 = A2.begin(); + auto a_ref = A_ref.begin(); + + for (; a1 != A1.end(); ++a1, ++a2, ++a_ref) + { + Assert(std::abs(a1->value() - a_ref->value()) < 1e-6, + ExcNotImplemented()); + Assert(std::abs(a2->value() - a_ref->value()) < 1e-6, + ExcNotImplemented()); + } + } } void @@ -137,7 +207,19 @@ public: } } - const MatrixFree matrix_free; + void + cell_function(FEEvaluation &phi) const + { + this->cell_operation(phi); + } + + const MatrixFree &matrix_free; + const AffineConstraints & constraints; const std::function