Merge pull request #628 from tigercosmos/neon_simd_1113

add the memory alignment check for NEON SIMD

Author: yungyuc

cpp/modmesh/simd/neon/neon.hpp

@@ -59,6 +59,11 @@ const T * check_between(T const * start, T const * end, T const & min_val, T con
     using cmpvec_t = type::vector_t<uint64_t>;
     constexpr size_t N_lane = type::vector_lane<T>;
 
+#ifndef NDEBUG
+    constexpr size_t alignment = get_recommended_alignment();
+    detail::check_alignment(start, alignment, "check_between start");
+#endif
+
     vec_t max_vec = vdupq(max_val);
     vec_t min_vec = vdupq(min_val);
     vec_t data_vec = {};
@@ -117,6 +122,14 @@ void add(T * dest, T const * dest_end, T const * src1, T const * src2)
     {
         using vec_t = type::vector_t<T>;
         constexpr size_t N_lane = type::vector_lane<T>;
+
+#ifndef NDEBUG
+        constexpr size_t alignment = get_recommended_alignment();
+        detail::check_alignment(dest, alignment, "add dest");
+        detail::check_alignment(src1, alignment, "add src1");
+        detail::check_alignment(src2, alignment, "add src2");
+#endif
+
         vec_t src1_vec;
         vec_t src2_vec;
         vec_t res_vec;
@@ -146,6 +159,14 @@ void sub(T * dest, T const * dest_end, T const * src1, T const * src2)
     {
         using vec_t = type::vector_t<T>;
         constexpr size_t N_lane = type::vector_lane<T>;
+
+#ifndef NDEBUG
+        constexpr size_t alignment = get_recommended_alignment();
+        detail::check_alignment(dest, alignment, "sub dest");
+        detail::check_alignment(src1, alignment, "sub src1");
+        detail::check_alignment(src2, alignment, "sub src2");
+#endif
+
         vec_t src1_vec;
         vec_t src2_vec;
         vec_t res_vec;
@@ -175,6 +196,14 @@ void mul(T * dest, T const * dest_end, T const * src1, T const * src2)
     {
         using vec_t = type::vector_t<T>;
         constexpr size_t N_lane = type::vector_lane<T>;
+
+#ifndef NDEBUG
+        constexpr size_t alignment = get_recommended_alignment();
+        detail::check_alignment(dest, alignment, "mul dest");
+        detail::check_alignment(src1, alignment, "mul src1");
+        detail::check_alignment(src2, alignment, "mul src2");
+#endif
+
         vec_t src1_vec;
         vec_t src2_vec;
         vec_t res_vec;
@@ -204,6 +233,14 @@ void div(T * dest, T const * dest_end, T const * src1, T const * src2)
     {
         using vec_t = type::vector_t<T>;
         constexpr size_t N_lane = type::vector_lane<T>;
+
+#ifndef NDEBUG
+        constexpr size_t alignment = get_recommended_alignment();
+        detail::check_alignment(dest, alignment, "div dest");
+        detail::check_alignment(src1, alignment, "div src1");
+        detail::check_alignment(src2, alignment, "div src2");
+#endif
+
         vec_t src1_vec;
         vec_t src2_vec;
         vec_t res_vec;

cpp/modmesh/simd/simd.hpp

@@ -28,8 +28,8 @@
  * POSSIBILITY OF SUCH DAMAGE.
  */
 
-#include <modmesh/simd/simd_support.hpp>
 #include <modmesh/simd/simd_generic.hpp>
+#include <modmesh/simd/simd_support.hpp>
 
 #include <modmesh/simd/neon/neon.hpp>
 
@@ -39,6 +39,48 @@ namespace modmesh
 namespace simd
 {
 
+namespace detail
+{
+#ifndef NDEBUG
+template <typename T>
+bool is_aligned(T const * pointer, size_t alignment)
+{
+    return (reinterpret_cast<std::uintptr_t>(pointer) % alignment) == 0;
+}
+
+template <typename T>
+void check_alignment(T const * pointer, size_t required_alignment, const char * name)
+{
+    if (!is_aligned(pointer, required_alignment))
+    {
+        std::fprintf(stderr,
+                     "Warning: %s pointer %p is not aligned to %zu bytes. "
+                     "SIMD performance may be degraded.\n",
+                     name,
+                     static_cast<const void *>(pointer),
+                     required_alignment);
+    }
+}
+#endif
+
+// Get the recommended memory alignment for SIMD operations based on the detected SIMD instruction set.
+inline constexpr size_t get_recommended_alignment()
+{
+#if defined(__aarch64__) || defined(__arm__)
+    return 16;
+#elif defined(__AVX512F__)
+    return 64;
+#elif defined(__AVX__) || defined(__AVX2__)
+    return 32;
+#elif defined(__SSE__) || defined(__SSE2__) || defined(__SSE3__) || defined(__SSSE3__) || defined(__SSE4_1__) || defined(__SSE4_2__)
+    return 16;
+#else
+    return 0;
+#endif
+}
+
+} // namespace detail
+
 // Check if each element from start to end (excluded end) is within the range [min_val, max_val)
 template <typename T>
 const T * check_between(T const * start, T const * end, T const & min_val, T const & max_val)

tests/test_buffer.py

@@ -1220,6 +1220,249 @@ def test_alignment_with_different_shapes(self):
         self.assertEqual(64, array3d.alignment)
         self.assertEqual((2, 4, 4), array3d.shape)
 
+    def test_alignment_with_simd_operations(self):
+        size = 16
+        alignments = [16, 32, 64]
+
+        for alignment in alignments:
+            arr1 = modmesh.SimpleArrayFloat64((size,), alignment)
+            arr2 = modmesh.SimpleArrayFloat64((size,), alignment)
+
+            self.assertEqual(alignment, arr1.alignment)
+            self.assertEqual(alignment, arr2.alignment)
+
+            for index in range(size):
+                arr1[index] = index * 2.0
+                arr2[index] = index * 3.0
+
+            result_add = arr1.add_simd(arr2)
+            self.assertEqual(0, result_add.alignment)
+            for index in range(size):
+                expected = index * 2.0 + index * 3.0
+                self.assertAlmostEqual(expected, result_add[index])
+
+            result_sub = arr1.sub_simd(arr2)
+            self.assertEqual(0, result_sub.alignment)
+            for index in range(size):
+                expected = index * 2.0 - index * 3.0
+                self.assertAlmostEqual(expected, result_sub[index])
+
+            result_mul = arr1.mul_simd(arr2)
+            self.assertEqual(0, result_mul.alignment)
+            for index in range(size):
+                expected = index * 2.0 * index * 3.0
+                self.assertAlmostEqual(expected, result_mul[index])
+
+            result_div = arr2.div_simd(arr1)
+            self.assertEqual(0, result_div.alignment)
+            for index in range(1, size):
+                expected = index * 3.0 / (index * 2.0)
+                self.assertAlmostEqual(expected, result_div[index])
+
+    def test_alignment_with_simd_operations_multidimensional(self):
+        alignments = [16, 32, 64]
+
+        for alignment in alignments:
+            array1_2d = modmesh.SimpleArrayFloat64((4, 8), alignment)
+            array2_2d = modmesh.SimpleArrayFloat64((4, 8), alignment)
+
+            self.assertEqual(alignment, array1_2d.alignment)
+            self.assertEqual(alignment, array2_2d.alignment)
+            self.assertEqual((4, 8), array1_2d.shape)
+            self.assertEqual((4, 8), array2_2d.shape)
+
+            for i in range(4):
+                for j in range(8):
+                    array1_2d[i, j] = (i * 8 + j) * 2.0
+                    array2_2d[i, j] = (i * 8 + j) * 3.0
+
+            result_add_2d = array1_2d.add_simd(array2_2d)
+            self.assertEqual(0, result_add_2d.alignment)
+            self.assertEqual((4, 8), result_add_2d.shape)
+            for i in range(4):
+                for j in range(8):
+                    value = i * 8 + j
+                    expected = value * 2.0 + value * 3.0
+                    self.assertAlmostEqual(expected, result_add_2d[i, j])
+
+            result_sub_2d = array1_2d.sub_simd(array2_2d)
+            self.assertEqual(0, result_sub_2d.alignment)
+            self.assertEqual((4, 8), result_sub_2d.shape)
+            for i in range(4):
+                for j in range(8):
+                    value = i * 8 + j
+                    expected = value * 2.0 - value * 3.0
+                    self.assertAlmostEqual(expected, result_sub_2d[i, j])
+
+            result_mul_2d = array1_2d.mul_simd(array2_2d)
+            self.assertEqual(0, result_mul_2d.alignment)
+            self.assertEqual((4, 8), result_mul_2d.shape)
+            for i in range(4):
+                for j in range(8):
+                    value = i * 8 + j
+                    expected = value * 2.0 * value * 3.0
+                    self.assertAlmostEqual(expected, result_mul_2d[i, j])
+
+            result_div_2d = array2_2d.div_simd(array1_2d)
+            self.assertEqual(0, result_div_2d.alignment)
+            self.assertEqual((4, 8), result_div_2d.shape)
+            for i in range(4):
+                for j in range(8):
+                    value = i * 8 + j
+                    if value > 0:
+                        expected = value * 3.0 / (value * 2.0)
+                        self.assertAlmostEqual(expected, result_div_2d[i, j])
+
+            array1_3d = modmesh.SimpleArrayFloat64((2, 4, 4), alignment)
+            array2_3d = modmesh.SimpleArrayFloat64((2, 4, 4), alignment)
+
+            self.assertEqual(alignment, array1_3d.alignment)
+            self.assertEqual(alignment, array2_3d.alignment)
+            self.assertEqual((2, 4, 4), array1_3d.shape)
+            self.assertEqual((2, 4, 4), array2_3d.shape)
+
+            for i in range(2):
+                for j in range(4):
+                    for k in range(4):
+                        array1_3d[i, j, k] = (i * 16 + j * 4 + k) * 2.0
+                        array2_3d[i, j, k] = (i * 16 + j * 4 + k) * 3.0
+
+            result_add_3d = array1_3d.add_simd(array2_3d)
+            self.assertEqual(0, result_add_3d.alignment)
+            self.assertEqual((2, 4, 4), result_add_3d.shape)
+            for i in range(2):
+                for j in range(4):
+                    for k in range(4):
+                        value = i * 16 + j * 4 + k
+                        expected = value * 2.0 + value * 3.0
+                        self.assertAlmostEqual(expected,
+                                               result_add_3d[i, j, k])
+
+            result_sub_3d = array1_3d.sub_simd(array2_3d)
+            self.assertEqual(0, result_sub_3d.alignment)
+            self.assertEqual((2, 4, 4), result_sub_3d.shape)
+            for i in range(2):
+                for j in range(4):
+                    for k in range(4):
+                        value = i * 16 + j * 4 + k
+                        expected = value * 2.0 - value * 3.0
+                        self.assertAlmostEqual(expected,
+                                               result_sub_3d[i, j, k])
+
+            result_mul_3d = array1_3d.mul_simd(array2_3d)
+            self.assertEqual(0, result_mul_3d.alignment)
+            self.assertEqual((2, 4, 4), result_mul_3d.shape)
+            for i in range(2):
+                for j in range(4):
+                    for k in range(4):
+                        value = i * 16 + j * 4 + k
+                        expected = value * 2.0 * value * 3.0
+                        self.assertAlmostEqual(expected,
+                                               result_mul_3d[i, j, k])
+
+            result_div_3d = array2_3d.div_simd(array1_3d)
+            self.assertEqual(0, result_div_3d.alignment)
+            self.assertEqual((2, 4, 4), result_div_3d.shape)
+            for i in range(2):
+                for j in range(4):
+                    for k in range(4):
+                        value = i * 16 + j * 4 + k
+                        if value > 0:
+                            expected = value * 3.0 / (value * 2.0)
+                            self.assertAlmostEqual(expected,
+                                                   result_div_3d[i, j, k])
+
+    def test_alignment_size_validation_multidimensional(self):
+        with self.assertRaisesRegex(
+                ValueError,
+                "ConcreteBuffer::allocate: size .* must be a multiple of alignment 16"  # noqa E501
+        ):
+            modmesh.SimpleArrayFloat64((1, 3), 16)
+
+        with self.assertRaisesRegex(
+                ValueError,
+                "ConcreteBuffer::allocate: size .* must be a multiple of alignment 32"  # noqa E501
+        ):
+            modmesh.SimpleArrayFloat64((1, 3), 32)
+
+        with self.assertRaisesRegex(
+                ValueError,
+                "ConcreteBuffer::allocate: size .* must be a multiple of alignment 64"  # noqa E501
+        ):
+            modmesh.SimpleArrayFloat64((3, 3), 64)
+
+        with self.assertRaisesRegex(
+                ValueError,
+                "ConcreteBuffer::allocate: size .* must be a multiple of alignment 16"  # noqa E501
+        ):
+            modmesh.SimpleArrayFloat64((1, 1, 1), 16)
+
+        with self.assertRaisesRegex(
+                ValueError,
+                "ConcreteBuffer::allocate: size .* must be a multiple of alignment 32"  # noqa E501
+        ):
+            modmesh.SimpleArrayFloat64((1, 1, 1), 32)
+
+        with self.assertRaisesRegex(
+                ValueError,
+                "ConcreteBuffer::allocate: size .* must be a multiple of alignment 64"  # noqa E501
+        ):
+            modmesh.SimpleArrayFloat64((2, 3, 5), 64)
+
+    def test_alignment_with_unaligned_rows(self):
+        # 2D arrays that row is not aligned
+        array1_2d = modmesh.SimpleArrayFloat64((2, 3), 16)
+        array2_2d = modmesh.SimpleArrayFloat64((2, 3), 16)
+
+        self.assertEqual(16, array1_2d.alignment)
+        self.assertEqual(16, array2_2d.alignment)
+        self.assertEqual((2, 3), array1_2d.shape)
+
+        for i in range(2):
+            for j in range(3):
+                array1_2d[i, j] = (i * 3 + j) * 2.0
+                array2_2d[i, j] = (i * 3 + j) * 3.0
+
+        # SIMD ops must tolerate the unaligned row stride and still yield exact math.  # noqa: E501
+        result_add = array1_2d.add_simd(array2_2d)
+        result_sub = array1_2d.sub_simd(array2_2d)
+        result_mul = array1_2d.mul_simd(array2_2d)
+
+        for i in range(2):
+            for j in range(3):
+                value = i * 3 + j
+                self.assertAlmostEqual(value * 5.0, result_add[i, j])
+                self.assertAlmostEqual(value * -1.0, result_sub[i, j])
+                self.assertAlmostEqual(value * value * 6.0, result_mul[i, j])
+
+        # Repeat with 3D data that the innermost dimension is unaligned.
+        array1_3d = modmesh.SimpleArrayFloat64((2, 2, 2), 32)
+        array2_3d = modmesh.SimpleArrayFloat64((2, 2, 2), 32)
+
+        self.assertEqual(32, array1_3d.alignment)
+        self.assertEqual(32, array2_3d.alignment)
+
+        for i in range(2):
+            for j in range(2):
+                for k in range(2):
+                    array1_3d[i, j, k] = (i * 4 + j * 2 + k) * 2.0
+                    array2_3d[i, j, k] = (i * 4 + j * 2 + k) * 3.0
+
+        result_add_3d = array1_3d.add_simd(array2_3d)
+        result_sub_3d = array1_3d.sub_simd(array2_3d)
+        result_mul_3d = array1_3d.mul_simd(array2_3d)
+
+        for i in range(2):
+            for j in range(2):
+                for k in range(2):
+                    value = i * 4 + j * 2 + k
+                    self.assertAlmostEqual(value * 5.0,
+                                           result_add_3d[i, j, k])
+                    self.assertAlmostEqual(value * -1.0,
+                                           result_sub_3d[i, j, k])
+                    self.assertAlmostEqual(value * value * 6.0,
+                                           result_mul_3d[i, j, k])
+
 
 class SimpleArrayCalculatorsTC(unittest.TestCase):
 
@@ -3154,12 +3397,6 @@ def test_alignment_with_push_back(self):
         # Verify alignment is maintained
         self.assertEqual(16, ct.alignment)
 
-    def test_alignment_with_simd_operations(self):
-        # TODO: implement tests for SIMD operations if applicable.
-        # It requires more effort on `as_array` and ConcreteBuffer side.
-        # See more: https://github.com/solvcon/modmesh/issues/620
-        pass
-
     def test_alignment_preserved_in_as_array(self):
         ct = modmesh.SimpleCollectorFloat64(16, 32)
         self.assertEqual(32, ct.alignment)