Vectorize NonMatching::MappingInfo #15137

bergbauer · 2023-04-25T14:39:03Z

This PR proposes to vectorize NonMatching::MappingInfo depending on the newly introduced template parameter Number. As this class interacts with FEPointEvaluation, which is vectorized internally also for scalar Number template arguments, the unit points are always stored vectorized to make them easily accessible.

To enable vectorized storage, separate AlignedVector objects are created for Jacobians, ... and a temporary internal::FEValuesImplementation::MappingRelatedData local object is used to use the interfaces which Mapping provides to compute mapping data.

The result is a 100% - 10% speedup in throughput of FEPointEvaluation depending on the polynomial degree and the dimension of the problem in my Laplace type benchmark.

peterrum · 2023-04-25T14:58:06Z

The result is a 100% - 10% speedup in throughput of FEPointEvaluation depending on the polynomial degree and the dimension of the problem in my Laplace type benchmark.

Very nice!

include/deal.II/matrix_free/fe_point_evaluation.h

include/deal.II/non_matching/mapping_info.h

kronbichler · 2023-04-25T15:59:39Z

include/deal.II/non_matching/mapping_info.h

     */
-    std::vector<unsigned int> unit_points_index;
+    AlignedVector<unsigned int> unit_points_index;


Does this need to be AlignedVector? In general, we've been trying to use AlignedVector only where necessary, i.e., large data structures (where we might profit from first-touch and similar) and for VectorizedArray fields. (The latter will be resolved once we require C++17, because then std::vector can respect alignment of types > 16 bytes, i.e., VectorizedArray.)

This is inspired by MappingInfoStorage, see

dealii/include/deal.II/matrix_free/mapping_info_storage.h

Line 299 in 4ec36db

AlignedVector<unsigned int> quadrature_point_offsets;

I can use std::vector if it aligns the memory properly for unsigned int.

include/deal.II/non_matching/mapping_info.h

kronbichler · 2023-04-25T16:31:38Z

/rebuild

peterrum · 2023-04-25T20:40:27Z

source/fe/mapping_q.cc

+          aux.resize(dim - 1);
+          aux[0].resize(n_original_q_points);
+          if (dim > 2)
+            aux[1].resize(n_original_q_points);


How is this better?

@jh66637 this?

I am not sure if this was the initial intention for this change, but with these changes, empty quadratures can be handled. I think in case of an empty quadrature rule a shortcut could be made, i.e., if (quadrature.size() == 0) return; in the first line of the function (or sth. similar).

@bergbauer @peterrum What do you think about this?

In my opinion the shortcut should not be in this function but in the reinit* functions of MappingInfo. Above change is necessary to resize the inner AlignedVectors correctly which is necessary if the InternalData object is reused.

include/deal.II/matrix_free/fe_point_evaluation.h

include/deal.II/non_matching/mapping_info.h

kronbichler

I think this goes to the right direction. I have some small comments left.

include/deal.II/base/vectorization.h

include/deal.II/matrix_free/fe_point_evaluation.h

kronbichler · 2023-04-28T15:47:41Z

include/deal.II/matrix_free/fe_point_evaluation.h

-            internal::FEPointEvaluation::
-              EvaluatorTypeTraits<dim, n_components, Number>::set_value(
-                val_and_grad.first, v, values[qb * stride + v]);
+          for (unsigned int v = 0; v < stride && q + v < n_q_points_scalar; ++v)


Is the compiler able to see that the second condition is a subset of the outer loop in line 1574 and hence removes the if statement altogether in vectorized mode with stride = 1?

include/deal.II/non_matching/mapping_info.h

peterrum

Looks good.

peterrum · 2023-05-01T08:28:13Z

include/deal.II/matrix_free/fe_point_evaluation.h

-          mapping_data.jacobians[point_index * n_lanes + v][d][s];
-  return vectorized_jacobian;
+  AssertIndexRange(point_index, n_q_points);
+  Assert(jacobian_ptr != nullptr, ExcMessage("jacobian_ptr is not set!"));


The message does not really help. I think in the case of FEEval the message tells you what to do:

dealii/include/deal.II/matrix_free/fe_evaluation.h

Lines 4522 to 4524 in aef0f70

Assert(this->jacobian != nullptr,

internal::ExcMatrixFreeAccessToUninitializedMappingField(

"update_gradients"));

Are you suggesting declaring an own exception with a verbose message what to do or to improve the current messages in place?

We have a verbose exception now which tells the user what to do!

peterrum · 2023-05-01T08:30:20Z

include/deal.II/matrix_free/fe_point_evaluation.h

+      read_value(const ScalarNumber vector_entry,
                 const unsigned int component,
-                 value_type &       result)
+                 scalar_value_type &result)


This looks odd. Different styles ScalarNumber vs. scalar_value_type. Could we write ScalarValueType. This is not really something what the user should use!?

This style dates back to the first implementation of this class and is therefore not the scope of this PR. I agree that it should be ScalarValueType. This should be done in a follow-up. The user should not use this but we expose it publicly with the using statements at the top of FEPointEvaluation.

peterrum · 2023-05-01T08:37:25Z

include/deal.II/non_matching/mapping_info.h

+        mapping_data.initialize(quadrature.size(), update_flags_mapping);
+
+        // reuse internal_mapping_data for MappingQ to avoid memory allocations
+        if (const MappingQ<dim, spacedim> *mapping_q =


Suggested change

if (const MappingQ<dim, spacedim> *mapping_q =

if (const MappingQ<dim, spacedim> * =

peterrum · 2023-05-02T21:02:10Z

Undefined symbols for architecture x86_64:
  "dealii::internal::VectorizedArrayTrait<dealii::VectorizedArray<double, 2ul> >::width", referenced from:
      dealii::internal::VectorizedArrayTrait<dealii::VectorizedArray<double, 2ul> >::get(dealii::VectorizedArray<double, 2ul>&, unsigned int) in data_out_resample.cc.o
  "dealii::internal::VectorizedArrayTrait<double>::width", referenced from:
      dealii::internal::VectorizedArrayTrait<double>::get(double&, unsigned int) in data_out_resample.cc.o

peterrum self-requested a review April 25, 2023 14:57

kronbichler approved these changes Apr 25, 2023

View reviewed changes

kronbichler added Matrix-free Mapping ready to test labels Apr 25, 2023

peterrum reviewed Apr 25, 2023

View reviewed changes

kronbichler reviewed Apr 28, 2023

View reviewed changes

peterrum reviewed May 1, 2023

View reviewed changes

bergbauer force-pushed the vectorized_fe_point_and_mapping_info branch from 9e42a0d to 0f1cad6 Compare May 2, 2023 16:07

kronbichler approved these changes May 2, 2023

View reviewed changes

peterrum approved these changes May 2, 2023

View reviewed changes

Vectorize MappingInfo

17a4b28

bergbauer force-pushed the vectorized_fe_point_and_mapping_info branch from 0f1cad6 to 17a4b28 Compare May 3, 2023 07:56

kronbichler merged commit e8d7c68 into dealii:master May 3, 2023
14 checks passed

This was referenced May 5, 2023

Avoid missing symbols and use constexpr functions instead #15181

Merged

Reduce overhead in calls to tensor product value function #15182

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Vectorize NonMatching::MappingInfo #15137

Vectorize NonMatching::MappingInfo #15137

bergbauer commented Apr 25, 2023

peterrum commented Apr 25, 2023

kronbichler Apr 25, 2023

bergbauer Apr 25, 2023

kronbichler commented Apr 25, 2023

peterrum Apr 25, 2023

bergbauer Apr 26, 2023

jh66637 Apr 26, 2023 •

edited

bergbauer Apr 27, 2023

kronbichler left a comment

kronbichler Apr 28, 2023

peterrum left a comment

peterrum May 1, 2023

bergbauer May 2, 2023

bergbauer May 2, 2023

peterrum May 1, 2023

bergbauer May 2, 2023

peterrum May 1, 2023

bergbauer May 2, 2023

peterrum commented May 2, 2023

	Assert(this->jacobian != nullptr,
	internal::ExcMatrixFreeAccessToUninitializedMappingField(
	"update_gradients"));

	if (const MappingQ<dim, spacedim> *mapping_q =
	if (const MappingQ<dim, spacedim> * =

Vectorize NonMatching::MappingInfo #15137

Vectorize NonMatching::MappingInfo #15137

Conversation

bergbauer commented Apr 25, 2023

peterrum commented Apr 25, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kronbichler commented Apr 25, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jh66637 Apr 26, 2023 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kronbichler left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

peterrum left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

peterrum commented May 2, 2023

jh66637 Apr 26, 2023 •

edited