hw04

CMakeLists.txt

@@ -6,4 +6,8 @@ if (NOT CMAKE_BUILD_TYPE)
     set(CMAKE_BUILD_TYPE Release)
 endif()
 
+find_package(OpenMP REQUIRED)
+
 add_executable(main main.cpp)
+target_compile_options(main PUBLIC -ffast-math -march=native)
+target_link_libraries(main PUBLIC OpenMP::OpenMP_CXX)

answer.md

@@ -0,0 +1,76 @@
+# 优化结果
+
+1. 未优化
+
+```
+Initial energy: -13.414000
+Final energy: -13.356842
+Time elapsed: 1946 ms
+```
+
+2. 将全局命名空间的`sqrt()`替换为`std::sqrt()`
+   - 这样会默认使用`float`的重载版本而不是`double`
+
+```
+Initial energy: -13.414000
+Final energy: -13.356841
+Time elapsed: 1571 ms
+```
+
+3. 使用SOA + SIMD优化
+    - 这样做的话不同对象的相同属性会连续存储，在循环中连续访问的时候便于使用AVX指令集优化。
+
+步骤：
+  - 将Star结构体内的`float`改为`vector<float>`
+  - 将循环中的迭代器改为下标访问
+  - CMakeLists中添加OpenMP包，添加`-ffast-math -march=native`选项启用AVX指令集
+  - 使用`#pragma omp simd`避免pointer aliasing影响优化
+
+```
+Initial energy: -13.414010
+Final energy: -13.356913
+Time elapsed: 234 ms
+```
+
+4. 将循环中的不变量提取出来，这样可以避免重复计算
+
+对于`step()`其中的不变量为`G*dt`和`eps*eps`, 提取到全局变量后加快了速度
+
+```
+Initial energy: -13.414012
+Final energy: -13.356912
+Time elapsed: 217 ms
+```
+5. 进行某些代数化简
+
+```cpp
+float x = stars.mass[j] * Gdt / d2;
+stars.vx[i] += dx * x;
+stars.vy[i] += dy * x;
+stars.vz[i] += dz * x;
+```
+
+```
+Initial energy: -13.414010
+Final energy: -13.356913
+Time elapsed: 201 ms
+```
+
+6. 不使用STL容器，使用C数组
+
+有一定提升，但感觉不是很必要
+
+```cpp
+Initial energy: -13.414012
+Final energy: -13.356912
+Time elapsed: 170 ms
+```
+
+7. 尝试使用循环展开，然而无效果
+
+# 结果&总结
+
+最终结果：`1946ms -> 170ms` （11.4倍）
+
+- 效果最显著优化：开启SIMD指令
+- 其他优化：使用`std::sqrt()`, 进行代数化简、循环无关量提出，不使用STL容器

main.cpp

@@ -8,59 +8,66 @@ float frand() {
     return (float)rand() / RAND_MAX * 2 - 1;
 }
 
-struct Star {
-    float px, py, pz;
-    float vx, vy, vz;
-    float mass;
-};
+constexpr int size = 48;
 
-std::vector<Star> stars;
+struct Stars {
+    float px[size], py[size], pz[size];
+    float vx[size], vy[size], vz[size];
+    float mass[size];
+};
+Stars stars;
 
 void init() {
-    for (int i = 0; i < 48; i++) {
-        stars.push_back({
-            frand(), frand(), frand(),
-            frand(), frand(), frand(),
-            frand() + 1,
-        });
+    for (int i = 0; i < size; i++) {
+        stars.px[i] = frand();
+        stars.py[i] = frand();
+        stars.pz[i] = frand();
+        stars.vx[i] = frand();
+        stars.vy[i] = frand();
+        stars.vz[i] = frand();
+        stars.mass[i] = frand()+1;
     }
 }
 
 float G = 0.001;
 float eps = 0.001;
 float dt = 0.01;
+float Gdt = G*dt;
+float epssqure = eps*eps;
 
 void step() {
-    for (auto &star: stars) {
-        for (auto &other: stars) {
-            float dx = other.px - star.px;
-            float dy = other.py - star.py;
-            float dz = other.pz - star.pz;
-            float d2 = dx * dx + dy * dy + dz * dz + eps * eps;
-            d2 *= sqrt(d2);
-            star.vx += dx * other.mass * G * dt / d2;
-            star.vy += dy * other.mass * G * dt / d2;
-            star.vz += dz * other.mass * G * dt / d2;
+#pragma omp simd
+    for (size_t i = 0; i != size; i++) {
+        for (size_t j = 0; j != size; j++) {
+            float dx = stars.px[j] - stars.px[i];
+            float dy = stars.py[j] - stars.py[i];
+            float dz = stars.pz[j] - stars.pz[i];
+            float d2 = dx * dx + dy * dy + dz * dz + epssqure;
+            d2 *= std::sqrt(d2);
+            float x = stars.mass[j] * Gdt / d2;
+            stars.vx[i] += dx * x;
+            stars.vy[i] += dy * x;
+            stars.vz[i] += dz * x;
         }
     }
-    for (auto &star: stars) {
-        star.px += star.vx * dt;
-        star.py += star.vy * dt;
-        star.pz += star.vz * dt;
+    for (size_t i = 0; i != size; i++) {
+        stars.px[i] += stars.vx[i] * dt;
+        stars.py[i] += stars.vy[i] * dt;
+        stars.pz[i] += stars.vz[i] * dt;
     }
 }
 
 float calc() {
     float energy = 0;
-    for (auto &star: stars) {
-        float v2 = star.vx * star.vx + star.vy * star.vy + star.vz * star.vz;
-        energy += star.mass * v2 / 2;
-        for (auto &other: stars) {
-            float dx = other.px - star.px;
-            float dy = other.py - star.py;
-            float dz = other.pz - star.pz;
+    for (size_t i = 0; i != size; i++) {
+        float v2 = stars.vx[i] * stars.vx[i] + stars.vy[i] * stars.vy[i] + stars.vz[i] * stars.vz[i];
+        energy += stars.mass[i] * v2 / 2;
+        for (size_t j = 0; j != size; j++) {
+            float dx = stars.px[j] - stars.px[i];
+            float dy = stars.py[j] - stars.py[i];
+            float dz = stars.pz[j] - stars.pz[i];
             float d2 = dx * dx + dy * dy + dz * dz + eps * eps;
-            energy -= other.mass * star.mass * G / sqrt(d2) / 2;
+            energy -= stars.mass[j] * stars.mass[i] * G / std::sqrt(d2) / 2;
         }
     }
     return energy;