From 4361121eb8629b42e3cc2b0e702c8e2b9a59711c Mon Sep 17 00:00:00 2001
From: yanmulin <yans1996@outlook.com>
Date: Tue, 18 Jan 2022 12:43:10 -0500
Subject: [PATCH 1/2] finish hw04

---
 CMakeLists.txt     |   5 ++-
 initial_main.cpp   |  88 +++++++++++++++++++++++++++++++++++++
 opt1_math_main.cpp |  88 +++++++++++++++++++++++++++++++++++++
 opt2_soa_main.cpp  | 107 +++++++++++++++++++++++++++++++++++++++++++++
 4 files changed, 287 insertions(+), 1 deletion(-)
 create mode 100644 initial_main.cpp
 create mode 100644 opt1_math_main.cpp
 create mode 100644 opt2_soa_main.cpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 29b152c..0ddd7a9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -6,4 +6,7 @@ if (NOT CMAKE_BUILD_TYPE)
     set(CMAKE_BUILD_TYPE Release)
 endif()
 
-add_executable(main main.cpp)
+add_executable(main opt2_soa_main.cpp)
+
+target_compile_options(main PUBLIC -ffast-math)
+target_compile_options(main PUBLIC -march=native)
\ No newline at end of file
diff --git a/initial_main.cpp b/initial_main.cpp
new file mode 100644
index 0000000..601767f
--- /dev/null
+++ b/initial_main.cpp
@@ -0,0 +1,88 @@
+#include <cstdio>
+#include <cstdlib>
+#include <vector>
+#include <chrono>
+#include <cmath>
+
+float frand() {
+    return (float)rand() / RAND_MAX * 2 - 1;
+}
+
+struct Star {
+    float px, py, pz;
+    float vx, vy, vz;
+    float mass;
+};
+
+std::vector<Star> stars;
+
+void init() {
+    for (int i = 0; i < 48; i++) {
+        stars.push_back({
+            frand(), frand(), frand(),
+            frand(), frand(), frand(),
+            frand() + 1,
+        });
+    }
+}
+
+float G = 0.001;
+float eps = 0.001;
+float dt = 0.01;
+
+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;
+        }
+    }
+    for (auto &star: stars) {
+        star.px += star.vx * dt;
+        star.py += star.vy * dt;
+        star.pz += star.vz * 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;
+            float d2 = dx * dx + dy * dy + dz * dz + eps * eps;
+            energy -= other.mass * star.mass * G / sqrt(d2) / 2;
+        }
+    }
+    return energy;
+}
+
+template <class Func>
+long benchmark(Func const &func) {
+    auto t0 = std::chrono::steady_clock::now();
+    func();
+    auto t1 = std::chrono::steady_clock::now();
+    auto dt = std::chrono::duration_cast<std::chrono::milliseconds>(t1 - t0);
+    return dt.count();
+}
+
+int main() {
+    init();
+    printf("Initial energy: %f\n", calc());
+    auto dt = benchmark([&] {
+        for (int i = 0; i < 100000; i++)
+            step();
+    });
+    printf("Final energy: %f\n", calc());
+    printf("Time elapsed: %ld ms\n", dt);
+    return 0;
+}
diff --git a/opt1_math_main.cpp b/opt1_math_main.cpp
new file mode 100644
index 0000000..95d01d0
--- /dev/null
+++ b/opt1_math_main.cpp
@@ -0,0 +1,88 @@
+#include <cstdio>
+#include <cstdlib>
+#include <vector>
+#include <chrono>
+#include <cmath>
+
+float frand() {
+    return (float)rand() / RAND_MAX * 2 - 1;
+}
+
+struct Star {
+    float px, py, pz;
+    float vx, vy, vz;
+    float mass;
+};
+
+std::vector<Star> stars;
+
+void init() {
+    for (int i = 0; i < 48; i++) {
+        stars.push_back({
+            frand(), frand(), frand(),
+            frand(), frand(), frand(),
+            frand() + 1,
+        });
+    }
+}
+
+float G = 0.001;
+float eps = 0.001;
+float dt = 0.01;
+
+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 *= std::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;
+        }
+    }
+    for (auto &star: stars) {
+        star.px += star.vx * dt;
+        star.py += star.vy * dt;
+        star.pz += star.vz * 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;
+            float d2 = dx * dx + dy * dy + dz * dz + eps * eps;
+            energy -= other.mass * star.mass * G / std::sqrt(d2) / 2;
+        }
+    }
+    return energy;
+}
+
+template <class Func>
+long benchmark(Func const &func) {
+    auto t0 = std::chrono::steady_clock::now();
+    func();
+    auto t1 = std::chrono::steady_clock::now();
+    auto dt = std::chrono::duration_cast<std::chrono::milliseconds>(t1 - t0);
+    return dt.count();
+}
+
+int main() {
+    init();
+    printf("Initial energy: %f\n", calc());
+    auto dt = benchmark([&] {
+        for (int i = 0; i < 100000; i++)
+            step();
+    });
+    printf("Final energy: %f\n", calc());
+    printf("Time elapsed: %ld ms\n", dt);
+    return 0;
+}
diff --git a/opt2_soa_main.cpp b/opt2_soa_main.cpp
new file mode 100644
index 0000000..eafe867
--- /dev/null
+++ b/opt2_soa_main.cpp
@@ -0,0 +1,107 @@
+#include <cstdio>
+#include <cstdlib>
+#include <array>
+#include <chrono>
+#include <cmath>
+
+constexpr size_t g_len = 48;
+
+float G = 0.001;
+float eps = 0.001;
+float dt = 0.01;
+
+struct StarVec {
+    std::array<float, g_len> px, py, pz;
+    std::array<float, g_len> vx, vy, vz;
+    std::array<float, g_len> mass;
+};
+
+StarVec stars;
+
+float frand() {
+    return (float)rand() / RAND_MAX * 2 - 1;
+}
+
+void init() {
+    printf("SOA mode\n");
+    for (int i = 0; i < g_len; 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;
+    }
+}
+
+void step() {
+    float eps2 = eps * eps;
+    for (size_t i=0;i<g_len;i++) {
+        float px = stars.px[i], py = stars.py[i], pz = stars.pz[i];
+        float vx = .0f, vy = .0f, vz = .0f;
+
+        for (size_t j=0;j<g_len;j++) {
+            float dx = stars.px[j] - px;
+            float dy = stars.py[j] - py;
+            float dz = stars.pz[j] - pz;
+
+            float d2 = 1.0 / (dx * dx + dy * dy + dz * dz + eps2);
+            d2 *= std::sqrt(d2);
+            float dt_G_mass = dt * G * stars.mass[j];
+            
+            vx += dx * d2 * dt_G_mass;
+            vy += dy * d2 * dt_G_mass;
+            vz += dz * d2 * dt_G_mass;
+        }
+
+        stars.vx[i] += vx;
+        stars.vy[i] += vy;
+        stars.vz[i] += vz;
+    }
+
+    for (size_t i=0;i<g_len;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 (size_t i=0;i<g_len;i++) {
+        float px = stars.px[i], py = stars.py[i], pz = stars.pz[i];
+        float vx = stars.vx[i], vy = stars.vy[i], vz = stars.vz[i];
+        float v2 = vx * vx + vy * vy + vz * vz;
+        energy += stars.mass[i] * v2 / 2;
+        for (size_t j=0;j<g_len;j++) {
+            float dx = stars.px[j] - px;
+            float dy = stars.py[j] - py;
+            float dz = stars.pz[j] - pz;
+            float d2 = dx * dx + dy * dy + dz * dz + eps * eps;
+            energy -= stars.mass[j] * stars.mass[i] * G / std::sqrt(d2) / 2;
+        }
+    }
+    return energy;
+}
+
+template <class Func>
+long benchmark(Func const &func) {
+    auto t0 = std::chrono::steady_clock::now();
+    func();
+    auto t1 = std::chrono::steady_clock::now();
+    auto dt = std::chrono::duration_cast<std::chrono::milliseconds>(t1 - t0);
+    return dt.count();
+}
+
+int main() {
+    init();
+    printf("Initial energy: %f\n", calc());
+    auto dt = benchmark([&] {
+        for (int i = 0; i < 100000; i++)
+            step();
+    });
+    printf("Final energy: %f\n", calc());
+    printf("Time elapsed: %ld ms\n", dt);
+    return 0;
+}
\ No newline at end of file

From f29a34aae6cccf019efcaffda6be4bf5db6072ea Mon Sep 17 00:00:00 2001
From: yanmulin <yans1996@outlook.com>
Date: Tue, 18 Jan 2022 12:51:31 -0500
Subject: [PATCH 2/2] clean code

---
 opt2_soa_main.cpp | 1 -
 1 file changed, 1 deletion(-)

diff --git a/opt2_soa_main.cpp b/opt2_soa_main.cpp
index eafe867..aa67bdb 100644
--- a/opt2_soa_main.cpp
+++ b/opt2_soa_main.cpp
@@ -23,7 +23,6 @@ float frand() {
 }
 
 void init() {
-    printf("SOA mode\n");
     for (int i = 0; i < g_len; i++) {
         stars.px[i] = frand();
         stars.py[i] = frand();