Expand multipole construction benchmarks, optimize for similar perfor…

…mance to manual initialization
JacksonCampolattaro · Oct 20, 2023 · a7a9ac1 · a7a9ac1
1 parent e4db894
commit a7a9ac1
Show file tree

Hide file tree

Showing 3 changed files with 95 additions and 23 deletions.
diff --git a/benchmarks/multipole.cpp b/benchmarks/multipole.cpp
@@ -29,35 +29,95 @@ TEST_CASE("benchmark: Multipole construction from a vector", "[Multipole]") {
 
         BENCHMARK("Octupole constructor") { return Octupole3f{vector}; };
         BENCHMARK("Octupole manual construction") {
-                                                        return std::tuple<
-                                                                SymmetricTensor3f<1>,
-                                                                SymmetricTensor3f<2>,
-                                                                SymmetricTensor3f<3>
-                                                        >{
-                                                                vector,
-                                                                SymmetricTensor3f<2>::CartesianPower(vector),
-                                                                SymmetricTensor3f<3>::CartesianPower(vector)
-                                                        };
-                                                    };
-    }
-
-    // Hexadecupole
-    {
-
-        BENCHMARK("Hexadecupole constructor") { return Hexadecupole3f{vector}; };
-        BENCHMARK("Hexadecupole manual construction") {
                                                       return std::tuple<
                                                               SymmetricTensor3f<1>,
                                                               SymmetricTensor3f<2>,
-                                                              SymmetricTensor3f<3>,
-                                                              SymmetricTensor3f<4>
+                                                              SymmetricTensor3f<3>
                                                       >{
                                                               vector,
                                                               SymmetricTensor3f<2>::CartesianPower(vector),
-                                                              SymmetricTensor3f<3>::CartesianPower(vector),
-                                                              SymmetricTensor3f<4>::CartesianPower(vector)
+                                                              SymmetricTensor3f<3>::CartesianPower(vector)
                                                       };
                                                   };
     }
 
+    // Hexadecupole
+    {
+
+        BENCHMARK("Hexadecupole constructor") { return Hexadecupole3f{vector}; };
+        BENCHMARK("Hexadecupole manual construction")
+                    {
+                        return std::tuple<
+                                SymmetricTensor3f<1>,
+                                SymmetricTensor3f<2>,
+                                SymmetricTensor3f<3>,
+                                SymmetricTensor3f<4>
+                        >{
+                                vector,
+                                SymmetricTensor3f<2>::CartesianPower(vector),
+                                SymmetricTensor3f<3>::CartesianPower(vector),
+                                SymmetricTensor3f<4>::CartesianPower(vector)
+                        };
+                    };
+        BENCHMARK("Hexadecupole cascading manual construction")
+                    {
+                        auto st2 = SymmetricTensor3f<2>::CartesianProduct(vector, vector);
+                        auto st3 = SymmetricTensor3f<3>::CartesianProduct(st2, vector);
+                        auto st4 = SymmetricTensor3f<4>::CartesianProduct(st3, vector);
+                        return std::tuple<
+                                SymmetricTensor3f<1>,
+                                SymmetricTensor3f<2>,
+                                SymmetricTensor3f<3>,
+                                SymmetricTensor3f<4>
+                        >{
+                                vector,
+                                st2,
+                                st3,
+                                st4
+                        };
+                    };
+    }
+
+    // Triacontadyupole
+    {
+
+        BENCHMARK("Triacontadyupole constructor") { return Triacontadyupole3f{vector}; };
+        BENCHMARK("Triacontadyupole manual construction")
+                    {
+                        return std::tuple<
+                                SymmetricTensor3f<1>,
+                                SymmetricTensor3f<2>,
+                                SymmetricTensor3f<3>,
+                                SymmetricTensor3f<4>,
+                                SymmetricTensor3f<5>
+                        >{
+                                vector,
+                                SymmetricTensor3f<2>::CartesianPower(vector),
+                                SymmetricTensor3f<3>::CartesianPower(vector),
+                                SymmetricTensor3f<4>::CartesianPower(vector),
+                                SymmetricTensor3f<5>::CartesianPower(vector)
+                        };
+                    };
+        BENCHMARK("Triacontadyupole cascading manual construction")
+                    {
+                        auto st2 = SymmetricTensor3f<2>::CartesianProduct(vector, vector);
+                        auto st3 = SymmetricTensor3f<3>::CartesianProduct(st2, vector);
+                        auto st4 = SymmetricTensor3f<4>::CartesianProduct(st3, vector);
+                        auto st5 = SymmetricTensor3f<5>::CartesianProduct(st4, vector);
+                        return std::tuple<
+                                SymmetricTensor3f<1>,
+                                SymmetricTensor3f<2>,
+                                SymmetricTensor3f<3>,
+                                SymmetricTensor3f<4>,
+                                SymmetricTensor3f<5>
+                        >{
+                                vector,
+                                st2,
+                                st3,
+                                st4,
+                                st5
+                        };
+                    };
+    }
+
 }
diff --git a/include/symtensor/Multipole.h b/include/symtensor/Multipole.h
@@ -65,6 +65,12 @@ namespace symtensor {
     /// 3d hexadecupole with floating point elements
     using Hexadecupole3f = Multipole3f<4>;
 
+    /// 2d triacontadyupole with floating point elements
+    using Triacontadyupole2f = Multipole2f<5>;
+
+    /// 3d triacontadyupole with floating point elements
+    using Triacontadyupole3f = Multipole3f<5>;
+
 }
 
 #endif //SYMTENSOR_MULTIPOLE_H
diff --git a/include/symtensor/MultipoleBase.h b/include/symtensor/MultipoleBase.h
@@ -153,9 +153,15 @@ namespace symtensor {
         explicit inline constexpr MultipoleBase(const Vector &vector) {
             tensor<1>() = vector;
             [&]<std::size_t... i>(std::index_sequence<i...>) {
-                ((tensor<i + 2>() = TensorType<i + 2>::CartesianProduct(tensor<i + 1>(), vector)), ...);
+                ((tensor<i + 2>() = TensorType<i + 2>::CartesianPower(vector)), ...);
             }(std::make_index_sequence<Order - 1>());
-            // todo: something like this might be better:
+
+            // todo: this seems to under-perform the "naive" approach
+            //            [&]<std::size_t... i>(std::index_sequence<i...>) {
+            //                ((tensor<i + 2>() = TensorType<i + 2>::CartesianProduct(tensor<i + 1>(), vector)), ...);
+            //            }(std::make_index_sequence<Order - 1>());
+
+            // todo: something like this might be more readable:
             //            [&]<std::size_t _0, std::size_t _1, std::size_t... r>(std::index_sequence<_0, _1, r...>) {
             //                ((tensor<r>() = TensorType<r>::CartesianProduct(tensor<r - 1>(), vector)), ...);
             //            }(std::make_index_sequence<Order - 1>());