From c91d75e24e26b6d0730458f3055ed9564225fec8 Mon Sep 17 00:00:00 2001 From: Dominik Still Date: Wed, 6 Nov 2024 15:13:47 +0100 Subject: [PATCH] Enable batched matrix-vector product --- .../deal.II/matrix_free/evaluation_kernels.h | 253 +++++++-- .../matrix_free/evaluation_kernels_face.h | 248 +++++++-- .../matrix_free/tensor_product_kernels.h | 424 ++++++++++++--- .../matrix_free/tensor_product_kernels_01.cc | 134 +++++ .../tensor_product_kernels_01.output | 509 ++++++++++++++++++ tests/simplex/matrix_free_05.cc | 352 ++++++++++++ tests/simplex/matrix_free_05.output | 3 + 7 files changed, 1755 insertions(+), 168 deletions(-) create mode 100644 tests/matrix_free/tensor_product_kernels_01.cc create mode 100644 tests/matrix_free/tensor_product_kernels_01.output create mode 100644 tests/simplex/matrix_free_05.cc create mode 100644 tests/simplex/matrix_free_05.output diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 41644d20f5..8edf57f5cb 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -637,17 +637,41 @@ namespace internal auto *out = fe_eval.begin_values(); const auto *in = values_dofs_actual; - for (unsigned int c = 0; c < n_components; ++c) + for (unsigned int c = 0; c < n_components; c += 3) { - apply_matrix_vector_product( - shape_values, in, out, n_dofs, n_q_points, 1, 1); - - out += n_q_points; - in += n_dofs; + if (c + 1 == n_components) + apply_matrix_vector_product( + shape_values, in, out, n_dofs, n_q_points, 1, 1); + else if (c + 2 == n_components) + apply_matrix_vector_product( + shape_values, in, out, n_dofs, n_q_points, 1, 1); + else + apply_matrix_vector_product( + shape_values, in, out, n_dofs, n_q_points, 1, 1); + + out += 3 * n_q_points; + in += 3 * n_dofs; } } @@ -658,16 +682,41 @@ namespace internal auto *out = fe_eval.begin_gradients(); const auto *in = values_dofs_actual; - for (unsigned int c = 0; c < n_components; ++c) + for (unsigned int c = 0; c < n_components; c += 3) { - apply_matrix_vector_product( - shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); - out += n_q_points * dim; - in += n_dofs; + if (c + 1 == n_components) + apply_matrix_vector_product( + shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + else if (c + 2 == n_components) + apply_matrix_vector_product( + shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + else + apply_matrix_vector_product( + shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + + out += 3 * n_q_points * dim; + in += 3 * n_dofs; } } } @@ -705,25 +754,76 @@ namespace internal auto *in = fe_eval.begin_values(); auto *out = values_dofs_actual; - for (unsigned int c = 0; c < n_components; ++c) + for (unsigned int c = 0; c < n_components; c += 3) { if (add_into_values_array == false) - apply_matrix_vector_product( - shape_values, in, out, n_dofs, n_q_points, 1, 1); + { + if (c + 1 == n_components) + apply_matrix_vector_product( + shape_values, in, out, n_dofs, n_q_points, 1, 1); + else if (c + 2 == n_components) + apply_matrix_vector_product( + shape_values, in, out, n_dofs, n_q_points, 1, 1); + else + apply_matrix_vector_product( + shape_values, in, out, n_dofs, n_q_points, 1, 1); + } else - apply_matrix_vector_product( - shape_values, in, out, n_dofs, n_q_points, 1, 1); - - in += n_q_points; - out += n_dofs; + { + if (c + 1 == n_components) + apply_matrix_vector_product( + shape_values, in, out, n_dofs, n_q_points, 1, 1); + else if (c + 2 == n_components) + apply_matrix_vector_product( + shape_values, in, out, n_dofs, n_q_points, 1, 1); + else + apply_matrix_vector_product( + shape_values, in, out, n_dofs, n_q_points, 1, 1); + } + out += 3 * n_dofs; + in += 3 * n_q_points; } } @@ -734,26 +834,77 @@ namespace internal auto *in = fe_eval.begin_gradients(); auto *out = values_dofs_actual; - for (unsigned int c = 0; c < n_components; ++c) + for (unsigned int c = 0; c < n_components; c += 3) { if (add_into_values_array == false && !(integration_flag & EvaluationFlags::values)) - apply_matrix_vector_product( - shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + { + if (c + 1 == n_components) + apply_matrix_vector_product( + shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + else if (c + 2 == n_components) + apply_matrix_vector_product( + shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + else + apply_matrix_vector_product( + shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + } else - apply_matrix_vector_product( - shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); - - in += n_q_points * dim; - out += n_dofs; + { + if (c + 1 == n_components) + apply_matrix_vector_product( + shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + else if (c + 2 == n_components) + apply_matrix_vector_product( + shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + else + apply_matrix_vector_product( + shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + } + out += 3 * n_dofs; + in += 3 * n_q_points * dim; } } } diff --git a/include/deal.II/matrix_free/evaluation_kernels_face.h b/include/deal.II/matrix_free/evaluation_kernels_face.h index e41dad1b80..e3f2e118bd 100644 --- a/include/deal.II/matrix_free/evaluation_kernels_face.h +++ b/include/deal.II/matrix_free/evaluation_kernels_face.h @@ -1475,17 +1475,41 @@ namespace internal auto *out = fe_eval.begin_values(); auto *in = values_dofs; - for (unsigned int c = 0; c < n_components; ++c) + for (unsigned int c = 0; c < n_components; c += 3) { - apply_matrix_vector_product( - shape_values, in, out, n_dofs, n_q_points, 1, 1); - - out += n_q_points; - in += n_dofs; + if (c + 1 == n_components) + apply_matrix_vector_product( + shape_values, in, out, n_dofs, n_q_points, 1, 1); + else if (c + 2 == n_components) + apply_matrix_vector_product( + shape_values, in, out, n_dofs, n_q_points, 1, 1); + else + apply_matrix_vector_product( + shape_values, in, out, n_dofs, n_q_points, 1, 1); + + out += 3 * n_q_points; + in += 3 * n_dofs; } } @@ -1497,16 +1521,40 @@ namespace internal const auto *const shape_gradients = &shape_data.shape_gradients_face(face_no, face_orientation, 0); - for (unsigned int c = 0; c < n_components; ++c) + for (unsigned int c = 0; c < n_components; c += 3) { - apply_matrix_vector_product( - shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); - out += n_q_points * dim; - in += n_dofs; + if (c + 1 == n_components) + apply_matrix_vector_product( + shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + else if (c + 2 == n_components) + apply_matrix_vector_product( + shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + else + apply_matrix_vector_product( + shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + out += 3 * n_q_points * dim; + in += 3 * n_dofs; } } @@ -1884,24 +1932,76 @@ namespace internal auto *in = fe_eval.begin_values(); auto *out = values_dofs; - for (unsigned int c = 0; c < n_components; ++c) + for (unsigned int c = 0; c < n_components; c += 3) { if (sum_into_values) - apply_matrix_vector_product( - shape_values, in, out, n_dofs, n_q_points, 1, 1); + { + if (c + 1 == n_components) + apply_matrix_vector_product( + shape_values, in, out, n_dofs, n_q_points, 1, 1); + else if (c + 2 == n_components) + apply_matrix_vector_product( + shape_values, in, out, n_dofs, n_q_points, 1, 1); + else + apply_matrix_vector_product( + shape_values, in, out, n_dofs, n_q_points, 1, 1); + } else - apply_matrix_vector_product( - shape_values, in, out, n_dofs, n_q_points, 1, 1); - in += n_q_points; - out += n_dofs; + { + if (c + 1 == n_components) + apply_matrix_vector_product( + shape_values, in, out, n_dofs, n_q_points, 1, 1); + else if (c + 2 == n_components) + apply_matrix_vector_product( + shape_values, in, out, n_dofs, n_q_points, 1, 1); + else + apply_matrix_vector_product( + shape_values, in, out, n_dofs, n_q_points, 1, 1); + } + in += 3 * n_q_points; + out += 3 * n_dofs; } } @@ -1917,21 +2017,73 @@ namespace internal { if (!sum_into_values && !(integration_flag & EvaluationFlags::values)) - apply_matrix_vector_product( - shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + { + if (c + 1 == n_components) + apply_matrix_vector_product( + shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + else if (c + 2 == n_components) + apply_matrix_vector_product( + shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + else + apply_matrix_vector_product( + shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + } else - apply_matrix_vector_product( - shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); - in += n_q_points * dim; - out += n_dofs; + { + if (c + 1 == n_components) + apply_matrix_vector_product( + shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + else if (c + 2 == n_components) + apply_matrix_vector_product( + shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + else + apply_matrix_vector_product( + shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + } + in += 3 * n_q_points * dim; + out += 3 * n_dofs; } } diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index cc7c78f36e..1676a97e93 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -195,7 +195,8 @@ namespace internal bool add, bool consider_strides, typename Number, - typename Number2> + typename Number2, + int n_components = 1> std::enable_if_t<(variant == evaluate_general), void> apply_matrix_vector_product(const Number2 *matrix, const Number *in, @@ -222,12 +223,15 @@ namespace internal const int stride_in = consider_strides ? stride_in_given : 1; const int stride_out = consider_strides ? stride_out_given : 1; + static_assert(n_components > 0 && n_components < 4, + "Invalid number of components"); + // specialization for n_rows = 2 that manually unrolls the innermost loop // to make the operation perform better (not completely as good as the // templated one, but much better than the generic version down below, // because the loop over col can be more effectively unrolled by the // compiler) - if (transpose_matrix && n_rows == 2) + if (transpose_matrix && n_rows == 2 && n_components == 1) { const Number2 *matrix_1 = matrix + n_columns; const Number x0 = in[0], x1 = in[stride_in]; @@ -240,7 +244,7 @@ namespace internal out[stride_out * col] = result; } } - else if (transpose_matrix && n_rows == 3) + else if (transpose_matrix && n_rows == 3 && n_components == 1) { const Number2 *matrix_1 = matrix + n_columns; const Number2 *matrix_2 = matrix_1 + n_columns; @@ -255,7 +259,8 @@ namespace internal out[stride_out * col] = result; } } - else if (std::abs(in - out) < std::min(stride_out * nn, stride_in * mm)) + else if (std::abs(in - out) < std::min(stride_out * nn, stride_in * mm) && + n_components == 1) { Assert(mm <= 128, ExcNotImplemented("For large sizes, arrays may not overlap")); @@ -286,24 +291,70 @@ namespace internal } else { + const Number *in0 = in; + const Number *in1 = n_components > 1 ? in + mm : nullptr; + const Number *in2 = n_components > 2 ? in + 2 * mm : nullptr; + + Number *out0 = out; + Number *out1 = n_components > 1 ? out + nn : nullptr; + Number *out2 = n_components > 2 ? out + 2 * nn : nullptr; + int nn_regular = (nn / 4) * 4; for (int col = 0; col < nn_regular; col += 4) { - Number res0, res1, res2, res3; + Number res[12]; if (transpose_matrix == true) { const Number2 *matrix_ptr = matrix + col; - res0 = matrix_ptr[0] * in[0]; - res1 = matrix_ptr[1] * in[0]; - res2 = matrix_ptr[2] * in[0]; - res3 = matrix_ptr[3] * in[0]; + const Number a = in0[0]; + res[0] = matrix_ptr[0] * a; + res[1] = matrix_ptr[1] * a; + res[2] = matrix_ptr[2] * a; + res[3] = matrix_ptr[3] * a; + + if (n_components > 1) + { + const Number b = in1[0]; + res[4] = matrix_ptr[0] * b; + res[5] = matrix_ptr[1] * b; + res[6] = matrix_ptr[2] * b; + res[7] = matrix_ptr[3] * b; + } + + if (n_components > 2) + { + const Number c = in2[0]; + res[8] = matrix_ptr[0] * c; + res[9] = matrix_ptr[1] * c; + res[10] = matrix_ptr[2] * c; + res[11] = matrix_ptr[3] * c; + } + matrix_ptr += n_columns; for (int i = 1; i < mm; ++i, matrix_ptr += n_columns) { - res0 += matrix_ptr[0] * in[stride_in * i]; - res1 += matrix_ptr[1] * in[stride_in * i]; - res2 += matrix_ptr[2] * in[stride_in * i]; - res3 += matrix_ptr[3] * in[stride_in * i]; + const Number a = in0[stride_in * i]; + res[0] += matrix_ptr[0] * a; + res[1] += matrix_ptr[1] * a; + res[2] += matrix_ptr[2] * a; + res[3] += matrix_ptr[3] * a; + + if (n_components > 1) + { + const Number b = in1[stride_in * i]; + res[4] += matrix_ptr[0] * b; + res[5] += matrix_ptr[1] * b; + res[6] += matrix_ptr[2] * b; + res[7] += matrix_ptr[3] * b; + } + if (n_components > 2) + { + const Number c = in2[stride_in * i]; + res[8] += matrix_ptr[0] * c; + res[9] += matrix_ptr[1] * c; + res[10] += matrix_ptr[2] * c; + res[11] += matrix_ptr[3] * c; + } } } else @@ -313,49 +364,144 @@ namespace internal const Number2 *matrix_2 = matrix + (col + 2) * n_columns; const Number2 *matrix_3 = matrix + (col + 3) * n_columns; - res0 = matrix_0[0] * in[0]; - res1 = matrix_1[0] * in[0]; - res2 = matrix_2[0] * in[0]; - res3 = matrix_3[0] * in[0]; + const Number a = in0[0]; + res[0] = matrix_0[0] * a; + res[1] = matrix_1[0] * a; + res[2] = matrix_2[0] * a; + res[3] = matrix_3[0] * a; + + if (n_components > 1) + { + const Number b = in1[0]; + res[4] = matrix_0[0] * b; + res[5] = matrix_1[0] * b; + res[6] = matrix_2[0] * b; + res[7] = matrix_3[0] * b; + } + + if (n_components > 2) + { + const Number c = in2[0]; + res[8] = matrix_0[0] * c; + res[9] = matrix_1[0] * c; + res[10] = matrix_2[0] * c; + res[11] = matrix_3[0] * c; + } + for (int i = 1; i < mm; ++i) { - res0 += matrix_0[i] * in[stride_in * i]; - res1 += matrix_1[i] * in[stride_in * i]; - res2 += matrix_2[i] * in[stride_in * i]; - res3 += matrix_3[i] * in[stride_in * i]; + const Number a = in0[stride_in * i]; + res[0] += matrix_0[i] * a; + res[1] += matrix_1[i] * a; + res[2] += matrix_2[i] * a; + res[3] += matrix_3[i] * a; + + if (n_components > 1) + { + const Number b = in1[stride_in * i]; + res[4] += matrix_0[i] * b; + res[5] += matrix_1[i] * b; + res[6] += matrix_2[i] * b; + res[7] += matrix_3[i] * b; + } + + if (n_components > 2) + { + const Number c = in2[stride_in * i]; + res[8] += matrix_0[i] * c; + res[9] += matrix_1[i] * c; + res[10] += matrix_2[i] * c; + res[11] += matrix_3[i] * c; + } } } if (add) { - out[0] += res0; - out[stride_out] += res1; - out[2 * stride_out] += res2; - out[3 * stride_out] += res3; + out0[0] += res[0]; + out0[stride_out] += res[1]; + out0[2 * stride_out] += res[2]; + out0[3 * stride_out] += res[3]; + if (n_components > 1) + { + out1[0] += res[4]; + out1[stride_out] += res[5]; + out1[2 * stride_out] += res[6]; + out1[3 * stride_out] += res[7]; + } + if (n_components > 2) + { + out2[0] += res[8]; + out2[stride_out] += res[9]; + out2[2 * stride_out] += res[10]; + out2[3 * stride_out] += res[11]; + } } else { - out[0] = res0; - out[stride_out] = res1; - out[2 * stride_out] = res2; - out[3 * stride_out] = res3; + out0[0] = res[0]; + out0[stride_out] = res[1]; + out0[2 * stride_out] = res[2]; + out0[3 * stride_out] = res[3]; + if (n_components > 1) + { + out1[0] = res[4]; + out1[stride_out] = res[5]; + out1[2 * stride_out] = res[6]; + out1[3 * stride_out] = res[7]; + } + if (n_components > 2) + { + out2[0] = res[8]; + out2[stride_out] = res[9]; + out2[2 * stride_out] = res[10]; + out2[3 * stride_out] = res[11]; + } } - out += 4 * stride_out; + out0 += 4 * stride_out; + if (n_components > 1) + out1 += 4 * stride_out; + if (n_components > 2) + out2 += 4 * stride_out; } if (nn - nn_regular == 3) { - Number res0, res1, res2; + Number res0, res1, res2, res3, res4, res5, res6, res7, res8; if (transpose_matrix == true) { const Number2 *matrix_ptr = matrix + nn_regular; - res0 = matrix_ptr[0] * in[0]; - res1 = matrix_ptr[1] * in[0]; - res2 = matrix_ptr[2] * in[0]; + res0 = matrix_ptr[0] * in0[0]; + res1 = matrix_ptr[1] * in0[0]; + res2 = matrix_ptr[2] * in0[0]; + if (n_components > 1) + { + res3 = matrix_ptr[0] * in1[0]; + res4 = matrix_ptr[1] * in1[0]; + res5 = matrix_ptr[2] * in1[0]; + } + if (n_components > 2) + { + res6 = matrix_ptr[0] * in2[0]; + res7 = matrix_ptr[1] * in2[0]; + res8 = matrix_ptr[2] * in2[0]; + } matrix_ptr += n_columns; for (int i = 1; i < mm; ++i, matrix_ptr += n_columns) { - res0 += matrix_ptr[0] * in[stride_in * i]; - res1 += matrix_ptr[1] * in[stride_in * i]; - res2 += matrix_ptr[2] * in[stride_in * i]; + res0 += matrix_ptr[0] * in0[stride_in * i]; + res1 += matrix_ptr[1] * in0[stride_in * i]; + res2 += matrix_ptr[2] * in0[stride_in * i]; + if (n_components > 1) + { + res3 += matrix_ptr[0] * in1[stride_in * i]; + res4 += matrix_ptr[1] * in1[stride_in * i]; + res5 += matrix_ptr[2] * in1[stride_in * i]; + } + if (n_components > 2) + { + res6 += matrix_ptr[0] * in2[stride_in * i]; + res7 += matrix_ptr[1] * in2[stride_in * i]; + res8 += matrix_ptr[2] * in2[stride_in * i]; + } } } else @@ -364,42 +510,110 @@ namespace internal const Number2 *matrix_1 = matrix + (nn_regular + 1) * n_columns; const Number2 *matrix_2 = matrix + (nn_regular + 2) * n_columns; - res0 = matrix_0[0] * in[0]; - res1 = matrix_1[0] * in[0]; - res2 = matrix_2[0] * in[0]; + res0 = matrix_0[0] * in0[0]; + res1 = matrix_1[0] * in0[0]; + res2 = matrix_2[0] * in0[0]; + if (n_components > 1) + { + res3 = matrix_0[0] * in1[0]; + res4 = matrix_1[0] * in1[0]; + res5 = matrix_2[0] * in1[0]; + } + if (n_components > 2) + { + res6 = matrix_0[0] * in2[0]; + res7 = matrix_1[0] * in2[0]; + res8 = matrix_2[0] * in2[0]; + } for (int i = 1; i < mm; ++i) { - res0 += matrix_0[i] * in[stride_in * i]; - res1 += matrix_1[i] * in[stride_in * i]; - res2 += matrix_2[i] * in[stride_in * i]; + res0 += matrix_0[i] * in0[stride_in * i]; + res1 += matrix_1[i] * in0[stride_in * i]; + res2 += matrix_2[i] * in0[stride_in * i]; + if (n_components > 1) + { + res3 += matrix_0[i] * in1[stride_in * i]; + res4 += matrix_1[i] * in1[stride_in * i]; + res5 += matrix_2[i] * in1[stride_in * i]; + } + if (n_components > 2) + { + res6 += matrix_0[i] * in2[stride_in * i]; + res7 += matrix_1[i] * in2[stride_in * i]; + res8 += matrix_2[i] * in2[stride_in * i]; + } } } if (add) { - out[0] += res0; - out[stride_out] += res1; - out[2 * stride_out] += res2; + out0[0] += res0; + out0[stride_out] += res1; + out0[2 * stride_out] += res2; + if (n_components > 1) + { + out1[0] += res3; + out1[stride_out] += res4; + out1[2 * stride_out] += res5; + } + if (n_components > 2) + { + out2[0] += res6; + out2[stride_out] += res7; + out2[2 * stride_out] += res8; + } } else { - out[0] = res0; - out[stride_out] = res1; - out[2 * stride_out] = res2; + out0[0] = res0; + out0[stride_out] = res1; + out0[2 * stride_out] = res2; + if (n_components > 1) + { + out1[0] = res3; + out1[stride_out] = res4; + out1[2 * stride_out] = res5; + } + if (n_components > 2) + { + out2[0] = res6; + out2[stride_out] = res7; + out2[2 * stride_out] = res8; + } } } else if (nn - nn_regular == 2) { - Number res0, res1; + Number res0, res1, res2, res3, res4, res5; if (transpose_matrix == true) { const Number2 *matrix_ptr = matrix + nn_regular; - res0 = matrix_ptr[0] * in[0]; - res1 = matrix_ptr[1] * in[0]; + res0 = matrix_ptr[0] * in0[0]; + res1 = matrix_ptr[1] * in0[0]; + if (n_components > 1) + { + res2 = matrix_ptr[0] * in1[0]; + res3 = matrix_ptr[1] * in1[0]; + } + if (n_components > 2) + { + res4 = matrix_ptr[0] * in2[0]; + res5 = matrix_ptr[1] * in2[0]; + } matrix_ptr += n_columns; for (int i = 1; i < mm; ++i, matrix_ptr += n_columns) { - res0 += matrix_ptr[0] * in[stride_in * i]; - res1 += matrix_ptr[1] * in[stride_in * i]; + res0 += matrix_ptr[0] * in0[stride_in * i]; + res1 += matrix_ptr[1] * in0[stride_in * i]; + if (n_components > 1) + { + res2 += matrix_ptr[0] * in1[stride_in * i]; + res3 += matrix_ptr[1] * in1[stride_in * i]; + } + if (n_components > 2) + { + res4 += matrix_ptr[0] * in2[stride_in * i]; + res5 += matrix_ptr[1] * in2[stride_in * i]; + } } } else @@ -407,47 +621,119 @@ namespace internal const Number2 *matrix_0 = matrix + nn_regular * n_columns; const Number2 *matrix_1 = matrix + (nn_regular + 1) * n_columns; - res0 = matrix_0[0] * in[0]; - res1 = matrix_1[0] * in[0]; + res0 = matrix_0[0] * in0[0]; + res1 = matrix_1[0] * in0[0]; + if (n_components > 1) + { + res2 = matrix_0[0] * in1[0]; + res3 = matrix_1[0] * in1[0]; + } + if (n_components > 2) + { + res4 = matrix_0[0] * in2[0]; + res5 = matrix_1[0] * in2[0]; + } for (int i = 1; i < mm; ++i) { - res0 += matrix_0[i] * in[stride_in * i]; - res1 += matrix_1[i] * in[stride_in * i]; + res0 += matrix_0[i] * in0[stride_in * i]; + res1 += matrix_1[i] * in0[stride_in * i]; + if (n_components > 1) + { + res2 += matrix_0[i] * in1[stride_in * i]; + res3 += matrix_1[i] * in1[stride_in * i]; + } + if (n_components > 2) + { + res4 += matrix_0[i] * in2[stride_in * i]; + res5 += matrix_1[i] * in2[stride_in * i]; + } } } if (add) { - out[0] += res0; - out[stride_out] += res1; + out0[0] += res0; + out0[stride_out] += res1; + if (n_components > 1) + { + out1[0] += res2; + out1[stride_out] += res3; + } + if (n_components > 2) + { + out2[0] += res4; + out2[stride_out] += res5; + } } else { - out[0] = res0; - out[stride_out] = res1; + out0[0] = res0; + out0[stride_out] = res1; + if (n_components > 1) + { + out1[0] = res2; + out1[stride_out] = res3; + } + if (n_components > 2) + { + out2[0] = res4; + out2[stride_out] = res5; + } } } else if (nn - nn_regular == 1) { - Number res0; + Number res0, res1, res2; if (transpose_matrix == true) { const Number2 *matrix_ptr = matrix + nn_regular; - res0 = matrix_ptr[0] * in[0]; + res0 = matrix_ptr[0] * in0[0]; + if (n_components > 1) + res1 = matrix_ptr[0] * in1[0]; + if (n_components > 2) + res2 = matrix_ptr[0] * in2[0]; matrix_ptr += n_columns; for (int i = 1; i < mm; ++i, matrix_ptr += n_columns) - res0 += matrix_ptr[0] * in[stride_in * i]; + { + res0 += matrix_ptr[0] * in0[stride_in * i]; + if (n_components > 1) + res1 += matrix_ptr[0] * in1[stride_in * i]; + if (n_components > 2) + res2 += matrix_ptr[0] * in2[stride_in * i]; + } } else { const Number2 *matrix_ptr = matrix + nn_regular * n_columns; - res0 = matrix_ptr[0] * in[0]; + res0 = matrix_ptr[0] * in0[0]; + if (n_components > 1) + res1 = matrix_ptr[0] * in1[0]; + if (n_components > 2) + res2 = matrix_ptr[0] * in2[0]; for (int i = 1; i < mm; ++i) - res0 += matrix_ptr[i] * in[stride_in * i]; + { + res0 += matrix_ptr[i] * in0[stride_in * i]; + if (n_components > 1) + res1 += matrix_ptr[i] * in1[stride_in * i]; + if (n_components > 2) + res2 += matrix_ptr[i] * in2[stride_in * i]; + } } if (add) - out[0] += res0; + { + out0[0] += res0; + if (n_components > 1) + out1[0] += res1; + if (n_components > 2) + out2[0] += res2; + } else - out[0] = res0; + { + out0[0] = res0; + if (n_components > 1) + out1[0] = res1; + if (n_components > 2) + out2[0] = res2; + } } } } diff --git a/tests/matrix_free/tensor_product_kernels_01.cc b/tests/matrix_free/tensor_product_kernels_01.cc new file mode 100644 index 0000000000..a6427e6a6c --- /dev/null +++ b/tests/matrix_free/tensor_product_kernels_01.cc @@ -0,0 +1,134 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2020 - 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + + +// Tests the evaluate_general variant of the apply_matrix_vector_product +// function where either one, two or three vectors are used at the same time + +#include + +#include "../tests.h" + +template +void +test() +{ + const unsigned int n_columns_raw = 20; + const unsigned int n_rows_raw = 105; + + const unsigned int n_columns = transpose ? n_rows_raw : n_columns_raw; + const unsigned int n_rows = transpose ? n_columns_raw : n_rows_raw; + + double matrix[n_rows * n_columns]; + for (unsigned int i = 0; i < n_columns * n_rows; ++i) + matrix[i] = 0.1 * (i + 1); + + double in[n_columns]; + double in_batched_2[2 * n_columns]; + double in_batched[3 * n_columns]; + for (unsigned int i = 0; i < n_columns; ++i) + { + in[i] = i + 1; + in_batched[i] = i + 1; + in_batched[n_columns + i] = i + 1; + in_batched_2[i] = i + 1; + in_batched_2[n_columns + i] = i + 1; + in_batched[2 * n_columns + i] = i + 1; + } + + + // Reference solutions + double out[n_rows]; + double out_batched_2[2 * n_rows]; + double out_batched[3 * n_rows]; + + for (unsigned int i = 0; i < n_rows; ++i) + { + out[i] = i + 1; + out_batched[i] = i + 1; + out_batched[n_rows + i] = i + 1; + out_batched_2[i] = i + 1; + out_batched_2[n_rows + i] = i + 1; + out_batched[2 * n_rows + i] = i + 1; + } + + dealii::internal::apply_matrix_vector_product< + dealii::internal::EvaluatorVariant::evaluate_general, + dealii::internal::EvaluatorQuantity::value, + /*n_rows*/ n_rows_raw, + /*n_columns*/ n_columns_raw, + /*strid_in*/ 1, + /*strid_out*/ 1, + /*transpose_matrix*/ transpose, + /*add*/ add, + double, + double>(matrix, in, out); + + dealii::internal::apply_matrix_vector_product< + dealii::internal::EvaluatorVariant::evaluate_general, + dealii::internal::EvaluatorQuantity::value, + /*transpose_matrix*/ transpose, + /*add*/ add, + /*consider_strides*/ false, + double, + double, + /*n_components*/ 3>( + matrix, in_batched, out_batched, n_rows_raw, n_columns_raw, 1, 1); + + dealii::internal::apply_matrix_vector_product< + dealii::internal::EvaluatorVariant::evaluate_general, + dealii::internal::EvaluatorQuantity::value, + /*transpose_matrix*/ transpose, + /*add*/ add, + /*consider_strides*/ false, + double, + double, + /*n_components*/ 2>( + matrix, in_batched_2, out_batched_2, n_rows_raw, n_columns_raw, 1, 1); + + + // Check if outputs are the same + deallog << "Check results for transpose = " << transpose + << " and add = " << add << std::endl; + for (unsigned int i = 0; i < n_rows; ++i) + { + if ((out[i] != out_batched[i]) || (out[i] != out_batched[n_rows + i]) || + (out[i] != out_batched[2 * n_rows + i])) + deallog << "Error at entry " << i << std::endl; + else + deallog << "Correct line " << i << std::endl; + } + for (unsigned int i = 0; i < n_rows; ++i) + { + if ((out[i] != out_batched_2[i]) || (out[i] != out_batched_2[n_rows + i])) + deallog << "Error at entry " << i << std::endl; + else + deallog << "Correct line " << i << std::endl; + } + deallog << std::endl; +} + + + +int +main() +{ + initlog(); + deallog << std::setprecision(9); + test(); + test(); + test(); + test(); +} diff --git a/tests/matrix_free/tensor_product_kernels_01.output b/tests/matrix_free/tensor_product_kernels_01.output new file mode 100644 index 0000000000..64bbbc0251 --- /dev/null +++ b/tests/matrix_free/tensor_product_kernels_01.output @@ -0,0 +1,509 @@ + +DEAL::Check results for transpose = 0 and add = 0 +DEAL::Correct line 0 +DEAL::Correct line 1 +DEAL::Correct line 2 +DEAL::Correct line 3 +DEAL::Correct line 4 +DEAL::Correct line 5 +DEAL::Correct line 6 +DEAL::Correct line 7 +DEAL::Correct line 8 +DEAL::Correct line 9 +DEAL::Correct line 10 +DEAL::Correct line 11 +DEAL::Correct line 12 +DEAL::Correct line 13 +DEAL::Correct line 14 +DEAL::Correct line 15 +DEAL::Correct line 16 +DEAL::Correct line 17 +DEAL::Correct line 18 +DEAL::Correct line 19 +DEAL::Correct line 20 +DEAL::Correct line 21 +DEAL::Correct line 22 +DEAL::Correct line 23 +DEAL::Correct line 24 +DEAL::Correct line 25 +DEAL::Correct line 26 +DEAL::Correct line 27 +DEAL::Correct line 28 +DEAL::Correct line 29 +DEAL::Correct line 30 +DEAL::Correct line 31 +DEAL::Correct line 32 +DEAL::Correct line 33 +DEAL::Correct line 34 +DEAL::Correct line 35 +DEAL::Correct line 36 +DEAL::Correct line 37 +DEAL::Correct line 38 +DEAL::Correct line 39 +DEAL::Correct line 40 +DEAL::Correct line 41 +DEAL::Correct line 42 +DEAL::Correct line 43 +DEAL::Correct line 44 +DEAL::Correct line 45 +DEAL::Correct line 46 +DEAL::Correct line 47 +DEAL::Correct line 48 +DEAL::Correct line 49 +DEAL::Correct line 50 +DEAL::Correct line 51 +DEAL::Correct line 52 +DEAL::Correct line 53 +DEAL::Correct line 54 +DEAL::Correct line 55 +DEAL::Correct line 56 +DEAL::Correct line 57 +DEAL::Correct line 58 +DEAL::Correct line 59 +DEAL::Correct line 60 +DEAL::Correct line 61 +DEAL::Correct line 62 +DEAL::Correct line 63 +DEAL::Correct line 64 +DEAL::Correct line 65 +DEAL::Correct line 66 +DEAL::Correct line 67 +DEAL::Correct line 68 +DEAL::Correct line 69 +DEAL::Correct line 70 +DEAL::Correct line 71 +DEAL::Correct line 72 +DEAL::Correct line 73 +DEAL::Correct line 74 +DEAL::Correct line 75 +DEAL::Correct line 76 +DEAL::Correct line 77 +DEAL::Correct line 78 +DEAL::Correct line 79 +DEAL::Correct line 80 +DEAL::Correct line 81 +DEAL::Correct line 82 +DEAL::Correct line 83 +DEAL::Correct line 84 +DEAL::Correct line 85 +DEAL::Correct line 86 +DEAL::Correct line 87 +DEAL::Correct line 88 +DEAL::Correct line 89 +DEAL::Correct line 90 +DEAL::Correct line 91 +DEAL::Correct line 92 +DEAL::Correct line 93 +DEAL::Correct line 94 +DEAL::Correct line 95 +DEAL::Correct line 96 +DEAL::Correct line 97 +DEAL::Correct line 98 +DEAL::Correct line 99 +DEAL::Correct line 100 +DEAL::Correct line 101 +DEAL::Correct line 102 +DEAL::Correct line 103 +DEAL::Correct line 104 +DEAL::Correct line 0 +DEAL::Correct line 1 +DEAL::Correct line 2 +DEAL::Correct line 3 +DEAL::Correct line 4 +DEAL::Correct line 5 +DEAL::Correct line 6 +DEAL::Correct line 7 +DEAL::Correct line 8 +DEAL::Correct line 9 +DEAL::Correct line 10 +DEAL::Correct line 11 +DEAL::Correct line 12 +DEAL::Correct line 13 +DEAL::Correct line 14 +DEAL::Correct line 15 +DEAL::Correct line 16 +DEAL::Correct line 17 +DEAL::Correct line 18 +DEAL::Correct line 19 +DEAL::Correct line 20 +DEAL::Correct line 21 +DEAL::Correct line 22 +DEAL::Correct line 23 +DEAL::Correct line 24 +DEAL::Correct line 25 +DEAL::Correct line 26 +DEAL::Correct line 27 +DEAL::Correct line 28 +DEAL::Correct line 29 +DEAL::Correct line 30 +DEAL::Correct line 31 +DEAL::Correct line 32 +DEAL::Correct line 33 +DEAL::Correct line 34 +DEAL::Correct line 35 +DEAL::Correct line 36 +DEAL::Correct line 37 +DEAL::Correct line 38 +DEAL::Correct line 39 +DEAL::Correct line 40 +DEAL::Correct line 41 +DEAL::Correct line 42 +DEAL::Correct line 43 +DEAL::Correct line 44 +DEAL::Correct line 45 +DEAL::Correct line 46 +DEAL::Correct line 47 +DEAL::Correct line 48 +DEAL::Correct line 49 +DEAL::Correct line 50 +DEAL::Correct line 51 +DEAL::Correct line 52 +DEAL::Correct line 53 +DEAL::Correct line 54 +DEAL::Correct line 55 +DEAL::Correct line 56 +DEAL::Correct line 57 +DEAL::Correct line 58 +DEAL::Correct line 59 +DEAL::Correct line 60 +DEAL::Correct line 61 +DEAL::Correct line 62 +DEAL::Correct line 63 +DEAL::Correct line 64 +DEAL::Correct line 65 +DEAL::Correct line 66 +DEAL::Correct line 67 +DEAL::Correct line 68 +DEAL::Correct line 69 +DEAL::Correct line 70 +DEAL::Correct line 71 +DEAL::Correct line 72 +DEAL::Correct line 73 +DEAL::Correct line 74 +DEAL::Correct line 75 +DEAL::Correct line 76 +DEAL::Correct line 77 +DEAL::Correct line 78 +DEAL::Correct line 79 +DEAL::Correct line 80 +DEAL::Correct line 81 +DEAL::Correct line 82 +DEAL::Correct line 83 +DEAL::Correct line 84 +DEAL::Correct line 85 +DEAL::Correct line 86 +DEAL::Correct line 87 +DEAL::Correct line 88 +DEAL::Correct line 89 +DEAL::Correct line 90 +DEAL::Correct line 91 +DEAL::Correct line 92 +DEAL::Correct line 93 +DEAL::Correct line 94 +DEAL::Correct line 95 +DEAL::Correct line 96 +DEAL::Correct line 97 +DEAL::Correct line 98 +DEAL::Correct line 99 +DEAL::Correct line 100 +DEAL::Correct line 101 +DEAL::Correct line 102 +DEAL::Correct line 103 +DEAL::Correct line 104 +DEAL:: +DEAL::Check results for transpose = 1 and add = 0 +DEAL::Correct line 0 +DEAL::Correct line 1 +DEAL::Correct line 2 +DEAL::Correct line 3 +DEAL::Correct line 4 +DEAL::Correct line 5 +DEAL::Correct line 6 +DEAL::Correct line 7 +DEAL::Correct line 8 +DEAL::Correct line 9 +DEAL::Correct line 10 +DEAL::Correct line 11 +DEAL::Correct line 12 +DEAL::Correct line 13 +DEAL::Correct line 14 +DEAL::Correct line 15 +DEAL::Correct line 16 +DEAL::Correct line 17 +DEAL::Correct line 18 +DEAL::Correct line 19 +DEAL::Correct line 0 +DEAL::Correct line 1 +DEAL::Correct line 2 +DEAL::Correct line 3 +DEAL::Correct line 4 +DEAL::Correct line 5 +DEAL::Correct line 6 +DEAL::Correct line 7 +DEAL::Correct line 8 +DEAL::Correct line 9 +DEAL::Correct line 10 +DEAL::Correct line 11 +DEAL::Correct line 12 +DEAL::Correct line 13 +DEAL::Correct line 14 +DEAL::Correct line 15 +DEAL::Correct line 16 +DEAL::Correct line 17 +DEAL::Correct line 18 +DEAL::Correct line 19 +DEAL:: +DEAL::Check results for transpose = 0 and add = 1 +DEAL::Correct line 0 +DEAL::Correct line 1 +DEAL::Correct line 2 +DEAL::Correct line 3 +DEAL::Correct line 4 +DEAL::Correct line 5 +DEAL::Correct line 6 +DEAL::Correct line 7 +DEAL::Correct line 8 +DEAL::Correct line 9 +DEAL::Correct line 10 +DEAL::Correct line 11 +DEAL::Correct line 12 +DEAL::Correct line 13 +DEAL::Correct line 14 +DEAL::Correct line 15 +DEAL::Correct line 16 +DEAL::Correct line 17 +DEAL::Correct line 18 +DEAL::Correct line 19 +DEAL::Correct line 20 +DEAL::Correct line 21 +DEAL::Correct line 22 +DEAL::Correct line 23 +DEAL::Correct line 24 +DEAL::Correct line 25 +DEAL::Correct line 26 +DEAL::Correct line 27 +DEAL::Correct line 28 +DEAL::Correct line 29 +DEAL::Correct line 30 +DEAL::Correct line 31 +DEAL::Correct line 32 +DEAL::Correct line 33 +DEAL::Correct line 34 +DEAL::Correct line 35 +DEAL::Correct line 36 +DEAL::Correct line 37 +DEAL::Correct line 38 +DEAL::Correct line 39 +DEAL::Correct line 40 +DEAL::Correct line 41 +DEAL::Correct line 42 +DEAL::Correct line 43 +DEAL::Correct line 44 +DEAL::Correct line 45 +DEAL::Correct line 46 +DEAL::Correct line 47 +DEAL::Correct line 48 +DEAL::Correct line 49 +DEAL::Correct line 50 +DEAL::Correct line 51 +DEAL::Correct line 52 +DEAL::Correct line 53 +DEAL::Correct line 54 +DEAL::Correct line 55 +DEAL::Correct line 56 +DEAL::Correct line 57 +DEAL::Correct line 58 +DEAL::Correct line 59 +DEAL::Correct line 60 +DEAL::Correct line 61 +DEAL::Correct line 62 +DEAL::Correct line 63 +DEAL::Correct line 64 +DEAL::Correct line 65 +DEAL::Correct line 66 +DEAL::Correct line 67 +DEAL::Correct line 68 +DEAL::Correct line 69 +DEAL::Correct line 70 +DEAL::Correct line 71 +DEAL::Correct line 72 +DEAL::Correct line 73 +DEAL::Correct line 74 +DEAL::Correct line 75 +DEAL::Correct line 76 +DEAL::Correct line 77 +DEAL::Correct line 78 +DEAL::Correct line 79 +DEAL::Correct line 80 +DEAL::Correct line 81 +DEAL::Correct line 82 +DEAL::Correct line 83 +DEAL::Correct line 84 +DEAL::Correct line 85 +DEAL::Correct line 86 +DEAL::Correct line 87 +DEAL::Correct line 88 +DEAL::Correct line 89 +DEAL::Correct line 90 +DEAL::Correct line 91 +DEAL::Correct line 92 +DEAL::Correct line 93 +DEAL::Correct line 94 +DEAL::Correct line 95 +DEAL::Correct line 96 +DEAL::Correct line 97 +DEAL::Correct line 98 +DEAL::Correct line 99 +DEAL::Correct line 100 +DEAL::Correct line 101 +DEAL::Correct line 102 +DEAL::Correct line 103 +DEAL::Correct line 104 +DEAL::Correct line 0 +DEAL::Correct line 1 +DEAL::Correct line 2 +DEAL::Correct line 3 +DEAL::Correct line 4 +DEAL::Correct line 5 +DEAL::Correct line 6 +DEAL::Correct line 7 +DEAL::Correct line 8 +DEAL::Correct line 9 +DEAL::Correct line 10 +DEAL::Correct line 11 +DEAL::Correct line 12 +DEAL::Correct line 13 +DEAL::Correct line 14 +DEAL::Correct line 15 +DEAL::Correct line 16 +DEAL::Correct line 17 +DEAL::Correct line 18 +DEAL::Correct line 19 +DEAL::Correct line 20 +DEAL::Correct line 21 +DEAL::Correct line 22 +DEAL::Correct line 23 +DEAL::Correct line 24 +DEAL::Correct line 25 +DEAL::Correct line 26 +DEAL::Correct line 27 +DEAL::Correct line 28 +DEAL::Correct line 29 +DEAL::Correct line 30 +DEAL::Correct line 31 +DEAL::Correct line 32 +DEAL::Correct line 33 +DEAL::Correct line 34 +DEAL::Correct line 35 +DEAL::Correct line 36 +DEAL::Correct line 37 +DEAL::Correct line 38 +DEAL::Correct line 39 +DEAL::Correct line 40 +DEAL::Correct line 41 +DEAL::Correct line 42 +DEAL::Correct line 43 +DEAL::Correct line 44 +DEAL::Correct line 45 +DEAL::Correct line 46 +DEAL::Correct line 47 +DEAL::Correct line 48 +DEAL::Correct line 49 +DEAL::Correct line 50 +DEAL::Correct line 51 +DEAL::Correct line 52 +DEAL::Correct line 53 +DEAL::Correct line 54 +DEAL::Correct line 55 +DEAL::Correct line 56 +DEAL::Correct line 57 +DEAL::Correct line 58 +DEAL::Correct line 59 +DEAL::Correct line 60 +DEAL::Correct line 61 +DEAL::Correct line 62 +DEAL::Correct line 63 +DEAL::Correct line 64 +DEAL::Correct line 65 +DEAL::Correct line 66 +DEAL::Correct line 67 +DEAL::Correct line 68 +DEAL::Correct line 69 +DEAL::Correct line 70 +DEAL::Correct line 71 +DEAL::Correct line 72 +DEAL::Correct line 73 +DEAL::Correct line 74 +DEAL::Correct line 75 +DEAL::Correct line 76 +DEAL::Correct line 77 +DEAL::Correct line 78 +DEAL::Correct line 79 +DEAL::Correct line 80 +DEAL::Correct line 81 +DEAL::Correct line 82 +DEAL::Correct line 83 +DEAL::Correct line 84 +DEAL::Correct line 85 +DEAL::Correct line 86 +DEAL::Correct line 87 +DEAL::Correct line 88 +DEAL::Correct line 89 +DEAL::Correct line 90 +DEAL::Correct line 91 +DEAL::Correct line 92 +DEAL::Correct line 93 +DEAL::Correct line 94 +DEAL::Correct line 95 +DEAL::Correct line 96 +DEAL::Correct line 97 +DEAL::Correct line 98 +DEAL::Correct line 99 +DEAL::Correct line 100 +DEAL::Correct line 101 +DEAL::Correct line 102 +DEAL::Correct line 103 +DEAL::Correct line 104 +DEAL:: +DEAL::Check results for transpose = 1 and add = 1 +DEAL::Correct line 0 +DEAL::Correct line 1 +DEAL::Correct line 2 +DEAL::Correct line 3 +DEAL::Correct line 4 +DEAL::Correct line 5 +DEAL::Correct line 6 +DEAL::Correct line 7 +DEAL::Correct line 8 +DEAL::Correct line 9 +DEAL::Correct line 10 +DEAL::Correct line 11 +DEAL::Correct line 12 +DEAL::Correct line 13 +DEAL::Correct line 14 +DEAL::Correct line 15 +DEAL::Correct line 16 +DEAL::Correct line 17 +DEAL::Correct line 18 +DEAL::Correct line 19 +DEAL::Correct line 0 +DEAL::Correct line 1 +DEAL::Correct line 2 +DEAL::Correct line 3 +DEAL::Correct line 4 +DEAL::Correct line 5 +DEAL::Correct line 6 +DEAL::Correct line 7 +DEAL::Correct line 8 +DEAL::Correct line 9 +DEAL::Correct line 10 +DEAL::Correct line 11 +DEAL::Correct line 12 +DEAL::Correct line 13 +DEAL::Correct line 14 +DEAL::Correct line 15 +DEAL::Correct line 16 +DEAL::Correct line 17 +DEAL::Correct line 18 +DEAL::Correct line 19 +DEAL:: diff --git a/tests/simplex/matrix_free_05.cc b/tests/simplex/matrix_free_05.cc new file mode 100644 index 0000000000..8af025dcde --- /dev/null +++ b/tests/simplex/matrix_free_05.cc @@ -0,0 +1,352 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2020 - 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + +// Solve Poisson problem problem on a simplex mesh with DG and MatrixFree and +// dim components. + +#include + +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include "../tests.h" + +#include "./simplex_grids.h" + + + +double +get_penalty_parameter(const unsigned int i, + const unsigned int j, + const unsigned int degree) +{ + if (degree == 1) + { + if (i != j) + return 8.0; + if (i == 0) + return 8.0; + if (i == 1) + return 64.0; + } + else if (degree == 2) + { + if (i != j) + return 32.0; + if (i == 0) + return 32.0; + if (i == 1) + return 64.0; + } + + DEAL_II_NOT_IMPLEMENTED(); + + return 0.0; +} + + + +template +class PoissonOperator +{ +public: + using VectorType = LinearAlgebra::distributed::Vector; + using number = double; + + using FECellIntegrator = FEEvaluation; + using FEFaceIntegrator = FEFaceEvaluation; + + PoissonOperator(const MatrixFree &matrix_free, + const unsigned int degree) + : matrix_free(matrix_free) + , degree(degree) + {} + + void + initialize_dof_vector(VectorType &vec) + { + matrix_free.initialize_dof_vector(vec); + } + + void + rhs(VectorType &vec) const + { + const int dummy = 0; + + matrix_free.template cell_loop( + [&](const auto &data, auto &dst, const auto &, const auto range) { + FECellIntegrator phi(matrix_free, range); + Tensor<1, dim> value; + for (unsigned int d = 0; d < dim; d++) + value[d] = 0.0; + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + phi.reinit(cell); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_value(value, q); + + phi.integrate_scatter(EvaluationFlags::values, dst); + } + }, + vec, + dummy, + true); + } + + void + vmult(VectorType &dst, const VectorType &src) const + { + matrix_free.template loop( + [&](const auto &data, auto &dst, const auto &src, const auto range) { + FECellIntegrator phi(matrix_free, range); + + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + phi.reinit(cell); + phi.gather_evaluate(src, EvaluationFlags::gradients); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(phi.get_gradient(q), q); + phi.integrate_scatter(EvaluationFlags::gradients, dst); + } + }, + [&](const auto &data, auto &dst, const auto &src, const auto range) { + FEFaceIntegrator fe_eval(data, range, true); + FEFaceIntegrator fe_eval_neighbor(data, range, false); + + for (unsigned int face = range.first; face < range.second; ++face) + { + fe_eval.reinit(face); + fe_eval_neighbor.reinit(face); + + fe_eval.gather_evaluate(src, + EvaluationFlags::values | + EvaluationFlags::gradients); + fe_eval_neighbor.gather_evaluate(src, + EvaluationFlags::values | + EvaluationFlags::gradients); + VectorizedArray sigmaF = + get_penalty_parameter(data.get_face_active_fe_index(range, true), + data.get_face_active_fe_index(range, false), + degree); + + for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) + { + Tensor<1, dim, VectorizedArray> average_value = + (fe_eval.get_value(q) - fe_eval_neighbor.get_value(q)) * 0.5; + Tensor<1, dim, VectorizedArray> average_valgrad = + fe_eval.get_normal_derivative(q) + + fe_eval_neighbor.get_normal_derivative(q); + average_valgrad = + average_value * 2. * sigmaF - average_valgrad * 0.5; + fe_eval.submit_normal_derivative(-average_value, q); + fe_eval_neighbor.submit_normal_derivative(-average_value, q); + fe_eval.submit_value(average_valgrad, q); + fe_eval_neighbor.submit_value(-average_valgrad, q); + } + fe_eval.integrate_scatter(EvaluationFlags::values | + EvaluationFlags::gradients, + dst); + fe_eval_neighbor.integrate_scatter(EvaluationFlags::values | + EvaluationFlags::gradients, + dst); + } + }, + [&](const auto &data, auto &dst, const auto &src, const auto range) { + FEFaceIntegrator fe_eval(data, range, true); + + for (unsigned int face = range.first; face < range.second; ++face) + { + fe_eval.reinit(face); + fe_eval.gather_evaluate(src, + EvaluationFlags::values | + EvaluationFlags::gradients); + VectorizedArray sigmaF = + get_penalty_parameter(data.get_face_active_fe_index(range), + data.get_face_active_fe_index(range), + degree); + + for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) + { + Tensor<1, dim, VectorizedArray> average_value = + fe_eval.get_value(q); + Tensor<1, dim, VectorizedArray> average_valgrad = + -fe_eval.get_normal_derivative(q); + average_valgrad += average_value * sigmaF; + fe_eval.submit_normal_derivative(-average_value, q); + fe_eval.submit_value(average_valgrad, q); + } + + fe_eval.integrate_scatter(EvaluationFlags::values | + EvaluationFlags::gradients, + dst); + } + }, + dst, + src, + true); + } + +private: + const MatrixFree &matrix_free; + const unsigned int degree; +}; + +template +void +test(const unsigned int degree) +{ + Triangulation tria; + + unsigned int subdivisions = degree == 1 ? 16 : 8; + + + GridGenerator::subdivided_hyper_cube_with_simplices(tria, subdivisions); + + FESystem fe(FE_SimplexDGP(degree), dim); + QGaussSimplex quad(degree + 1); + MappingFE mapping(FE_SimplexP(1)); + + + DoFHandler dof_handler(tria); + + dof_handler.distribute_dofs(fe); + + AffineConstraints constraints; + constraints.close(); + + const auto solve_and_postprocess = + [&](const auto &poisson_operator, + auto &x, + auto &b) -> std::pair { + ReductionControl reduction_control(3000, 1e-7, 1e-2); + SolverCG> solver(reduction_control); + + solver.solve(poisson_operator, x, b, PreconditionIdentity()); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + printf("Solved in %d iterations.\n", reduction_control.last_step()); + + constraints.distribute(x); + +#if 0 + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + x.update_ghost_values(); + data_out.add_data_vector(dof_handler, x, "solution"); + data_out.build_patches(mapping, degree); + data_out.write_vtu_with_pvtu_record("./", + "result-" + std::to_string(dim) + "-" + + std::to_string(degree) + "-" + + std::to_string(version), + 0, + MPI_COMM_WORLD); +#endif + + Vector difference(tria.n_active_cells()); + + deallog << "dim=" << dim << ' '; + deallog << "degree=" << degree << ' '; + + VectorTools::integrate_difference(mapping, + dof_handler, + x, + Functions::ZeroFunction(dim), + difference, + quad, + VectorTools::L2_norm); + + const double error = + VectorTools::compute_global_error(tria, difference, VectorTools::L2_norm); + + if (error < 0.042) + deallog << "OK" << std::endl; + else + deallog << "FAIL" << std::endl; + + return {reduction_control.last_step(), reduction_control.last_value()}; + }; + + const auto mf_algo = [&]() { + typename MatrixFree::AdditionalData additional_data; + additional_data.mapping_update_flags = update_gradients | update_values; + additional_data.mapping_update_flags_inner_faces = + update_gradients | update_values; + additional_data.mapping_update_flags_boundary_faces = + update_gradients | update_values; + additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + + MatrixFree matrix_free; + matrix_free.reinit( + mapping, dof_handler, constraints, quad, additional_data); + + PoissonOperator poisson_operator(matrix_free, degree); + + LinearAlgebra::distributed::Vector x, b; + poisson_operator.initialize_dof_vector(x); + poisson_operator.initialize_dof_vector(b); + + poisson_operator.rhs(b); + + return solve_and_postprocess(poisson_operator, x, b); + }; + + mf_algo(); +} + + +int +main(int argc, char **argv) +{ + initlog(); + + deallog.depth_file(1); + + Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); + test<2>(1); + test<3>(1); +} diff --git a/tests/simplex/matrix_free_05.output b/tests/simplex/matrix_free_05.output new file mode 100644 index 0000000000..9c0cc1041c --- /dev/null +++ b/tests/simplex/matrix_free_05.output @@ -0,0 +1,3 @@ + +DEAL::dim=2 degree=1 OK +DEAL::dim=3 degree=1 OK -- 2.39.5