From: Peter Munch Date: Sat, 1 Oct 2022 16:53:01 +0000 (+0200) Subject: Precompute inverse diagonal (FDM) X-Git-Tag: v9.5.0-rc1~193^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c1daba4287b68a2461775418a2b926dda80018d7;p=dealii.git Precompute inverse diagonal (FDM) --- diff --git a/include/deal.II/lac/tensor_product_matrix.h b/include/deal.II/lac/tensor_product_matrix.h index 47e6a1aa77..dcb542fffb 100644 --- a/include/deal.II/lac/tensor_product_matrix.h +++ b/include/deal.II/lac/tensor_product_matrix.h @@ -376,12 +376,18 @@ public: /** * Constructor. */ - AdditionalData(const bool compress_matrices = true); + AdditionalData(const bool compress_matrices = true, + const bool precompute_inverse_diagonal = true); /** * Try to compress internal matrices. Default: true. */ bool compress_matrices; + + /** + * Precompute inverse diagonal. + */ + bool precompute_inverse_diagonal; }; /** @@ -445,6 +451,11 @@ private: */ const bool compress_matrices; + /** + * Precompute inverse diagonal. + */ + const bool precompute_inverse_diagonal; + /** * Container used to collect 1d matrices if no compression is * requested. The memory is freed during finalize(). @@ -490,6 +501,11 @@ private: */ AlignedVector eigenvalues; + /** + * Vector of inverted eigenvalues. + */ + AlignedVector inverted_eigenvalues; + /** * Pointer into mass_matrices, derivative_matrices, and eigenvalues. */ @@ -499,6 +515,11 @@ private: * Pointer into mass_matrices, derivative_matrices, and eigenvalues. */ std::vector matrix_ptr; + + /** + * Number of rows in 1 of each cell. + */ + std::vector vector_n_rows_1d; }; @@ -810,7 +831,8 @@ namespace internal AlignedVector &tmp, const unsigned int n_rows_1d_non_templated, const std::array &eigenvectors, - const std::array &eigenvalues) + const std::array &eigenvalues, + const Number *inverted_eigenvalues = nullptr) { const unsigned int n_rows_1d = n_rows_1d_templated == 0 ? n_rows_1d_non_templated : @@ -836,8 +858,13 @@ namespace internal { const Number *S = eigenvectors[0]; eval.template apply<0, true, false>(S, src, t); + for (unsigned int i = 0; i < n_rows_1d; ++i) - t[i] /= eigenvalues[0][i]; + if (inverted_eigenvalues) + t[i] *= inverted_eigenvalues[i]; + else + t[i] /= eigenvalues[0][i]; + eval.template apply<0, false, false>(S, t, dst); } @@ -847,9 +874,14 @@ namespace internal const Number *S1 = eigenvectors[1]; eval.template apply<0, true, false>(S0, src, t); eval.template apply<1, true, false>(S1, t, dst); + for (unsigned int i1 = 0, c = 0; i1 < n_rows_1d; ++i1) for (unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c) - dst[c] /= (eigenvalues[1][i1] + eigenvalues[0][i0]); + if (inverted_eigenvalues) + dst[c] *= inverted_eigenvalues[c]; + else + dst[c] /= (eigenvalues[1][i1] + eigenvalues[0][i0]); + eval.template apply<0, false, false>(S0, dst, t); eval.template apply<1, false, false>(S1, t, dst); } @@ -862,11 +894,16 @@ namespace internal eval.template apply<0, true, false>(S0, src, t); eval.template apply<1, true, false>(S1, t, dst); eval.template apply<2, true, false>(S2, dst, t); + for (unsigned int i2 = 0, c = 0; i2 < n_rows_1d; ++i2) for (unsigned int i1 = 0; i1 < n_rows_1d; ++i1) for (unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c) - t[c] /= (eigenvalues[2][i2] + eigenvalues[1][i1] + - eigenvalues[0][i0]); + if (inverted_eigenvalues) + t[c] *= inverted_eigenvalues[c]; + else + t[c] /= (eigenvalues[2][i2] + eigenvalues[1][i1] + + eigenvalues[0][i0]); + eval.template apply<0, false, false>(S0, t, dst); eval.template apply<1, false, false>(S1, dst, t); eval.template apply<2, false, false>(S2, t, dst); @@ -896,7 +933,8 @@ namespace internal AlignedVector & tmp, const unsigned int n_rows_1d, const std::array &eigenvectors, - const std::array &eigenvalues); + const std::array &eigenvalues, + const Number *inverted_eigenvalues = nullptr); } // namespace TensorProductMatrixSymmetricSum } // namespace internal @@ -1073,8 +1111,10 @@ TensorProductMatrixSymmetricSum::reinit( template TensorProductMatrixSymmetricSumCollection:: - AdditionalData::AdditionalData(const bool compress_matrices) + AdditionalData::AdditionalData(const bool compress_matrices, + const bool precompute_inverse_diagonal) : compress_matrices(compress_matrices) + , precompute_inverse_diagonal(precompute_inverse_diagonal) {} @@ -1084,6 +1124,7 @@ TensorProductMatrixSymmetricSumCollection:: TensorProductMatrixSymmetricSumCollection( const AdditionalData &additional_data) : compress_matrices(additional_data.compress_matrices) + , precompute_inverse_diagonal(additional_data.precompute_inverse_diagonal) {} @@ -1248,7 +1289,7 @@ TensorProductMatrixSymmetricSumCollection::finalize() } else { - // case 2) compression requested but none possible + // case 3) compress this->vector_ptr.resize(cache.size() + 1); this->matrix_ptr.resize(cache.size() + 1); @@ -1277,6 +1318,137 @@ TensorProductMatrixSymmetricSumCollection::finalize() cache.clear(); } + + if (precompute_inverse_diagonal) + { + if (dim == 1) + { + // 1D case: simply invert 1D eigenvalues + for (unsigned int i = 0; i < this->eigenvalues.size(); ++i) + this->eigenvalues[i] = Number(1.0) / this->eigenvalues[i]; + std::swap(this->inverted_eigenvalues, eigenvalues); + } + else + { + // 2D and 3D case: we have 2 or 3 1d eigenvalues so that we + // need to combine these + + // step 1) if eigenvalues/eigenvectors are compressed, we + // need to compress the diagonal (the combination of ev + // indices) as well. This is an optional step. + std::vector indices_ev; + + if (indices.size() > 0) + { + // 1a) create cache (ev indics -> diag index) + const unsigned int n_cells = indices.size() / dim; + std::map, unsigned int> cache_ev; + std::vector cache_ev_idx(n_cells); + + for (unsigned int i = 0, c = 0; i < n_cells; ++i) + { + std::array id; + + for (unsigned int d = 0; d < dim; ++d, ++c) + id[d] = indices[c]; + + const auto id_ptr = cache_ev.find(id); + + if (id_ptr == cache_ev.end()) + { + const auto size = cache_ev.size(); + cache_ev_idx[i] = size; + cache_ev[id] = size; + } + else + { + cache_ev_idx[i] = id_ptr->second; + } + } + + // 1b) store diagonal indices for each cell + std::vector new_indices; + new_indices.reserve(indices.size() / dim * (dim + 1)); + + for (unsigned int i = 0, c = 0; i < n_cells; ++i) + { + for (unsigned int d = 0; d < dim; ++d, ++c) + new_indices.push_back(indices[c]); + new_indices.push_back(cache_ev_idx[i]); + } + + // 1c) transpose cache (diag index -> ev indices) + indices_ev.resize(cache_ev.size() * dim); + for (const auto &entry : cache_ev) + for (unsigned int d = 0; d < dim; ++d) + indices_ev[entry.second * dim + d] = entry.first[d]; + + std::swap(this->indices, new_indices); + } + + // step 2) allocate memory and set pointers + const unsigned int n_diag = + ((indices_ev.size() > 0) ? indices_ev.size() : + (matrix_ptr.size() - 1)) / + dim; + + std::vector new_vector_ptr(n_diag + 1, 0); + std::vector new_vector_n_rows_1d(n_diag, 0); + + for (unsigned int i = 0; i < n_diag; ++i) + { + const unsigned int c = (indices_ev.size() > 0) ? + indices_ev[dim * i + 0] : + (dim * i + 0); + + const unsigned int n_rows = vector_ptr[c + 1] - vector_ptr[c]; + + new_vector_n_rows_1d[i] = n_rows; + new_vector_ptr[i + 1] = Utilities::pow(n_rows, dim); + } + + for (unsigned int i = 0; i < n_diag; ++i) + new_vector_ptr[i + 1] += new_vector_ptr[i]; + + this->inverted_eigenvalues.resize(new_vector_ptr.back()); + + // step 3) loop over all unique diagonal entries and invert + for (unsigned int i = 0; i < n_diag; ++i) + { + std::array evs; + + for (unsigned int d = 0; d < dim; ++d) + evs[d] = + &this + ->eigenvalues[this->vector_ptr[(indices_ev.size() > 0) ? + indices_ev[dim * i + d] : + (dim * i + d)]]; + + const unsigned int mm = new_vector_n_rows_1d[i]; + if (dim == 2) + { + for (unsigned int i1 = 0, c = 0; i1 < mm; ++i1) + for (unsigned int i0 = 0; i0 < mm; ++i0, ++c) + this->inverted_eigenvalues[new_vector_ptr[i] + c] = + Number(1.0) / (evs[1][i1] + evs[0][i0]); + } + else + { + for (unsigned int i2 = 0, c = 0; i2 < mm; ++i2) + for (unsigned int i1 = 0; i1 < mm; ++i1) + for (unsigned int i0 = 0; i0 < mm; ++i0, ++c) + this->inverted_eigenvalues[new_vector_ptr[i] + c] = + Number(1.0) / (evs[2][i2] + evs[1][i1] + evs[0][i0]); + } + } + + // step 4) clean up + std::swap(this->vector_ptr, new_vector_ptr); + std::swap(this->vector_n_rows_1d, new_vector_n_rows_1d); + } + + this->eigenvalues.clear(); + } } @@ -1292,28 +1464,92 @@ TensorProductMatrixSymmetricSumCollection:: Number * dst = dst_in.begin(); const Number *src = src_in.begin(); - std::array eigenvectors, eigenvalues; - unsigned int n_rows_1d_non_templated = 0; - - for (unsigned int d = 0; d < dim; ++d) + if (this->eigenvalues.empty() == false) { - const unsigned int translated_index = - (indices.size() > 0) ? indices[dim * index + d] : (dim * index + d); - - eigenvectors[d] = - this->eigenvectors.data() + matrix_ptr[translated_index]; - eigenvalues[d] = this->eigenvalues.data() + vector_ptr[translated_index]; - n_rows_1d_non_templated = - vector_ptr[translated_index + 1] - vector_ptr[translated_index]; - } + std::array eigenvectors; + std::array eigenvalues; + unsigned int n_rows_1d_non_templated = 0; - if (n_rows_1d != -1) - internal::TensorProductMatrixSymmetricSum::apply_inverse< - n_rows_1d == -1 ? 0 : n_rows_1d>( - dst, src, tmp_array, n_rows_1d_non_templated, eigenvectors, eigenvalues); + for (unsigned int d = 0; d < dim; ++d) + { + const unsigned int translated_index = + (indices.size() > 0) ? indices[dim * index + d] : (dim * index + d); + + eigenvectors[d] = + this->eigenvectors.data() + matrix_ptr[translated_index]; + eigenvalues[d] = + this->eigenvalues.data() + vector_ptr[translated_index]; + n_rows_1d_non_templated = + vector_ptr[translated_index + 1] - vector_ptr[translated_index]; + } + + if (n_rows_1d != -1) + internal::TensorProductMatrixSymmetricSum::apply_inverse< + n_rows_1d == -1 ? 0 : n_rows_1d>(dst, + src, + tmp_array, + n_rows_1d_non_templated, + eigenvectors, + eigenvalues); + else + internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>( + dst, + src, + tmp_array, + n_rows_1d_non_templated, + eigenvectors, + eigenvalues); + } else - internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>( - dst, src, tmp_array, n_rows_1d_non_templated, eigenvectors, eigenvalues); + { + std::array eigenvectors; + const Number * inverted_eigenvalues = nullptr; + unsigned int n_rows_1d_non_templated = 0; + + for (unsigned int d = 0; d < dim; ++d) + { + const unsigned int translated_index = + (indices.size() > 0) ? + indices[((dim == 1) ? 1 : (dim + 1)) * index + d] : + (dim * index + d); + + eigenvectors[d] = + this->eigenvectors.data() + matrix_ptr[translated_index]; + } + + { + const unsigned int translated_index = + ((indices.size() > 0) && (dim != 1)) ? + indices[(dim + 1) * index + dim] : + index; + + inverted_eigenvalues = + this->inverted_eigenvalues.data() + vector_ptr[translated_index]; + n_rows_1d_non_templated = + (dim == 1) ? + (vector_ptr[translated_index + 1] - vector_ptr[translated_index]) : + vector_n_rows_1d[translated_index]; + } + + if (n_rows_1d != -1) + internal::TensorProductMatrixSymmetricSum::apply_inverse< + n_rows_1d == -1 ? 0 : n_rows_1d>(dst, + src, + tmp_array, + n_rows_1d_non_templated, + eigenvectors, + {}, + inverted_eigenvalues); + else + internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>( + dst, + src, + tmp_array, + n_rows_1d_non_templated, + eigenvectors, + {}, + inverted_eigenvalues); + } } diff --git a/include/deal.II/lac/tensor_product_matrix.templates.h b/include/deal.II/lac/tensor_product_matrix.templates.h index c7b026fea0..90a16f2e32 100644 --- a/include/deal.II/lac/tensor_product_matrix.templates.h +++ b/include/deal.II/lac/tensor_product_matrix.templates.h @@ -69,16 +69,34 @@ namespace internal AlignedVector & tmp, const unsigned int n_rows_1d, const std::array &eigenvectors, - const std::array &eigenvalues) + const std::array &eigenvalues, + const Number *inverted_eigenvalues) { if (n_rows_1d_templated == n_rows_1d) - apply_inverse( - dst, src, tmp, n_rows_1d, eigenvectors, eigenvalues); + apply_inverse(dst, + src, + tmp, + n_rows_1d, + eigenvectors, + eigenvalues, + inverted_eigenvalues); else if (n_rows_1d_templated < FDM_N_ROWS_MAX) select_apply_inverse( - dst, src, tmp, n_rows_1d, eigenvectors, eigenvalues); + dst, + src, + tmp, + n_rows_1d, + eigenvectors, + eigenvalues, + inverted_eigenvalues); else - apply_inverse<0>(dst, src, tmp, n_rows_1d, eigenvectors, eigenvalues); + apply_inverse<0>(dst, + src, + tmp, + n_rows_1d, + eigenvectors, + eigenvalues, + inverted_eigenvalues); } } // namespace TensorProductMatrixSymmetricSum } // namespace internal diff --git a/source/lac/tensor_product_matrix.inst.in b/source/lac/tensor_product_matrix.inst.in index 6607de8aca..ec5c1c41ca 100644 --- a/source/lac/tensor_product_matrix.inst.in +++ b/source/lac/tensor_product_matrix.inst.in @@ -35,7 +35,8 @@ for (deal_II_dimension : DIMENSIONS; const std::array &eigenvectors, const std::array - &eigenvalues); + & eigenvalues, + const deal_II_scalar_vectorized *inverted_eigenvalues); } for (deal_II_dimension : DIMENSIONS; deal_II_scalar : REAL_SCALARS) @@ -55,5 +56,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_scalar : REAL_SCALARS) AlignedVector & tmp, const unsigned int n_rows, const std::array &eigenvectors, - const std::array &eigenvalues); + const std::array &eigenvalues, + const deal_II_scalar *inverted_eigenvalues); } diff --git a/tests/lac/tensor_product_matrix_08.cc b/tests/lac/tensor_product_matrix_08.cc index 508869c3b9..bbe162923a 100644 --- a/tests/lac/tensor_product_matrix_08.cc +++ b/tests/lac/tensor_product_matrix_08.cc @@ -53,11 +53,15 @@ do_test_mesh(const Mapping &mapping, const Triangulation &tria) const auto harmonic_patch_extent = GridTools::compute_harmonic_patch_extent(mapping, tria, quadrature_face); - FDM collection_0(typename FDM::AdditionalData(true)); - FDM collection_1(typename FDM::AdditionalData(false)); + FDM collection_0(typename FDM::AdditionalData(true, false)); + FDM collection_1(typename FDM::AdditionalData(false, false)); + FDM collection_2(typename FDM::AdditionalData(true, true)); + FDM collection_3(typename FDM::AdditionalData(false, true)); collection_0.reserve(tria.n_active_cells()); collection_1.reserve(tria.n_active_cells()); + collection_2.reserve(tria.n_active_cells()); + collection_3.reserve(tria.n_active_cells()); for (const auto &cell : tria.active_cell_iterators()) { @@ -80,10 +84,18 @@ do_test_mesh(const Mapping &mapping, const Triangulation &tria) collection_1.insert(cell->active_cell_index(), M_and_K.first, M_and_K.second); + collection_2.insert(cell->active_cell_index(), + M_and_K.first, + M_and_K.second); + collection_3.insert(cell->active_cell_index(), + M_and_K.first, + M_and_K.second); } collection_0.finalize(); collection_1.finalize(); + collection_2.finalize(); + collection_3.finalize(); deallog << "Storage sizes: " << collection_0.storage_size() << " " << collection_1.storage_size() << std::endl; @@ -93,6 +105,8 @@ do_test_mesh(const Mapping &mapping, const Triangulation &tria) AlignedVector tmp; FullMatrix matrix_0(fe.n_dofs_per_cell(), fe.n_dofs_per_cell()); FullMatrix matrix_1(fe.n_dofs_per_cell(), fe.n_dofs_per_cell()); + FullMatrix matrix_2(fe.n_dofs_per_cell(), fe.n_dofs_per_cell()); + FullMatrix matrix_3(fe.n_dofs_per_cell(), fe.n_dofs_per_cell()); for (unsigned int cell = 0; cell < tria.n_active_cells(); ++cell) { @@ -114,14 +128,34 @@ do_test_mesh(const Mapping &mapping, const Triangulation &tria) tmp); for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j) matrix_1[j][i] = dst[j]; + + collection_2.apply_inverse(cell, + make_array_view(dst), + make_array_view(src), + tmp); + for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j) + matrix_2[j][i] = dst[j]; + + collection_3.apply_inverse(cell, + make_array_view(dst), + make_array_view(src), + tmp); + for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j) + matrix_3[j][i] = dst[j]; } - FloatingPointComparator comp(1e-5, false); + FloatingPointComparator comp(1e-5, true); Assert((comp.compare(matrix_0, matrix_1) == FloatingPointComparator::ComparisonResult::equal), ExcInternalError()); + Assert((comp.compare(matrix_0, matrix_2) == + FloatingPointComparator::ComparisonResult::equal), + ExcInternalError()); + Assert((comp.compare(matrix_0, matrix_3) == + FloatingPointComparator::ComparisonResult::equal), + ExcInternalError()); } deallog << "OK!" << std::endl;