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<evaluate_general,
- EvaluatorQuantity::value,
- /* transpose_matrix */ true,
- /* add */ false,
- /* consider_strides */ false>(
- 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<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ true,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 1>(
+ shape_values, in, out, n_dofs, n_q_points, 1, 1);
+ else if (c + 2 == n_components)
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ true,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 2>(
+ shape_values, in, out, n_dofs, n_q_points, 1, 1);
+ else
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ true,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 3>(
+ shape_values, in, out, n_dofs, n_q_points, 1, 1);
+
+ out += 3 * n_q_points;
+ in += 3 * n_dofs;
}
}
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<evaluate_general,
- EvaluatorQuantity::value,
- /* transpose_matrix */ true,
- /* add */ false,
- /* consider_strides */ false>(
- 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<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ true,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 1>(
+ shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+ else if (c + 2 == n_components)
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ true,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 2>(
+ shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+ else
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ true,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 3>(
+ shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+
+ out += 3 * n_q_points * dim;
+ in += 3 * n_dofs;
}
}
}
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<evaluate_general,
- EvaluatorQuantity::value,
- /* transpose_matrix */ false,
- /* add */ false,
- /* consider_strides */ false>(
- shape_values, in, out, n_dofs, n_q_points, 1, 1);
+ {
+ if (c + 1 == n_components)
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 1>(
+ shape_values, in, out, n_dofs, n_q_points, 1, 1);
+ else if (c + 2 == n_components)
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 2>(
+ shape_values, in, out, n_dofs, n_q_points, 1, 1);
+ else
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 3>(
+ shape_values, in, out, n_dofs, n_q_points, 1, 1);
+ }
else
- apply_matrix_vector_product<evaluate_general,
- EvaluatorQuantity::value,
- /* transpose_matrix */ false,
- /* add */ true,
- /* consider_strides */ false>(
- 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<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ true,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 1>(
+ shape_values, in, out, n_dofs, n_q_points, 1, 1);
+ else if (c + 2 == n_components)
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ true,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 2>(
+ shape_values, in, out, n_dofs, n_q_points, 1, 1);
+ else
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ true,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 3>(
+ shape_values, in, out, n_dofs, n_q_points, 1, 1);
+ }
+ out += 3 * n_dofs;
+ in += 3 * n_q_points;
}
}
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<evaluate_general,
- EvaluatorQuantity::value,
- /* transpose_matrix */ false,
- /* add */ false,
- /* consider_strides */ false>(
- shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+ {
+ if (c + 1 == n_components)
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 1>(
+ shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+ else if (c + 2 == n_components)
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 2>(
+ shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+ else
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 3>(
+ shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+ }
else
- apply_matrix_vector_product<evaluate_general,
- EvaluatorQuantity::value,
- /* transpose_matrix */ false,
- /* add */ true,
- /* consider_strides */ false>(
- 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<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ true,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 1>(
+ shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+ else if (c + 2 == n_components)
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ true,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 2>(
+ shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+ else
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ true,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 3>(
+ shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+ }
+ out += 3 * n_dofs;
+ in += 3 * n_q_points * dim;
}
}
}
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<evaluate_general,
- EvaluatorQuantity::value,
- /*transpose_matrix*/ true,
- /*add*/ false,
- /*consider_strides*/ false>(
- 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<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ true,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 1>(
+ shape_values, in, out, n_dofs, n_q_points, 1, 1);
+ else if (c + 2 == n_components)
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ true,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 2>(
+ shape_values, in, out, n_dofs, n_q_points, 1, 1);
+ else
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ true,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 3>(
+ shape_values, in, out, n_dofs, n_q_points, 1, 1);
+
+ out += 3 * n_q_points;
+ in += 3 * n_dofs;
}
}
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<evaluate_general,
- EvaluatorQuantity::value,
- /*transpose_matrix*/ true,
- /*add*/ false,
- /*consider_strides*/ false>(
- 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<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ true,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 1>(
+ shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+ else if (c + 2 == n_components)
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ true,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 2>(
+ shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+ else
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ true,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 3>(
+ shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+ out += 3 * n_q_points * dim;
+ in += 3 * n_dofs;
}
}
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<evaluate_general,
- EvaluatorQuantity::value,
- /*transpose_matrix*/ false,
- /*add*/ true,
- /*consider_strides*/ false>(
- shape_values, in, out, n_dofs, n_q_points, 1, 1);
+ {
+ if (c + 1 == n_components)
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ true,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 1>(
+ shape_values, in, out, n_dofs, n_q_points, 1, 1);
+ else if (c + 2 == n_components)
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ true,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 2>(
+ shape_values, in, out, n_dofs, n_q_points, 1, 1);
+ else
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ true,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 3>(
+ shape_values, in, out, n_dofs, n_q_points, 1, 1);
+ }
else
- apply_matrix_vector_product<evaluate_general,
- EvaluatorQuantity::value,
- /*transpose_matrix*/ false,
- /*add*/ false,
- /*consider_strides*/ false>(
- 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<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 1>(
+ shape_values, in, out, n_dofs, n_q_points, 1, 1);
+ else if (c + 2 == n_components)
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 2>(
+ shape_values, in, out, n_dofs, n_q_points, 1, 1);
+ else
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 3>(
+ shape_values, in, out, n_dofs, n_q_points, 1, 1);
+ }
+ in += 3 * n_q_points;
+ out += 3 * n_dofs;
}
}
{
if (!sum_into_values &&
!(integration_flag & EvaluationFlags::values))
- apply_matrix_vector_product<evaluate_general,
- EvaluatorQuantity::value,
- /*transpose_matrix*/ false,
- /*add*/ false,
- /*consider_strides*/ false>(
- shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+ {
+ if (c + 1 == n_components)
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 1>(
+ shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+ else if (c + 2 == n_components)
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 2>(
+ shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+ else
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ false,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 3>(
+ shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+ }
else
- apply_matrix_vector_product<evaluate_general,
- EvaluatorQuantity::value,
- /*transpose_matrix*/ false,
- /*add*/ true,
- /*consider_strides*/ false>(
- 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<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ true,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 1>(
+ shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+ else if (c + 2 == n_components)
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ true,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 2>(
+ shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+ else
+ apply_matrix_vector_product<evaluate_general,
+ EvaluatorQuantity::value,
+ /*transpose_matrix*/ false,
+ /*add*/ true,
+ /*consider_strides*/ false,
+ Number,
+ Number2,
+ /*n_components*/ 3>(
+ shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+ }
+ in += 3 * n_q_points * dim;
+ out += 3 * n_dofs;
}
}
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,
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];
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;
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"));
}
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
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
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
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;
+ }
}
}
}
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// 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 <deal.II/matrix_free/tensor_product_kernels.h>
+
+#include "../tests.h"
+
+template <bool transpose, bool add>
+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<false, false>();
+ test<true, false>();
+ test<false, true>();
+ test<true, true>();
+}
--- /dev/null
+
+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::
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// 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 <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_wedge_p.h>
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/sparse_matrix.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/meshworker/copy_data.h>
+#include <deal.II/meshworker/mesh_loop.h>
+#include <deal.II/meshworker/scratch_data.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#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 <int dim, int n_components = dim>
+class PoissonOperator
+{
+public:
+ using VectorType = LinearAlgebra::distributed::Vector<double>;
+ using number = double;
+
+ using FECellIntegrator = FEEvaluation<dim, -1, 0, n_components, number>;
+ using FEFaceIntegrator = FEFaceEvaluation<dim, -1, 0, n_components, number>;
+
+ PoissonOperator(const MatrixFree<dim, double> &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<VectorType, int>(
+ [&](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<VectorType, VectorType>(
+ [&](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<number> 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<number>> average_value =
+ (fe_eval.get_value(q) - fe_eval_neighbor.get_value(q)) * 0.5;
+ Tensor<1, dim, VectorizedArray<number>> 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<number> 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<number>> average_value =
+ fe_eval.get_value(q);
+ Tensor<1, dim, VectorizedArray<number>> 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<dim, double> &matrix_free;
+ const unsigned int degree;
+};
+
+template <int dim>
+void
+test(const unsigned int degree)
+{
+ Triangulation<dim> tria;
+
+ unsigned int subdivisions = degree == 1 ? 16 : 8;
+
+
+ GridGenerator::subdivided_hyper_cube_with_simplices(tria, subdivisions);
+
+ FESystem<dim> fe(FE_SimplexDGP<dim>(degree), dim);
+ QGaussSimplex<dim> quad(degree + 1);
+ MappingFE<dim> mapping(FE_SimplexP<dim>(1));
+
+
+ DoFHandler<dim> dof_handler(tria);
+
+ dof_handler.distribute_dofs(fe);
+
+ AffineConstraints<double> constraints;
+ constraints.close();
+
+ const auto solve_and_postprocess =
+ [&](const auto &poisson_operator,
+ auto &x,
+ auto &b) -> std::pair<unsigned int, double> {
+ ReductionControl reduction_control(3000, 1e-7, 1e-2);
+ SolverCG<std::remove_reference_t<decltype(x)>> 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<dim> 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<double> difference(tria.n_active_cells());
+
+ deallog << "dim=" << dim << ' ';
+ deallog << "degree=" << degree << ' ';
+
+ VectorTools::integrate_difference(mapping,
+ dof_handler,
+ x,
+ Functions::ZeroFunction<dim>(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<dim, double>::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<dim, double>::AdditionalData::none;
+
+ MatrixFree<dim, double> matrix_free;
+ matrix_free.reinit(
+ mapping, dof_handler, constraints, quad, additional_data);
+
+ PoissonOperator<dim> poisson_operator(matrix_free, degree);
+
+ LinearAlgebra::distributed::Vector<double> 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);
+}
--- /dev/null
+
+DEAL::dim=2 degree=1 OK
+DEAL::dim=3 degree=1 OK
projects / dealii.git / commitdiff