From 89df99b73085dfdf2306dd356caa25e3614e3fc9 Mon Sep 17 00:00:00 2001 From: AvvAL Date: Sun, 10 Dec 2023 17:45:17 +0300 Subject: [PATCH] done --- src/cl/clion_defines.cl | 1 + src/cl/lbvh.cl | 193 +++++++++++++- src/cl/nbody.cl | 27 +- src/main_lbvh.cpp | 569 +++++++++++++++++++++------------------- 4 files changed, 503 insertions(+), 287 deletions(-) diff --git a/src/cl/clion_defines.cl b/src/cl/clion_defines.cl index c2934ba0..501700fc 100644 --- a/src/cl/clion_defines.cl +++ b/src/cl/clion_defines.cl @@ -21,6 +21,7 @@ gentype clamp (gentype x, float minval, float maxval); gentype degrees (gentype radians); gentype max (gentype x, gentype y); gentype min (gentype x, gentype y); +gentype sqrt (gentype x); gentype mix (gentype x, gentype y, gentype a); gentype radians (gentype degrees); gentype sign (gentype x); diff --git a/src/cl/lbvh.cl b/src/cl/lbvh.cl index e0423d86..4b088573 100644 --- a/src/cl/lbvh.cl +++ b/src/cl/lbvh.cl @@ -130,16 +130,15 @@ morton_t zOrder(float fx, float fy, int i){ if (x < 0 || x >= (1 << NBITS_PER_DIM)) { printf("098245490432590890\n"); -// return 0; + //return 0; } if (y < 0 || y >= (1 << NBITS_PER_DIM)) { printf("432764328764237823\n"); -// return 0; + //return 0; } + morton_t morton_code = (spreadBits(y) << 1) | spreadBits(x); - // TODO - - return 0; + return (morton_code << 32) | i; } __kernel void generateMortonCodes(__global const float *pxs, __global const float *pys, @@ -191,25 +190,132 @@ void __kernel merge(__global const morton_t *as, __global morton_t *as_sorted, u int findSplit(__global const morton_t *codes, int i_begin, int i_end, int bit_index) { - // TODO + if (getBit(codes[i_begin], bit_index) == getBit(codes[i_end - 1], bit_index)) { + return -1; + } + + int start = i_begin; + while (i_begin < i_end) { + int mid = (i_begin + i_end) / 2; + if (getBit(codes[start], bit_index) != getBit(codes[mid], bit_index)) + i_end = mid; + else + i_begin = mid + 1; + } + return i_begin; } void findRegion(int *i_begin, int *i_end, int *bit_index, __global const morton_t *codes, int N, int i_node) { - // TODO + if (i_node < 1 || i_node > N - 2) { + printf("842384298293482\n"); + return; + } + int dir = 0; + int i_bit = NBITS - 1; + for (; i_bit >= 0; --i_bit) { + if (getBit(codes[i_node - 1], i_bit) == 0 && getBit(codes[i_node + 1], i_bit) == 1) { + dir = getBit(codes[i_node], i_bit) == 0 ? -1 : 1; + break; + } + } + + if (dir == 0) { + printf("8923482374983\n"); + return; + } + + int K = NBITS - i_bit; + morton_t pref0 = getBits(codes[i_node], i_bit, K); + + + int start = i_node, end = dir == 1 ? N : -1; + while (start + dir != end) { + int mid = start + (end - start) / 2; + if (getBits(codes[mid], i_bit, K) == pref0) { + start = mid; + } else { + end = mid; + } + } + int i_node_end = start; + *bit_index = i_bit - 1; + + if (dir > 0) { + *i_begin = i_node; + *i_end = i_node_end + 1; + } else { + *i_begin = i_node_end; + *i_end = i_node + 1; + } } void initLBVHNode(__global struct Node *nodes, int i_node, __global const morton_t *codes, int N, __global const float *pxs, __global const float *pys, __global const float *mxs) { - // TODO + __global struct Node* node = &nodes[i_node]; + clear(&node->bbox); + node->mass = 0; + node->cmsx = 0; + node->cmsy = 0; + + if (i_node >= N - 1) { + node->child_left = -1; + node->child_right = -1; + int i_point = i_node - (N - 1); + + int idx = getIndex(codes[i_point]); + float center_mass_x = pxs[idx], center_mass_y = pys[idx]; + float mass = mxs[idx]; + + growPoint(&node->bbox, center_mass_x, center_mass_y); + node->cmsx = center_mass_x; + node->cmsy = center_mass_y; + node->mass = mass; + + return; + } + + int i_begin = 0, i_end = N, bit_index = NBITS - 1; + if (i_node) { + findRegion(&i_begin, &i_end, &bit_index, codes, N, i_node); + } + + bool found = false; + for (int i_bit = bit_index; i_bit >= 0; --i_bit) { + int split = findSplit(codes, i_begin, i_end, i_bit); + if (split < 0) + continue; + + if (split < 1) { + printf("043204230042342\n"); + return; + } + + // проинициализировать nodes[i_node].child_left, nodes[i_node].child_right на основе i_begin, i_end, split + node->child_left = split - i_begin < 2 ? split - 1 + (N - 1) : split - 1; + node->child_right = i_end - split < 2 ? split + (N - 1) : split; + // не забудьте на N-1 сдвинуть индексы, указывающие на листья + + found = true; + break; + } + + if (!found) { + printf("54356549645\n"); + } } __kernel void buidLBVH(__global const float *pxs, __global const float *pys, __global const float *mxs, __global const morton_t *codes, __global struct Node *nodes, int N) { - // TODO + int tree_size = LBVHSize(N); + int i_node = get_global_id(0); + if (i_node >= tree_size) + return; + + initLBVHNode(nodes, i_node, codes, N, pxs, pys, mxs); } void initFlag(__global int *flags, int i_node, __global const struct Node *nodes, int level) @@ -299,7 +405,66 @@ bool barnesHutCondition(float x, float y, __global const struct Node *node) void calculateForce(float x0, float y0, float m0, __global const struct Node *nodes, __global float *force_x, __global float *force_y) { - // TODO + // основная идея ускорения - аггрегировать в узлах дерева веса и центры масс, + // и не спускаться внутрь, если точка запроса не пересекает ноду, а заменить на взаимодействие с ее центром масс + + int stack[2 * NBITS_PER_DIM]; + int stack_size = 0; + stack[stack_size++] = 0; + while (stack_size) { + int i_node = stack[--stack_size]; + __global const struct Node *node = &nodes[i_node]; + + if (isLeaf(node)) { + continue; + } + + // если запрос содержится и а левом и в правом ребенке - то они в одном пикселе + { + __global const struct Node *left = &nodes[node->child_left]; + __global const struct Node *right = &nodes[node->child_right]; + if (contains(&left->bbox, x0, y0) && contains(&right->bbox, x0, y0)) { + if (!equals(&left->bbox, &right->bbox)) { + printf("42357987645432456547\n"); + } + if (!equalsPoint(&left->bbox, x0, y0)) { + printf("5446456456435656\n"); + } + continue; + } + } + int idx_children[] = {node->child_left, node->child_right}; + for (int i = 0; i < 2; ++i) { + int i_child = idx_children[i]; + __global const struct Node *child = &nodes[i_child]; + // С точки зрения ббоксов заходить в ребенка, ббокс которого не пересекаем, не нужно (из-за того, что в листьях у нас точки и они не высовываются за свой регион пространства) + // Но, с точки зрения физики, замена гравитационного влияния всех точек в регионе на взаимодействие с суммарной массой в центре масс - это точное решение только в однородном поле (например, на поверхности земли) + // У нас поле неоднородное, и такая замена - лишь приближение. Чтобы оно было достаточно точным, будем спускаться внутрь ноды, пока она не станет похожа на точечное тело (маленький размер ее ббокса относительно нашего расстояния до центра масс ноды) + if (!contains(&child->bbox, x0, y0) && barnesHutCondition(x0, y0, child)) { + float dx = child->cmsx - x0; + float dy = child->cmsy - y0; + float dr2 = max(100.f, dx * dx + dy * dy); + + float dr2_inv = 1.f / dr2; + float dr_inv = sqrt(dr2_inv); + + float ex = dx * dr_inv; + float ey = dy * dr_inv; + + float fx = ex * dr2_inv * GRAVITATIONAL_FORCE; + float fy = ey * dr2_inv * GRAVITATIONAL_FORCE; + + *force_x += child->mass * fx; + *force_y += child->mass * fy; + } else { + stack[stack_size++] = i_child; + if (stack_size >= 2 * NBITS_PER_DIM) { + printf("0420392384283\n"); + return; + } + } + } + } } __kernel void calculateForces( @@ -311,7 +476,13 @@ __kernel void calculateForces( int N, int t) { - // TODO + int i = get_global_id(0); + float x0 = pxs[i]; + float y0 = pys[i]; + float m0 = mxs[i]; + __global float *force_x = dvx2d + t * N + i; + __global float *force_y = dvy2d + t * N + i; + calculateForce(x0, y0, m0, nodes, force_x, force_y); } __kernel void integrate( diff --git a/src/cl/nbody.cl b/src/cl/nbody.cl index c70482cc..a90c034a 100644 --- a/src/cl/nbody.cl +++ b/src/cl/nbody.cl @@ -26,7 +26,32 @@ __kernel void nbody_calculate_force_global( float y0 = pys[i]; float m0 = mxs[i]; - // TODO + for (int j = 0; j < N; ++j) { + + if (j == i) { + continue; + } + + float x1 = pxs[j]; + float y1 = pys[j]; + float m1 = mxs[j]; + + float dx = x1 - x0; + float dy = y1 - y0; + float dr2 = max(100.f, dx * dx + dy * dy); + + float dr2_inv = 1.f / dr2; + float dr_inv = sqrt(dr2_inv); + + float ex = dx * dr_inv; + float ey = dy * dr_inv; + + float fx = ex * dr2_inv * GRAVITATIONAL_FORCE; + float fy = ey * dr2_inv * GRAVITATIONAL_FORCE; + + dvx[i] += m1 * fx; + dvy[i] += m1 * fy; + } } __kernel void nbody_integrate( diff --git a/src/main_lbvh.cpp b/src/main_lbvh.cpp index fc1d7dc1..400a63e6 100644 --- a/src/main_lbvh.cpp +++ b/src/main_lbvh.cpp @@ -17,28 +17,28 @@ // может понадобиться поменять индекс локально чтобы выбрать GPU если у вас более одного девайса -#define OPENCL_DEVICE_INDEX 0 +#define OPENCL_DEVICE_INDEX 1 // TODO включить чтобы начали запускаться тесты -#define ENABLE_TESTING 0 +#define ENABLE_TESTING 1 // имеет смысл отключать при оффлайн симуляции больших N, но в итоговом решении стоит оставить #define EVALUATE_PRECISION 1 // удобно включить при локальном тестировании -#define ENABLE_GUI 0 +#define ENABLE_GUI 1 // сброс картинок симуляции на диск #define SAVE_IMAGES 0 // TODO на сервер лучше коммитить самую простую конфигурацию. Замеры по времени получатся нерелевантные, но зато быстрее отработает CI // TODO локально интересны замеры на самой сложной версии, которую получится дождаться -#define NBODY_INITIAL_STATE_COMPLEXITY 0 -//#define NBODY_INITIAL_STATE_COMPLEXITY 1 -//#define NBODY_INITIAL_STATE_COMPLEXITY 2 +// #define NBODY_INITIAL_STATE_COMPLEXITY 0 +// #define NBODY_INITIAL_STATE_COMPLEXITY 1 +#define NBODY_INITIAL_STATE_COMPLEXITY 2 // использовать lbvh для построения начального состояния. Нужно на очень больших N (>1000000) -#define ENABLE_LBVH_STATE_INITIALIZATION 0 +#define ENABLE_LBVH_STATE_INITIALIZATION 1 //////////////////////////////////////////////////////////// @@ -47,7 +47,8 @@ //////////////////////////////////////////////////////////// -struct Color { +struct Color +{ unsigned char r, g, b; }; @@ -62,22 +63,26 @@ const Color WHITE{255, 255, 255}; const double GRAVITATIONAL_FORCE = 0.0001; -struct DeltaState { +struct DeltaState +{ std::vector dvx2d; std::vector dvy2d; }; -struct State { +struct State +{ + + State() { + } - State() {} State(int N) : pxs(N) , pys(N) , vxs(N) , vys(N) , mxs(N) - , coord_shift(0) - {} + , coord_shift(0) { + } std::vector pxs; std::vector pys; @@ -90,7 +95,8 @@ struct State { int coord_shift; }; -struct Point { +struct Point +{ int x, y; bool operator==(const Point &rhs) const { @@ -101,30 +107,26 @@ struct Point { return !(rhs == *this); } - Point & operator+=(const Point &other) - { + Point &operator+=(const Point &other) { x += other.x; y += other.y; return *this; } - Point & operator-=(const Point &other) - { + Point &operator-=(const Point &other) { x -= other.x; y -= other.y; return *this; } - Point operator-(const Point &other) const - { + Point operator-(const Point &other) const { Point result = *this; result -= other; return result; } }; -void bresenham(std::vector &line_points, const Point &from, const Point &to) -{ +void bresenham(std::vector &line_points, const Point &from, const Point &to) { line_points.clear(); double x1 = from.x; @@ -177,11 +179,11 @@ void bresenham(std::vector &line_points, const Point &from, const Point & } //Преобразует 0xbbbb в 0x0b0b0b0b -int spreadBits(int word){ - word = (word ^ (word << 8 )) & 0x00ff00ff; - word = (word ^ (word << 4 )) & 0x0f0f0f0f; - word = (word ^ (word << 2 )) & 0x33333333; - word = (word ^ (word << 1 )) & 0x55555555; +int spreadBits(int word) { + word = (word ^ (word << 8)) & 0x00ff00ff; + word = (word ^ (word << 4)) & 0x0f0f0f0f; + word = (word ^ (word << 2)) & 0x33333333; + word = (word ^ (word << 1)) & 0x55555555; return word; } @@ -189,30 +191,30 @@ using morton_t = uint64_t; const int NBITS_PER_DIM = 16; const int NBITS = NBITS_PER_DIM /*x dimension*/ + NBITS_PER_DIM /*y dimension*/ + 32 /*index augmentation*/; //Convert xy coordinate to a 32 bit morton/z order code + 32 bit index augmentation for distinguishing between duplicates -morton_t zOrder(const Point &coord, int i){ - if (coord.x < 0 || coord.x >= (1 << NBITS_PER_DIM)) throw std::runtime_error("098245490432590890"); - if (coord.y < 0 || coord.y >= (1 << NBITS_PER_DIM)) throw std::runtime_error("432764328764237823"); - int x = coord.x; - int y = coord.y; - - throw std::runtime_error("not implemented"); -// morton_t morton_code = TODO -// -// // augmentation -// return (morton_code << 32) | i; +morton_t zOrder(const Point &coord, int i) { + if (coord.x < 0 || coord.x >= (1 << NBITS_PER_DIM)) + throw std::runtime_error("098245490432590890"); + if (coord.y < 0 || coord.y >= (1 << NBITS_PER_DIM)) + throw std::runtime_error("432764328764237823"); + + int x = spreadBits(coord.x); + int y = spreadBits(coord.y); + + morton_t morton_code = (y << 1) | x; + return (morton_code << 32) | i; } #pragma pack (push, 1) -struct BBox { +struct BBox +{ - BBox() - { + BBox() { clear(); } - explicit BBox(const Point &point) : BBox() - { + explicit BBox(const Point &point) + : BBox() { grow(point); } @@ -221,59 +223,50 @@ struct BBox { grow(fx, fy); } - void clear() - { + void clear() { minx = std::numeric_limits::max(); maxx = std::numeric_limits::lowest(); miny = minx; maxy = maxx; } - bool contains(const Point &point) const - { + bool contains(const Point &point) const { return point.x >= minx && point.x <= maxx && point.y >= miny && point.y <= maxy; } - bool contains(float fx, float fy) const - { + bool contains(float fx, float fy) const { int x = fx + 0.5; int y = fy + 0.5; return x >= minx && x <= maxx && y >= miny && y <= maxy; } - bool empty() const - { + bool empty() const { return minx > maxx; } - bool operator==(const BBox &other) const - { + bool operator==(const BBox &other) const { return minx == other.minx && maxx == other.maxx && miny == other.miny && maxy == other.maxy; } - bool operator!=(const BBox &other) const - { + bool operator!=(const BBox &other) const { return !(*this == other); } - void grow(const Point &point) - { + void grow(const Point &point) { minx = std::min(minx, point.x); maxx = std::max(maxx, point.x); miny = std::min(miny, point.y); maxy = std::max(maxy, point.y); } - void grow(float fx, float fy) - { + void grow(float fx, float fy) { minx = std::min(minx, int(fx + 0.5)); maxx = std::max(maxx, int(fx + 0.5)); miny = std::min(miny, int(fy + 0.5)); maxy = std::max(maxy, int(fy + 0.5)); } - void grow(const BBox &other) - { + void grow(const BBox &other) { grow(other.min()); grow(other.max()); } @@ -287,21 +280,25 @@ struct BBox { Point max() const { return Point{maxx, maxy}; } private: - int minx, maxx; int miny, maxy; }; -struct Node { +bool feq(float x, float y) { + float maxXY = std::max(1e-8f, std::max(std::fabs(x) , std::fabs(y))); + return std::fabs(x - y) < 1e-3f * maxXY; +} +struct Node +{ bool hasLeftChild() const { return child_left >= 0; } bool hasRightChild() const { return child_right >= 0; } bool isLeaf() const { return !hasLeftChild() && !hasRightChild(); } bool operator==(const Node &other) const { - return std::tie(child_left, child_right, bbox, mass, cmsx, cmsy) - == std::tie(other.child_left, other.child_right, other.bbox, other.mass, other.cmsx, other.cmsy); + return child_left == other.child_left && child_right == other.child_right && bbox == other.bbox && + feq(mass, other.mass) && feq(cmsx, other.cmsx) && feq(cmsy,other.cmsy); } bool operator!=(const Node &other) const { @@ -318,19 +315,16 @@ struct Node { }; #pragma pack (pop) -morton_t getBits(morton_t morton_code, int bit_index, int prefix_size) -{ +morton_t getBits(morton_t morton_code, int bit_index, int prefix_size) { return (morton_code >> bit_index) & ((1ull << prefix_size) - 1ull); } -int getBit(morton_t morton_code, int bit_index) -{ +int getBit(morton_t morton_code, int bit_index) { return (morton_code >> bit_index) & 1; } // из аугментированного мортоновского кода можно извлечь индекс в изначальном массиве -int getIndex(morton_t morton_code) -{ +int getIndex(morton_t morton_code) { morton_t mask = 1; mask = (mask << 32) - 1; return morton_code & mask; @@ -338,32 +332,33 @@ int getIndex(morton_t morton_code) // N листьев, N-1 внутренних нод int LBVHSize(int N) { - return N + N-1; + return N + N - 1; } using points_mass_functor = std::function(int)>; -Point makePoint(float x, float y) -{ + +Point makePoint(float x, float y) { if (x < 0.f || y < 0.f) { throw std::runtime_error("0432959435934534"); } return Point{int(x + 0.5f), int(y + 0.5f)}; } -void initLBVHNode(std::vector &nodes, int i_node, const std::vector &codes, const points_mass_functor &points_mass_array); +void initLBVHNode(std::vector &nodes, int i_node, const std::vector &codes, + const points_mass_functor &points_mass_array); void buildBBoxes(std::vector &nodes, std::vector &flags, int N, bool use_omp = false); void buildBBoxesRecursive(std::vector &nodes, Node &root); void drawLBVH(images::Image &canvas, const std::vector &nodes, int coord_shift = 0); -using interactive_callback_t = std::function &, const std::vector &, const std::vector&)>; +using interactive_callback_t = std::function &, const std::vector &, + const std::vector &)>; // https://en.wikipedia.org/wiki/Barnes%E2%80%93Hut_simulation -bool barnesHutCondition(float x, float y, const Node &node) -{ +bool barnesHutCondition(float x, float y, const Node &node) { float dx = x - node.cmsx; float dy = y - node.cmsy; float s = std::max(node.bbox.maxX() - node.bbox.minX(), node.bbox.maxY() - node.bbox.minY()); - float d2 = dx*dx + dy*dy; + float d2 = dx * dx + dy * dy; float thresh = 0.5; // возвращаем true, если находимся от ноды достаточно далеко, чтобы можно было ее считать примерно точечной @@ -371,18 +366,16 @@ bool barnesHutCondition(float x, float y, const Node &node) return s * s < d2 * thresh * thresh; } -void calculateForce(float x0, float y0, float m0, const std::vector &nodes, float *force_x, float *force_y) -{ +void calculateForce(float x0, float y0, float m0, const std::vector &nodes, float *force_x, float *force_y) { // основная идея ускорения - аггрегировать в узлах дерева веса и центры масс, // и не спускаться внутрь, если точка запроса не пересекает ноду, а заменить на взаимодействие с ее центром масс int stack[2 * NBITS_PER_DIM]; int stack_size = 0; - // TODO кладем корень на стек - throw std::runtime_error("not implemented"); - /* while (stack_size) { - // TODO берем ноду со стека - throw std::runtime_error("not implemented"); + stack[stack_size++] = 0; + while (stack_size) { + int i_node = stack[--stack_size]; + auto node = nodes[i_node]; if (node.isLeaf()) { continue; @@ -402,28 +395,39 @@ void calculateForce(float x0, float y0, float m0, const std::vector &nodes continue; } } - for (int i_child : {node.child_left, node.child_right}) { const Node &child = nodes[i_child]; // С точки зрения ббоксов заходить в ребенка, ббокс которого не пересекаем, не нужно (из-за того, что в листьях у нас точки и они не высовываются за свой регион пространства) // Но, с точки зрения физики, замена гравитационного влияния всех точек в регионе на взаимодействие с суммарной массой в центре масс - это точное решение только в однородном поле (например, на поверхности земли) // У нас поле неоднородное, и такая замена - лишь приближение. Чтобы оно было достаточно точным, будем спускаться внутрь ноды, пока она не станет похожа на точечное тело (маленький размер ее ббокса относительно нашего расстояния до центра масс ноды) if (!child.bbox.contains(x0, y0) && barnesHutCondition(x0, y0, child)) { - // TODO посчитать взаимодействие точки с центром масс ноды - throw std::runtime_error("not implemented"); + float dx = child.cmsx - x0; + float dy = child.cmsy - y0; + float dr2 = std::max(100.f, dx * dx + dy * dy); + + float dr2_inv = 1.f / dr2; + float dr_inv = sqrt(dr2_inv); + + float ex = dx * dr_inv; + float ey = dy * dr_inv; + + float fx = ex * dr2_inv * GRAVITATIONAL_FORCE; + float fy = ey * dr2_inv * GRAVITATIONAL_FORCE; + + *force_x += child.mass * fx; + *force_y += child.mass * fy; } else { - // TODO кладем ребенка на стек - throw std::runtime_error("not implemented"); + stack[stack_size++] = i_child; if (stack_size >= 2 * NBITS_PER_DIM) { throw std::runtime_error("0420392384283"); } } } - }*/ + } } -void integrate(int i, std::vector &pxs, std::vector &pys, std::vector &vxs, std::vector &vys, float *dvx, float *dvy, int coord_shift) -{ +void integrate(int i, std::vector &pxs, std::vector &pys, std::vector &vxs, + std::vector &vys, float *dvx, float *dvy, int coord_shift) { vxs[i] += dvx[i]; vys[i] += dvy[i]; pxs[i] += vxs[i]; @@ -449,8 +453,8 @@ void integrate(int i, std::vector &pxs, std::vector &pys, std::vec } // in: initial conditions, out: 2D array for integration -void nbody_cpu_lbvh(DeltaState &delta_state, State &initial_state, int N, int NT, const interactive_callback_t *interactive_callback = nullptr) -{ +void nbody_cpu_lbvh(DeltaState &delta_state, State &initial_state, int N, int NT, + const interactive_callback_t *interactive_callback = nullptr) { int NT_interactive = interactive_callback ? 1 : NT; delta_state.dvx2d.assign(N * NT_interactive, 0.0); @@ -474,11 +478,11 @@ void nbody_cpu_lbvh(DeltaState &delta_state, State &initial_state, int N, int NT for (int t = 0; t < NT; ++t) { int t_interactive = interactive_callback ? 0 : t; - float * dvx = &delta_state.dvx2d[t_interactive * N]; - float * dvy = &delta_state.dvy2d[t_interactive * N]; + float *dvx = &delta_state.dvx2d[t_interactive * N]; + float *dvy = &delta_state.dvy2d[t_interactive * N]; -// инициализируем мортоновские коды -#pragma omp parallel for + // инициализируем мортоновские коды + #pragma omp parallel for for (int i = 0; i < N; ++i) { codes[i] = zOrder(makePoint(pxs[i], pys[i]), i); } @@ -486,8 +490,8 @@ void nbody_cpu_lbvh(DeltaState &delta_state, State &initial_state, int N, int NT // упорядочиваем тела по z-curve std::sort(codes.begin(), codes.end()); -// строим LBVH -#pragma omp parallel for + // строим LBVH + #pragma omp parallel for for (int i_node = 0; i_node < tree_size; ++i_node) { initLBVHNode(nodes, i_node, codes, points_mass_array); } @@ -495,7 +499,7 @@ void nbody_cpu_lbvh(DeltaState &delta_state, State &initial_state, int N, int NT // инициализируем ббоксы и массы buildBBoxes(nodes, buffer, N, false/*omp here can cause stuttering on start of simulation*/); -#pragma omp parallel for + #pragma omp parallel for for (int i = 0; i < N; ++i) { float x0 = pxs[i]; float y0 = pys[i]; @@ -508,7 +512,7 @@ void nbody_cpu_lbvh(DeltaState &delta_state, State &initial_state, int N, int NT (*interactive_callback)(pxs, pys, nodes); } -#pragma omp parallel for + #pragma omp parallel for for (int i = 0; i < N; ++i) { integrate(i, pxs, pys, vxs, vys, dvx, dvy, initial_state.coord_shift); } @@ -521,8 +525,8 @@ void nbody_cpu_lbvh(DeltaState &delta_state, State &initial_state, int N, int NT } // in: initial conditions, out: 2D array for integration -void nbody_cpu(DeltaState &delta_state, State &initial_state, int N, int NT, const interactive_callback_t *interactive_callback = nullptr) -{ +void nbody_cpu(DeltaState &delta_state, State &initial_state, int N, int NT, + const interactive_callback_t *interactive_callback = nullptr) { int NT_interactive = interactive_callback ? 1 : NT; delta_state.dvx2d.assign(N * NT_interactive, 0.0); @@ -539,11 +543,11 @@ void nbody_cpu(DeltaState &delta_state, State &initial_state, int N, int NT, con for (int t = 0; t < NT; ++t) { int t_interactive = interactive_callback ? 0 : t; - float * dvx = &delta_state.dvx2d[t_interactive * N]; - float * dvy = &delta_state.dvy2d[t_interactive * N]; + float *dvx = &delta_state.dvx2d[t_interactive * N]; + float *dvy = &delta_state.dvy2d[t_interactive * N]; -// to be kernel 1 -#pragma omp parallel for + // to be kernel 1 + #pragma omp parallel for for (int i = 0; i < N; ++i) { float x0 = pxs[i]; float y0 = pys[i]; @@ -580,8 +584,8 @@ void nbody_cpu(DeltaState &delta_state, State &initial_state, int N, int NT, con (*interactive_callback)(pxs, pys, {}); } -// to be kernel 2 -#pragma omp parallel for + // to be kernel 2 + #pragma omp parallel for for (int i = 0; i < N; ++i) { vxs[i] += dvx[i]; vys[i] += dvy[i]; @@ -596,8 +600,8 @@ void nbody_cpu(DeltaState &delta_state, State &initial_state, int N, int NT, con } } -void nbody_gpu_lbvh(DeltaState &delta_state, State &initial_state, int N, int NT, const interactive_callback_t *interactive_callback = nullptr) -{ +void nbody_gpu_lbvh(DeltaState &delta_state, State &initial_state, int N, int NT, + const interactive_callback_t *interactive_callback = nullptr) { int tree_size = LBVHSize(N); std::vector nodes(tree_size); @@ -701,7 +705,7 @@ void nbody_gpu_lbvh(DeltaState &delta_state, State &initial_state, int N, int NT N, level); int n_updated; - flags_gpu.readN(&n_updated, 1, N-1); + flags_gpu.readN(&n_updated, 1, N - 1); if (!n_updated) break; @@ -736,8 +740,8 @@ void nbody_gpu_lbvh(DeltaState &delta_state, State &initial_state, int N, int NT dvy2d_gpu.readN(dvy2d.data(), N * NT_interactive); } -void nbody_gpu(DeltaState &delta_state, State &initial_state, int N, int NT, const interactive_callback_t *interactive_callback = nullptr) -{ +void nbody_gpu(DeltaState &delta_state, State &initial_state, int N, int NT, + const interactive_callback_t *interactive_callback = nullptr) { int NT_interactive = interactive_callback ? 1 : NT; delta_state.dvx2d.assign(N * NT_interactive, 0.0); @@ -813,8 +817,7 @@ void nbody_gpu(DeltaState &delta_state, State &initial_state, int N, int NT, con dvy2d_gpu.readN(dvy2d.data(), N * NT_interactive); } -State makeSimpleState() -{ +State makeSimpleState() { State state; state.pxs = {-25.f, 25.f, 100.f}; state.pys = {-0.f, 0.f, 0.f}; @@ -824,8 +827,7 @@ State makeSimpleState() return state; } -State makeRandomState(int N, int minx, int maxx, int miny, int maxy) -{ +State makeRandomState(int N, int minx, int maxx, int miny, int maxy) { State result(N); int w = maxx - minx; @@ -835,8 +837,10 @@ State makeRandomState(int N, int minx, int maxx, int miny, int maxy) result.pxs[i] = std::rand() % w + minx; result.pys[i] = std::rand() % h + miny; - result.vxs[i] = 0; std::rand() % 3 - 1; - result.vys[i] = 0; std::rand() % 3 - 1; + result.vxs[i] = 0; + std::rand() % 3 - 1; + result.vys[i] = 0; + std::rand() % 3 - 1; result.mxs[i] = std::rand() % 20 + 1; } @@ -844,16 +848,15 @@ State makeRandomState(int N, int minx, int maxx, int miny, int maxy) return result; } -State makeCircularState() -{ +State makeCircularState() { const int N = 1000; State result(N); - for (int i = 0; i < N-1; ++i) { + for (int i = 0; i < N - 1; ++i) { int r = 20; // circle sampling - float angle = 3.14159 * 2 / (N-1) * i; + float angle = 3.14159 * 2 / (N - 1) * i; int x = r * std::sin(angle); int y = r * std::cos(angle); result.pxs[i] = x; @@ -876,7 +879,7 @@ State makeCircularState() for (int i = 0; i < N - 1; ++i) { float rx = result.pxs[i]; float ry = result.pys[i]; - float r = std::sqrt(rx*rx + ry*ry); + float r = std::sqrt(rx * rx + ry * ry); { rx /= r; ry /= r; @@ -895,8 +898,7 @@ State makeCircularState() return result; } -State makeGalacticState(int N, int s0, int s1) -{ +State makeGalacticState(int N, int s0, int s1) { State result(N); int minx = -s1; @@ -932,8 +934,8 @@ State makeGalacticState(int N, int s0, int s1) DeltaState delta_state; auto tmp = result; -// для больших N (миллион) дождаться даже одного брутфорс шага проблематично -#if ENABLE_LBVH_STATE_INITIALIZATION + // для больших N (миллион) дождаться даже одного брутфорс шага проблематично + #if ENABLE_LBVH_STATE_INITIALIZATION { tmp.coord_shift = 1 << (NBITS_PER_DIM - 1); for (int i = 0; i < N; ++i) { @@ -942,14 +944,14 @@ State makeGalacticState(int N, int s0, int s1) } nbody_cpu_lbvh(delta_state, tmp, N, 1); } -#else + #else nbody_cpu(delta_state, tmp, N, 1); -#endif + #endif for (int i = 0; i < N - 1; ++i) { float rx = result.pxs[i]; float ry = result.pys[i]; - float r = std::sqrt(rx*rx + ry*ry); + float r = std::sqrt(rx * rx + ry * ry); { rx /= r; ry /= r; @@ -968,31 +970,35 @@ State makeGalacticState(int N, int s0, int s1) return result; } -int findSplit(const std::vector &codes, int i_begin, int i_end, int bit_index) -{ +int findSplit(const std::vector &codes, int i_begin, int i_end, int bit_index) { // Если биты в начале и в конце совпадают, то этот бит незначащий - if (getBit(codes[i_begin], bit_index) == getBit(codes[i_end-1], bit_index)) { + if (getBit(codes[i_begin], bit_index) == getBit(codes[i_end - 1], bit_index)) { return -1; } // наивная версия, линейный поиск, можно использовать для отладки бинпоиска + // int linear_i = 0; // for (int i = i_begin + 1; i < i_end; ++i) { - // int a = getBit(codes[i-1].first, bit_index); - // int b = getBit(codes[i].first, bit_index); + // int a = getBit(codes[i-1], bit_index); + // int b = getBit(codes[i], bit_index); // if (a < b) { - // return i; + // linear_i = i; // } // } - // TODO бинпоиск для нахождения разбиения области ответственности ноды - throw std::runtime_error("not implemented"); - - // избыточно, так как на входе в функцию проверили, что ответ существует, но приятно иметь sanity-check на случай если набагали - throw std::runtime_error("4932492039458209485"); + int start = i_begin; + while (i_begin < i_end) { + int mid = (i_begin + i_end) / 2; + if (getBit(codes[start], bit_index) != getBit(codes[mid], bit_index)) + i_end = mid; + else + i_begin = mid + 1; + } + return i_begin; } -void buildLBVHRecursive(std::vector &nodes, const std::vector &codes, const std::vector &points, int i_begin, int i_end, int bit_index) -{ +void buildLBVHRecursive(std::vector &nodes, const std::vector &codes, const std::vector &points, + int i_begin, int i_end, int bit_index) { int i_node = nodes.size(); nodes.emplace_back(); @@ -1010,7 +1016,8 @@ void buildLBVHRecursive(std::vector &nodes, const std::vector &c for (int i_bit = bit_index; i_bit >= 0; --i_bit) { int split = findSplit(codes, i_begin, i_end, i_bit); - if (split < 0) continue; + if (split < 0) + continue; nodes[i_node].child_left = nodes.size(); buildLBVHRecursive(nodes, codes, points, i_begin, split, i_bit - 1); @@ -1025,8 +1032,7 @@ void buildLBVHRecursive(std::vector &nodes, const std::vector &c throw std::runtime_error("043242304023: potentially found duplicate morton code"); } -void findRegion(int *i_begin, int *i_end, int *bit_index, const std::vector &codes, int i_node) -{ +void findRegion(int *i_begin, int *i_end, int *bit_index, const std::vector &codes, int i_node) { int N = codes.size(); if (i_node < 1 || i_node > N - 2) { throw std::runtime_error("842384298293482"); @@ -1036,10 +1042,12 @@ void findRegion(int *i_begin, int *i_end, int *bit_index, const std::vector= 0; --i_bit) { - // TODO найти dir и значащий бит - throw std::runtime_error("not implemented"); + if (getBit(codes[i_node - 1], i_bit) == 0 && getBit(codes[i_node + 1], i_bit) == 1) { + dir = getBit(codes[i_node], i_bit) == 0 ? -1 : 1; + break; + } } if (dir == 0) { @@ -1055,20 +1063,26 @@ void findRegion(int *i_begin, int *i_end, int *bit_index, const std::vector= 0 && i < int(codes.size()); i += dir) { - // if (getBits(codes[i], i_bit, K) == pref0) { - // i_node_end = i; - // } else { - // break; - // } - // } - // if (i_node_end == -1) { - // throw std::runtime_error("47248457284332098"); - // } - - // TODO бинпоиск зоны ответственности - throw std::runtime_error("not implemented"); - + // for (int i = i_node; i >= 0 && i < int(codes.size()); i += dir) { + // if (getBits(codes[i], i_bit, K) == pref0) { + // i_node_end = i; + // } else { + // break; + // } + // } + // if (i_node_end == -1) { + // throw std::runtime_error("47248457284332098"); + // } + int start = i_node, end = dir == 1 ? N : -1; + while (start + dir != end) { + int mid = start + (end - start) / 2; + if (getBits(codes[mid], i_bit, K) == pref0) { + start = mid; + } else { + end = mid; + } + } + i_node_end = start; *bit_index = i_bit - 1; if (dir > 0) { @@ -1080,8 +1094,8 @@ void findRegion(int *i_begin, int *i_end, int *bit_index, const std::vector &nodes, int i_node, const std::vector &codes, const points_mass_functor &points_mass_array) -{ +void initLBVHNode(std::vector &nodes, int i_node, const std::vector &codes, + const points_mass_functor &points_mass_array) { // инициализация ссылок на соседей для нод lbvh // если мы лист, то просто инициализируем минус единицами (нет детей), иначе ищем своб зону ответственности и запускаем на ней findSplit // можно заполнить пропуски в виде тудушек, можно реализовать с чистого листа самостоятельно, если так проще @@ -1096,10 +1110,10 @@ void initLBVHNode(std::vector &nodes, int i_node, const std::vector= N-1) { + if (i_node >= N - 1) { nodes[i_node].child_left = -1; nodes[i_node].child_right = -1; - int i_point = i_node - (N-1); + int i_point = i_node - (N - 1); float center_mass_x, center_mass_y; float mass; @@ -1115,32 +1129,28 @@ void initLBVHNode(std::vector &nodes, int i_node, const std::vector= 0; --i_bit) { - /* - int split = TODO - if (split < 0) continue; + int split = findSplit(codes, i_begin, i_end, i_bit); + if (split < 0) + continue; if (split < 1) { throw std::runtime_error("043204230042342"); } - */ - throw std::runtime_error("not implemented"); - // TODO проинициализировать nodes[i_node].child_left, nodes[i_node].child_right на основе i_begin, i_end, split + // проинициализировать nodes[i_node].child_left, nodes[i_node].child_right на основе i_begin, i_end, split + nodes[i_node].child_left = split - i_begin < 2 ? split - 1 + (N - 1) : split - 1; + nodes[i_node].child_right = i_end - split < 2 ? split + (N - 1) : split; // не забудьте на N-1 сдвинуть индексы, указывающие на листья - throw std::runtime_error("not implemented"); - - found = true; break; } @@ -1150,25 +1160,25 @@ void initLBVHNode(std::vector &nodes, int i_node, const std::vector &nodes, const std::vector &codes, const std::vector &points) -{ +void buildLBVH(std::vector &nodes, const std::vector &codes, const std::vector &points) { const int N = codes.size(); int tree_size = LBVHSize(N); nodes.resize(tree_size); - const points_mass_functor points_mass_array = [&](int i) { return std::make_tuple((float) points[i].x, (float) points[i].y, 1.f); }; + const points_mass_functor points_mass_array = [&](int i) { + return std::make_tuple((float) points[i].x, (float) points[i].y, 1.f); + }; // можно раскомментировать и будет работать, но для дебага удобнее оставить однопоточную версию - // #pragma omp parallel for + #pragma omp parallel for for (int i_node = 0; i_node < tree_size; ++i_node) { initLBVHNode(nodes, i_node, codes, points_mass_array); } } -void printMortonCodes(const std::vector &codes) -{ +void printMortonCodes(const std::vector &codes) { std::cout << "morton codes: \n"; - for (int bit_index = NBITS-1; bit_index >= 0; --bit_index) { + for (int bit_index = NBITS - 1; bit_index >= 0; --bit_index) { for (int i = 0; i < (int) codes.size(); ++i) { int bit = getBit(codes[i], bit_index); std::cout << bit << " "; @@ -1178,8 +1188,7 @@ void printMortonCodes(const std::vector &codes) std::cout << std::flush; } -void growNode(Node &root, std::vector &nodes) -{ +void growNode(Node &root, std::vector &nodes) { const Node &left = nodes[root.child_left]; const Node &right = nodes[root.child_right]; @@ -1199,9 +1208,9 @@ void growNode(Node &root, std::vector &nodes) root.cmsy = (left.cmsy * m0 + right.cmsy * m1) / root.mass; } -void buildBBoxesRecursive(std::vector &nodes, Node &root) -{ - if (root.isLeaf()) return; +void buildBBoxesRecursive(std::vector &nodes, Node &root) { + if (root.isLeaf()) + return; buildBBoxesRecursive(nodes, nodes[root.child_left]); buildBBoxesRecursive(nodes, nodes[root.child_right]); @@ -1209,8 +1218,7 @@ void buildBBoxesRecursive(std::vector &nodes, Node &root) growNode(root, nodes); } -void initFlag(std::vector &flags, int i_node, std::vector &nodes, int level) -{ +void initFlag(std::vector &flags, int i_node, std::vector &nodes, int level) { flags[i_node] = -1; Node &node = nodes[i_node]; @@ -1230,34 +1238,39 @@ void initFlag(std::vector &flags, int i_node, std::vector &nodes, int } } -void buildBBoxes(std::vector &nodes, std::vector &flags, int N, bool use_omp) -{ - flags.resize(N-1); +void buildBBoxes(std::vector &nodes, std::vector &flags, int N, bool use_omp) { + flags.resize(N - 1); // NBITS раз проходимся по всему дереву и инициализируем только те ноды, у которых проинициализированы ббоксы обоих детей // не самый оптимальный вариант (O(NlogN) вместо O(N)), зато легко переложить на GPU for (int level = 0; level < NBITS; ++level) { -// знаем, что листья располагаются после N-1 внутренних нод дерева, а для них уже ббоксы посчитаны -> можно пропустить -// не сработает для рекурсивно построенного дерева, там такого порядка не вводили + // знаем, что листья располагаются после N-1 внутренних нод дерева, а для них уже ббоксы посчитаны -> можно пропустить + // не сработает для рекурсивно построенного дерева, там такого порядка не вводили -// чтобы не было гонки в многопоточном режиме (и, по аналогии, потом на видеокарте), в первом проходе отметим ноды, которые нужно обновить, и только на втором проходе обновим -#pragma omp parallel for if(use_omp) - for (int i_node = 0; i_node < N-1; ++i_node) { + // чтобы не было гонки в многопоточном режиме (и, по аналогии, потом на видеокарте), в первом проходе отметим ноды, которые нужно обновить, и только на втором проходе обновим + #pragma omp parallel for if(use_omp) + for (int i_node = 0; i_node < N - 1; ++i_node) { initFlag(flags, i_node, nodes, level); } int n_updated = 0; -#pragma omp parallel for if(use_omp) reduction(+:n_updated) - for (int i_node = 0; i_node < N-1; ++i_node) { - // TODO если находимся на нужном уровне (нужный flag), проинициализируем ббокс и центр масс ноды -// if (TODO) { -// TODO -// ++n_updated; -// } - + #pragma omp parallel for if(use_omp) reduction(+:n_updated) + for (int i_node = 0; i_node < N - 1; ++i_node) { + // если находимся на нужном уровне (нужный flag), проинициализируем ббокс и центр масс ноды + if (flags[i_node] != level) + continue; + auto &node = nodes[i_node]; + auto &left_node = nodes[node.child_left], &right_node = nodes[node.child_right]; + node.bbox.clear(); + node.bbox.grow(left_node.bbox); + node.bbox.grow(right_node.bbox); + node.mass = left_node.mass + right_node.mass; + node.cmsx = (left_node.cmsx * left_node.mass + right_node.cmsx * right_node.mass) / node.mass; + node.cmsy = (left_node.cmsy * left_node.mass + right_node.cmsy * right_node.mass) / node.mass; + ++n_updated; } -// std::cout << "n updated: " << n_updated << std::endl; + // std::cout << "n updated: " << n_updated << std::endl; // если глубина небольшая, то раньше закончим if (!n_updated) { @@ -1266,13 +1279,12 @@ void buildBBoxes(std::vector &nodes, std::vector &flags, int N, bool } } -void drawLBVH(images::Image &canvas, const std::vector &nodes, int coord_shift) -{ -#pragma omp parallel for +void drawLBVH(images::Image &canvas, const std::vector &nodes, int coord_shift) { + #pragma omp parallel for for (int y = 0; y < (int) canvas.height; ++y) { for (int x = 0; x < (int) canvas.width; ++x) { - Point point{x+coord_shift, y+coord_shift}; + Point point{x + coord_shift, y + coord_shift}; int depth = 0; @@ -1320,11 +1332,10 @@ void drawLBVH(images::Image &canvas, const std::vector &nod } } -void checkLBVHInvariants(const std::vector &nodes, int N) -{ +void checkLBVHInvariants(const std::vector &nodes, int N) { // проверим количество нод в дереве - if (nodes.size() != N-1 + N /*N+1 inner nodes + N leaves*/) { - throw std::runtime_error("4923942304203423: " + std::to_string(nodes.size()) + " vs " + std::to_string(N-1)); + if (nodes.size() != N - 1 + N /*N+1 inner nodes + N leaves*/) { + throw std::runtime_error("4923942304203423: " + std::to_string(nodes.size()) + " vs " + std::to_string(N - 1)); } // у каждой ноды либо нет ни одного ребенка, тогда она лист, либо есть два ребенка @@ -1346,8 +1357,10 @@ void checkLBVHInvariants(const std::vector &nodes, int N) ++used[i_node]; ++total_visited; - if (nodes[i_node].hasLeftChild()) stack.push_back(nodes[i_node].child_left); - if (nodes[i_node].hasRightChild()) stack.push_back(nodes[i_node].child_right); + if (nodes[i_node].hasLeftChild()) + stack.push_back(nodes[i_node].child_left); + if (nodes[i_node].hasRightChild()) + stack.push_back(nodes[i_node].child_right); } // стек не пустой -> есть циклы if (!stack.empty()) { @@ -1361,24 +1374,23 @@ void checkLBVHInvariants(const std::vector &nodes, int N) } } -void nbody(bool interactive, bool evaluate_precision, int nbody_impl_index) -{ +void nbody(bool interactive, bool evaluate_precision, int nbody_impl_index) { // State initial_state = makeRandomState(N, -200, 200, -200, 200); // State initial_state = makeCircularState(); -#if NBODY_INITIAL_STATE_COMPLEXITY == 2 + #if NBODY_INITIAL_STATE_COMPLEXITY == 2 // конфигурация из гифки State initial_state = makeGalacticState(100000, 10, 300); const int canvas_size = 1500; -#elif NBODY_INITIAL_STATE_COMPLEXITY == 1 + #elif NBODY_INITIAL_STATE_COMPLEXITY == 1 // конфигурация полегче State initial_state = makeGalacticState(10000, 10, 100); const int canvas_size = 500; -#else + #else // вообще легкая конфигурация жесть State initial_state = makeGalacticState(1000, 5, 40); const int canvas_size = 150; -#endif + #endif const int N = initial_state.pxs.size(); @@ -1434,17 +1446,18 @@ void nbody(bool interactive, bool evaluate_precision, int nbody_impl_index) int i_frame = 0; timer tm_framerate; - interactive_callback_t interactive_callback = [&](const std::vector &pxs, const std::vector &pys, const std::vector &nodes) -> void { + interactive_callback_t interactive_callback = [&](const std::vector &pxs, const std::vector &pys, + const std::vector &nodes) -> void { canvas.fill(zero); if (nodes.size()) { drawLBVH(canvas, nodes, initial_state.coord_shift - canvas_size / 2); } -#pragma omp parallel for + #pragma omp parallel for for (int i = 0; i < N; ++i) { - float x = state.pxs[i] + canvas.width / 2 - initial_state.coord_shift; - float y = state.pys[i] + canvas.height / 2 - initial_state.coord_shift; + float x = state.pxs[i] + canvas.width / 2.0f - initial_state.coord_shift; + float y = state.pys[i] + canvas.height / 2.0f - initial_state.coord_shift; if (x < 0 || x >= canvas.width || y < 0 || y >= canvas.height) { continue; @@ -1452,7 +1465,7 @@ void nbody(bool interactive, bool evaluate_precision, int nbody_impl_index) Color c = colors[i % colors.size()]; -#pragma omp critical + #pragma omp critical { canvas(y, x, 0) = std::min(canvas(y, x, 0) + c.r * 0.2, 255); @@ -1463,7 +1476,8 @@ void nbody(bool interactive, bool evaluate_precision, int nbody_impl_index) int step = 100; if (interactive && ++i_frame % step == 0) { - std::cout << "simulated " << i_frame << " frames, N: " << N << ", framerate: " << (step / tm_framerate.elapsed()) << " fps, method: " << nbody_impl_name << std::endl; + std::cout << "simulated " << i_frame << " frames, N: " << N << ", framerate: " << ( + step / tm_framerate.elapsed()) << " fps, method: " << nbody_impl_name << std::endl; tm_framerate.restart(); } @@ -1487,7 +1501,8 @@ void nbody(bool interactive, bool evaluate_precision, int nbody_impl_index) nbody_implementation(delta_state, state, N, NT, interactive ? &interactive_callback : nullptr); - std::cout << "simulated " << NT << " frames, N: " << N << ", framerate: " << (NT / tm_framerate.elapsed()) << " fps, method: " << nbody_impl_name << std::endl; + std::cout << "simulated " << NT << " frames, N: " << N << ", framerate: " << (NT / tm_framerate.elapsed()) << + " fps, method: " << nbody_impl_name << std::endl; if (interactive) { return; @@ -1500,16 +1515,17 @@ void nbody(bool interactive, bool evaluate_precision, int nbody_impl_index) int NT_test = std::min(20, NT); nbody_cpu(delta_state_tmp, state_tmp, N, NT_test, nullptr); for (int t = 0; t < NT_test; ++t) { - float * dvx = &delta_state.dvx2d[t * N]; - float * dvy = &delta_state.dvy2d[t * N]; + float *dvx = &delta_state.dvx2d[t * N]; + float *dvy = &delta_state.dvy2d[t * N]; - float * dvx_tmp = &delta_state_tmp.dvx2d[t * N]; - float * dvy_tmp = &delta_state_tmp.dvy2d[t * N]; + float *dvx_tmp = &delta_state_tmp.dvx2d[t * N]; + float *dvy_tmp = &delta_state_tmp.dvy2d[t * N]; int n_good = 0; for (int i = 0; i < N; ++i) { double err = 0.1 * std::abs(dvx_tmp[i]); - if (std::abs(dvx[i] - dvx_tmp[i]) < err) n_good++; + if (std::abs(dvx[i] - dvx_tmp[i]) < err) + n_good++; } EXPECT_GE(n_good, 0.9 * N); } @@ -1519,8 +1535,8 @@ void nbody(bool interactive, bool evaluate_precision, int nbody_impl_index) state = initial_state; for (int t = 0; t < NT; ++t) { - float * dvx = &delta_state.dvx2d[t * N]; - float * dvy = &delta_state.dvy2d[t * N]; + float *dvx = &delta_state.dvx2d[t * N]; + float *dvy = &delta_state.dvy2d[t * N]; for (int i = 0; i < N; ++i) { state.vxs[i] += dvx[i]; @@ -1533,12 +1549,12 @@ void nbody(bool interactive, bool evaluate_precision, int nbody_impl_index) } } -void checkTreesEqual(const std::vector &nodes_recursive, const std::vector &nodes, const Node &root_recursive, const Node &root) -{ +void checkTreesEqual(const std::vector &nodes_recursive, const std::vector &nodes, + const Node &root_recursive, const Node &root) { EXPECT_EQ(root_recursive.bbox, root.bbox); - EXPECT_EQ(root_recursive.mass, root.mass); - EXPECT_EQ(root_recursive.cmsx, root.cmsx); - EXPECT_EQ(root_recursive.cmsy, root.cmsy); + EXPECT_NEAR(root_recursive.mass, root.mass, 0.0001); + EXPECT_NEAR(root_recursive.cmsx, root.cmsx, 0.0001); + EXPECT_NEAR(root_recursive.cmsy, root.cmsy, 0.0001); EXPECT_EQ(root_recursive.hasLeftChild(), root.hasLeftChild()); EXPECT_EQ(root_recursive.hasRightChild(), root.hasRightChild()); @@ -1550,8 +1566,7 @@ void checkTreesEqual(const std::vector &nodes_recursive, const std::vector } } -TEST (LBVH, CPU) -{ +TEST(LBVH, CPU) { if (!ENABLE_TESTING) return; @@ -1572,7 +1587,8 @@ TEST (LBVH, CPU) int N = 10000; std::vector points; std::vector codes; - points.reserve(N); codes.reserve(N); + points.reserve(N); + codes.reserve(N); for (int i = 0; i < N; ++i) { // circle sampling @@ -1607,7 +1623,7 @@ TEST (LBVH, CPU) { std::vector nodes_recursive; - buildLBVHRecursive(nodes_recursive, codes, points, 0, N, NBITS-1); + buildLBVHRecursive(nodes_recursive, codes, points, 0, N, NBITS - 1); buildBBoxesRecursive(nodes_recursive, nodes_recursive.front()); EXPECT_NO_THROW(checkTreesEqual(nodes_recursive, nodes, nodes_recursive.front(), nodes.front())); } @@ -1620,7 +1636,7 @@ TEST (LBVH, CPU) std::vector buffer; for (int i = 1; i < (int) points.size(); ++i) { buffer.clear(); - bresenham(buffer, points[getIndex(codes[i-1])], points[getIndex(codes[i])]); + bresenham(buffer, points[getIndex(codes[i - 1])], points[getIndex(codes[i])]); for (const auto &[x, y] : buffer) { canvas(y, x, 0) = 255; canvas(y, x, 1) = 255; @@ -1639,8 +1655,7 @@ TEST (LBVH, CPU) } } -TEST (LBVH, GPU) -{ +TEST(LBVH, GPU) { if (!ENABLE_TESTING) return; @@ -1797,7 +1812,7 @@ TEST (LBVH, GPU) N, level); int n_updated; - flags_gpu.readN(&n_updated, 1, N-1); + flags_gpu.readN(&n_updated, 1, N - 1); // std::cout << "n updated: " << n_updated << std::endl; @@ -1888,12 +1903,18 @@ TEST (LBVH, GPU) integrate(i, pxs_cpu, pys_cpu, vxs_cpu, vys_cpu, dvx_cpu.data(), dvy_cpu.data(), 0); double rel_eps_super_good = 1e-3; - if (std::abs(pxs[i] - pxs_cpu[i]) < rel_eps_super_good * std::abs(pxs_cpu[i])) n_super_good_pxs++; - if (std::abs(pys[i] - pys_cpu[i]) < rel_eps_super_good * std::abs(pys_cpu[i])) n_super_good_pys++; - if (std::abs(vxs[i] - vxs_cpu[i]) < rel_eps_super_good * std::abs(vxs_cpu[i])) n_super_good_vxs++; - if (std::abs(vys[i] - vys_cpu[i]) < rel_eps_super_good * std::abs(vys_cpu[i])) n_super_good_vys++; - if (std::abs(dvx[i] - dvx_cpu[i]) < rel_eps_super_good * std::abs(dvx_cpu[i])) n_super_good_dvx++; - if (std::abs(dvy[i] - dvy_cpu[i]) < rel_eps_super_good * std::abs(dvy_cpu[i])) n_super_good_dvy++; + if (std::abs(pxs[i] - pxs_cpu[i]) < rel_eps_super_good * std::abs(pxs_cpu[i])) + n_super_good_pxs++; + if (std::abs(pys[i] - pys_cpu[i]) < rel_eps_super_good * std::abs(pys_cpu[i])) + n_super_good_pys++; + if (std::abs(vxs[i] - vxs_cpu[i]) < rel_eps_super_good * std::abs(vxs_cpu[i])) + n_super_good_vxs++; + if (std::abs(vys[i] - vys_cpu[i]) < rel_eps_super_good * std::abs(vys_cpu[i])) + n_super_good_vys++; + if (std::abs(dvx[i] - dvx_cpu[i]) < rel_eps_super_good * std::abs(dvx_cpu[i])) + n_super_good_dvx++; + if (std::abs(dvy[i] - dvy_cpu[i]) < rel_eps_super_good * std::abs(dvy_cpu[i])) + n_super_good_dvy++; double rel_eps = 0.5; EXPECT_NEAR(pxs[i], pxs_cpu[i], rel_eps * std::abs(pxs_cpu[i])); @@ -1913,8 +1934,7 @@ TEST (LBVH, GPU) } } -TEST (LBVH, Nbody) -{ +TEST(LBVH, Nbody) { if (!ENABLE_TESTING) return; @@ -1925,16 +1945,15 @@ TEST (LBVH, Nbody) bool evaluate_precision = (NBODY_INITIAL_STATE_COMPLEXITY < 2) && EVALUATE_PRECISION; -#if NBODY_INITIAL_STATE_COMPLEXITY < 2 - nbody(false, evaluate_precision, 0); // cpu naive - nbody(false, evaluate_precision, 1); // gpu naive -#endif - nbody(false, evaluate_precision, 2); // cpu lbvh + #if NBODY_INITIAL_STATE_COMPLEXITY < 2 + nbody(false, evaluate_precision, 0);// cpu naive + nbody(false, evaluate_precision, 1);// gpu naive + #endif + nbody(false, evaluate_precision, 2);// cpu lbvh nbody(false, evaluate_precision, 3); // gpu lbvh } -TEST (LBVH, Nbody_meditation) -{ +TEST(LBVH, Nbody_meditation) { if (!ENABLE_TESTING) return; @@ -1946,5 +1965,5 @@ TEST (LBVH, Nbody_meditation) context.init(device.device_id_opencl); context.activate(); - nbody(true, false, 3); // gpu lbvh -} \ No newline at end of file + nbody(true, false, 3);// gpu lbvh +}