#include "all/deryabin_m_hoare_sort_simple_merge/i...
作成日: 2025年11月11日
使用モデル GPT-5 Thinking by Chat01
質問
#include "all/deryabin_m_hoare_sort_simple_merge/include/ops_all.hpp"
#include <oneapi/tbb/parallel_for.h>
#include <oneapi/tbb/parallel_invoke.h>
#include <algorithm>
#include <bit>
#include <boost/mpi/collectives/broadcast.hpp>
#include <cmath>
#include <cstddef>
#include <iterator>
#include <vector>
double deryabin_m_hoare_sort_simple_merge_mpi::PivotCalculation(std::vector<double>::iterator first,
std::vector<double>::iterator last) {
const auto mid = first + ((last - first) >> 1);
double pivot_value = NAN;
if (last - first < 199) {
pivot_value = *mid;
} else {
if (*first < *mid) {
if (*mid < *last) {
pivot_value = *mid;
} else {
pivot_value = std::max(*first, *last);
}
} else {
if (*first < *last) {
pivot_value = *first;
} else {
pivot_value = std::max(*mid, *last);
}
}
}
return pivot_value;
}
void deryabin_m_hoare_sort_simple_merge_mpi::SeqHoaraSort(std::vector<double>::iterator first,
std::vector<double>::iterator last) {
if (first >= last) {
return;
}
const double pivot_value = PivotCalculation(first, last);
auto left = first;
auto right = last;
do {
while (left < last && *left < pivot_value) {
left++;
}
while (right > first && *right > pivot_value) {
right--;
}
if (*left == *right && left != right) {
if (*left < *(left + 1)) {
left++;
} else {
right--;
}
}
std::iter_swap(left, right);
} while (left != right);
SeqHoaraSort(first, right);
SeqHoaraSort(left + 1, last);
}
void deryabin_m_hoare_sort_simple_merge_mpi::HoaraSort(std::vector<double>::iterator first,
std::vector<double>::iterator last) {
if (first >= last) {
return;
}
const double pivot_value = PivotCalculation(first, last);
auto left = first;
auto right = last;
do {
while (left < last && *left < pivot_value) {
left++;
}
while (right > first && *right > pivot_value) {
right--;
}
if (*left == *right && left != right) {
if (*left < *(left + 1)) {
left++;
} else {
right--;
}
}
std::iter_swap(left, right);
} while (left != right);
if (last - first < 199) {
HoaraSort(first, right);
HoaraSort(left + 1, last);
} else {
oneapi::tbb::parallel_invoke(&first, &right { HoaraSort(first, right); },
&left, &last { HoaraSort(left + 1, last); });
}
}
void deryabin_m_hoare_sort_simple_merge_mpi::MergeUnequalTwoParts(std::vector<double>::iterator first,
std::vector<double>::iterator mid,
std::vector<double>::iterator last) {
const auto left_end = std::upper_bound(first, mid, *mid);
const auto right_start = std::upper_bound(mid, last, (mid - 1));
const size_t overlap_len = std::min(std::distance(left_end, mid), std::distance(mid, right_start));
oneapi::tbb::parallel_for(oneapi::tbb::blocked_range<size_t>(0, overlap_len),
[&left_end, &mid](const oneapi::tbb::blocked_range<size_t>& r) {
for (size_t i = r.begin(); i < r.end(); ++i) {
auto left = left_end + static_cast<long>(i);
auto right = mid + static_cast<long>(i);
if (left > right) {
std::iter_swap(left, right);
}
}
});
size_t right_len = right_start - (mid - 1);
size_t left_len = mid - left_end;
if (right_len > left_len + 1) {
size_t delta = right_len - left_len;
auto base = &(right_start - static_cast<long>(delta));
oneapi::tbb::parallel_for(oneapi::tbb::blocked_range<size_t>(0, left_len),
[&base, delta](const oneapi::tbb::blocked_range<size_t>& r) {
for (size_t j = r.begin(); j < r.end(); ++j) {
double current = base - j;
for (size_t i = 0; i < delta - 1; ++i) {
if (current[i] > current[i + 1]) {
std::swap(current[i], current[i + 1]);
}
}
}
});
}
std::inplace_merge(left_end, mid, right_start);
}
bool deryabin_m_hoare_sort_simple_merge_mpi::HoareSortTaskSequential::PreProcessingImpl() {
input_array_A_ = reinterpret_cast<std::vector<double>*>(task_data->inputs[0])[0];
dimension_ = task_data->inputs_count[0];
chunk_count_ = task_data->inputs_count[1];
min_chunk_size_ = dimension_ / chunk_count_;
return true;
}
bool deryabin_m_hoare_sort_simple_merge_mpi::HoareSortTaskSequential::ValidationImpl() {
return task_data->inputs_count[0] > 2 && task_data->inputs_count[1] >= 2 &&
task_data->inputs_count[0] == task_data->outputs_count[0];
}
bool deryabin_m_hoare_sort_simple_merge_mpi::HoareSortTaskSequential::RunImpl() {
size_t count = 0;
while (count != chunk_count_) {
HoaraSort(input_array_A_.begin() + static_cast<long>(count * min_chunk_size_),
input_array_A_.begin() + static_cast<long>((count + 1) * min_chunk_size_) - 1);
count++;
}
size_t chunk_count = chunk_count_;
for (size_t i = 0; i < static_cast<size_t>(std::bit_width(chunk_count_)) - 1; i++) {
for (size_t j = 0; j < chunk_count; j++) {
std::inplace_merge(input_array_A_.begin() + static_cast<long>(j * min_chunk_size_ << (i + 1)),
input_array_A_.begin() + static_cast<long>((j << 1 | 1) * (min_chunk_size_ << i)),
input_array_A_.begin() + static_cast<long>((j + 1) * min_chunk_size_ << (i + 1)));
chunk_count--;
}
}
return true;
}
bool deryabin_m_hoare_sort_simple_merge_mpi::HoareSortTaskSequential::PostProcessingImpl() {
reinterpret_cast<std::vector<double>*>(task_data->outputs[0])[0] = input_array_A_;
return true;
}
bool deryabin_m_hoare_sort_simple_merge_mpi::HoareSortTaskMPI::PreProcessingImpl() {
if (world_.rank() == 0) {
input_array_A_ = reinterpret_cast<std::vector<double>*>(task_data->inputs[0])[0];
dimension_ = task_data->inputs_count[0];
chunk_count_ = world_.size();
min_chunk_size_ = dimension_ / chunk_count_;
rest_ = dimension_ % chunk_count_;
}
boost::mpi::broadcast(world_, dimension_, 0);
if (world_.rank() != 0) {
input_array_A_.resize(dimension_);
}
boost::mpi::broadcast(world_, input_array_A_.data(), static_cast<int>(dimension_), 0);
boost::mpi::broadcast(world_, chunk_count_, 0);
boost::mpi::broadcast(world_, min_chunk_size_, 0);
boost::mpi::broadcast(world_, rest_, 0);
return true;
}
bool deryabin_m_hoare_sort_simple_merge_mpi::HoareSortTaskMPI::ValidationImpl() {
if (world_.rank() == 0) {
return task_data->inputs_count[0] > 2 && task_data->inputs_count[0] == task_data->outputs_count[0];
}
return true;
}
bool deryabin_m_hoare_sort_simple_merge_mpi::HoareSortTaskMPI::HandleEvenCase(size_t world_rank, size_t step) {
if (world_rank == 0) {
return true;
}
size_t block_size = min_chunk_size_ * step;
size_t start_idx = (chunk_count_ - (world_rank + step)) * min_chunk_size_;
if (world_rank == chunk_count_ - step) {
block_size += rest_;
} else {
start_idx += rest_;
}
if (world_rank >= step) {
world_.send(static_cast<int>(world_rank - step), 0, input_array_A_.data() + start_idx,
static_cast<int>(block_size));
} else {
world_.send(0, 0, input_array_A_.data() + start_idx, static_cast<int>(block_size));
}
return true;
}
bool deryabin_m_hoare_sort_simple_merge_mpi::HoareSortTaskMPI::HandleOddCase(size_t world_rank, size_t step) {
const bool special_odd_case = (chunk_count_ & 1) != 0U && world_rank == 0;
size_t block_size = min_chunk_size_ * step;
size_t start_idx = special_odd_case ? (chunk_count_ - (2 * step - 1)) * min_chunk_size_
: (chunk_count_ - (world_rank + 2 * step)) * min_chunk_size_;
const size_t recv_rank = special_odd_case ? step - 1 : world_rank + step;
if (recv_rank == chunk_count_ - step) {
block_size += rest_;
} else {
start_idx += rest_;
}
world_.recv(static_cast<int>(recv_rank), 0, input_array_A_.data() + start_idx, static_cast<int>(block_size));
const auto end_idx = special_odd_case
? input_array_A_.end()
: input_array_A_.begin() + static_cast<long>(start_idx + (block_size * 2) - rest_);
if (end_idx - (input_array_A_.begin() + static_cast<long>(start_idx)) < 200) {
std::inplace_merge(input_array_A_.begin() + static_cast<long>(start_idx),
input_array_A_.begin() + static_cast<long>(start_idx + block_size), end_idx);
} else {
MergeUnequalTwoParts(input_array_A_.begin() + static_cast<long>(start_idx),
input_array_A_.begin() + static_cast<long>(start_idx + block_size), end_idx);
}
return true;
}
bool deryabin_m_hoare_sort_simple_merge_mpi::HoareSortTaskMPI::RunImpl() {
const size_t world_rank = world_.rank();
const auto start_iter =
input_array_A_.begin() + static_cast<long>(((chunk_count_ - world_rank - 1) * min_chunk_size_) +
(world_rank == chunk_count_ - 1 ? 0 : rest_));
const auto end_iter =
input_array_A_.begin() + static_cast<long>(((chunk_count_ - world_rank) * min_chunk_size_) + rest_) - 1;
HoaraSort(start_iter, end_iter);
const size_t iterations = std::bit_width(chunk_count_ - 1);
for (size_t i = 0; i < iterations; ++i) {
const size_t step = 1ULL << i;
if (((chunk_count_ - world_rank) & (step - 1)) == 0U || world_rank == 0) {
const bool is_even = ((chunk_count_ & 1) != 0U) ? (((chunk_count_ - world_rank + step - 1) / step) & 1) != 0U
: (((chunk_count_ - world_rank - 1) / step) & 1) == 0U;
if (is_even) {
HandleEvenCase(world_rank, step);
} else {
HandleOddCase(world_rank, step);
}
}
}
return true;
}
bool deryabin_m_hoare_sort_simple_merge_mpi::HoareSortTaskMPI::PostProcessingImpl() {
if (world_.rank() == 0) {
reinterpret_cast<std::vector<double>*>(task_data->outputs[0])[0] = input_array_A_;
}
return true;
} - внимательно изучи эту реализацию быстрой сортировки - особое внимание обрати на тактику слияния - здесь используется деление на несколько уровней и на каждом следующем уровне происходит слив результатов на сужающийся диапазон процессов - причём обрати внимание что с ростом уровней этот диапазон сужается - также обрати внимание что это работает не только для числа процессов степеней двойки - это работает не только например для 8 процессов - но также для 7 - из-за особой логики распределения ресурсов между процессами и особой логики выбора кому что отправлять на каждом уровне между процессами - ты должен очень внимательно и очень аккуратно применить эту логику рассылок и распределения к моему АГП, но сначала что я хочу - похоже я столкнулся с пределом пропускной способности сети MPI - поэтому я хочу использовать следующий подход - например на начальной стадии поиска (пусть стадия поиска будет определяться через прогресс - например по близости dmax к нулю - в алгоритме уже встречается много раз определения прогресса - можно взять один из них - например какой-нибудь экспоненциальный тип прогресса) и разделить этот прогресс на число уровней в зависимости от числа доступных процессов - например если доступно только 4 процесса - то будет только два уровня - на первом уровне например 3ий отправляет 2ому - 2ой только принимает, 1ый - нулевому - нулевой только принимает, на следующем уровне то что приняли процессы на предыдущем уровне уже начинают рассылать - например в данном примере на предыдущем уровне принимали 2ой и нулевой - у них самые актуальные данные - и на самом последнем уровне пусть происходит взаимный обмен между этими двумя процессами с самыми актуальными данными - чтобы выявить действительно лучшие данные и пусть финальный взаимный обмен лучшими значениями происходит в конце AGP_Manip2D а обмены до этого в agp_run_branch_mpi, это нужно применить к коду: // agp_algo.cpp
// Полностью обновлённый файл DLL с учётом всех требований.
// ВАЖНО: собирается как C++17/20, /arch:AVX2, /Oi, /O2, /fp:fast желательны.
#include "pch.h"
// ==================== ВАШИ МАКРО И ФУНКЦИИ (без изменений интерфейсов), мелкие доработки дальше ====================
#define XOR_RAND(state, result_var)
do {
uint32_t s = (state);
s ^= s << 13;
s ^= s >> 17;
s ^= s << 5;
(state) = s;
float tmp = (float)((double)(s) * (1.0/4294967296.0));
result_var = tmp;
} while (0)
#define XOR_RAND_GRSH(state, result_var)
do {
uint32_t s = (state);
s ^= s << 13;
s ^= s >> 17;
s ^= s << 5;
(state) = s;
result_var = fmaf((float)(int32_t)s, 0x1.0p-31f, -1.0f);
} while (0)
#define FABE13_COS(x, result_var)
do {
const float ax = fabsf(x);
float r = fmodf(ax, 6.28318530718f);
if (r > 3.14159265359f)
r = 6.28318530718f - r;
if (r < 1.57079632679f) {
const float t2 = r * r;
const float t4 = t2 * t2;
result_var = fmaf(t4, fmaf(t2, -0.0013888889f, 0.0416666667f), fmaf(t2, -0.5f, 1.0f));
} else {
r = 3.14159265359f - r;
const float t2 = r * r;
const float t4 = t2 * t2;
result_var = -fmaf(t4, fmaf(t2, -0.0013888889f, 0.0416666667f), fmaf(t2, -0.5f, 1.0f));
}
} while (0)
#define FABE13_SIN(x, result_var)
do {
const float x = (x);
const float ax = fabsf(x);
float r = fmodf(ax, 6.28318530718f);
bool sfl = r > 3.14159265359f;
if (sfl)
r = 6.28318530718f - r;
bool cfl = r > 1.57079632679f;
if (cfl)
r = 3.14159265359f - r;
const float t2 = r * r;
float _s = fmaf(t2, fmaf(t2, fmaf(t2, -0.0001984127f, 0.0083333333f), -0.16666666f), 1.0f) * r;
result_var = ((x < 0.0f) ^ sfl) ? -_s : _s;
} while (0)
#define FABE13_SINCOS(in, sin_out, cos_out, n)
do {
int i = 0;
const int limit = (n) & ~7;
if ((n) >= 8) {
static __declspec(align(32)) const __m256 VEC_TWOPI = _mm256_set1_ps(6.28318530718f);
static __declspec(align(32)) const __m256 VEC_PI = _mm256_set1_ps(3.14159265359f);
static __declspec(align(32)) const __m256 VEC_PI_2 = _mm256_set1_ps(1.57079632679f);
static __declspec(align(32)) const __m256 INV_TWOPI = _mm256_set1_ps(0.15915494309189535f);
static __declspec(align(32)) const __m256 BIAS = _mm256_set1_ps(12582912.0f);
static __declspec(align(32)) const __m256 VEC_COS_P5 = _mm256_set1_ps(-0.0013888889f);
static __declspec(align(32)) const __m256 VEC_COS_P3 = _mm256_set1_ps(0.0416666667f);
static __declspec(align(32)) const __m256 VEC_COS_P1 = _mm256_set1_ps(-0.5f);
static __declspec(align(32)) const __m256 VEC_COS_P0 = _mm256_set1_ps(1.0f);
static __declspec(align(32)) const __m256 VEC_SIN_P5 = _mm256_set1_ps(-0.0001984127f);
static __declspec(align(32)) const __m256 VEC_SIN_P3 = _mm256_set1_ps(0.0083333333f);
static __declspec(align(32)) const __m256 VEC_SIN_P1 = _mm256_set1_ps(-0.16666666f);
static __declspec(align(32)) const __m256 VEC_SIN_P0 = _mm256_set1_ps(1.0f);
static __declspec(align(32)) const __m256 VEC_ZERO = _mm256_setzero_ps();
while (i < limit) {
const __m256 vx = _mm256_load_ps(&(in)[i]);
const __m256 vax = _mm256_andnot_ps(_mm256_set1_ps(-0.0f), vx);
__m256 q = _mm256_fmadd_ps(vax, INV_TWOPI, BIAS);
q = _mm256_sub_ps(q, BIAS);
const __m256 r = _mm256_fnmadd_ps(VEC_TWOPI, q, vax);
const __m256 r1 = _mm256_min_ps(r, _mm256_sub_ps(VEC_TWOPI, r));
const __m256 r2 = _mm256_min_ps(r1, _mm256_sub_ps(VEC_PI, r1));
const __m256 t2 = _mm256_mul_ps(r2, r2);
const __m256 cosv = _mm256_fmadd_ps(t2, _mm256_fmadd_ps(t2, _mm256_fmadd_ps(t2, VEC_COS_P5, VEC_COS_P3), VEC_COS_P1), VEC_COS_P0);
const __m256 sinv = _mm256_mul_ps(_mm256_fmadd_ps(t2, _mm256_fmadd_ps(t2, _mm256_fmadd_ps(t2, VEC_SIN_P5, VEC_SIN_P3), VEC_SIN_P1), VEC_SIN_P0), r2);
const __m256 cflip = _mm256_cmp_ps(r1, VEC_PI_2, _CMP_GT_OQ);
const __m256 sflip = _mm256_xor_ps(_mm256_cmp_ps(vx, VEC_ZERO, _CMP_LT_OQ), _mm256_cmp_ps(r, VEC_PI, _CMP_GT_OQ));
_mm256_store_ps(&(cos_out)[i], _mm256_blendv_ps(cosv, _mm256_sub_ps(VEC_ZERO, cosv), cflip));
_mm256_store_ps(&(sin_out)[i], _mm256_blendv_ps(sinv, _mm256_sub_ps(VEC_ZERO, sinv), sflip));
i += 8;
}
}
while (i < (n)) {
const float x = (in)[i];
const float ax = fabsf(x);
float q = fmaf(ax, 0.15915494309189535f, 12582912.0f);
q -= 12582912.0f;
float r = fmaf(-6.28318530718f, q, ax);
const bool sflip = r > 3.14159265359f;
if (sflip)
r = 6.28318530718f - r;
const bool cflip = r > 1.57079632679f;
if (cflip)
r = 3.14159265359f - r;
const float t2 = r * r;
const float c = fmaf(t2, fmaf(t2, fmaf(t2, -0.0013888889f, 0.0416666667f), -0.5f), 1.0f);
const float s = fmaf(t2, fmaf(t2, fmaf(t2, -0.0001984127f, 0.0083333333f), -0.16666666f), 1.0f) * r;
(cos_out)[i] = cflip ? -c : c;
(sin_out)[i] = ((x < 0.0f) ^ sflip) ? -s : s;
++i;
}
} while (0)
// ---------------- Peano / Morton таблички и MPI-глобали (без изменения интерфейса) ----------------
enum List : uint8_t { Top = 0b00u, Down = 0b01u, Left = 0b10u, Right = 0b11u };
__declspec(align(4)) struct Step final { const uint8_t next, dx, dy; };
__declspec(align(4)) struct InvStep final { const uint8_t q, next; };
__declspec(align(64)) static const Step g_step_tbl[4][4] = {
{ { Right,0u,0u }, { Top,0u,1u }, { Top,1u,1u }, { Left,1u,0u } },
{ { Left,1u,1u }, { Down,1u,0u }, { Down,0u,0u }, { Right,0u,1u } },
{ { Down,1u,1u }, { Left,0u,1u }, { Left,0u,0u }, { Top,1u,0u } },
{ { Top,0u,0u }, { Right,1u,0u }, { Right,1u,1u }, { Down,0u,1u } }
};
__declspec(align(64)) static const InvStep g_inv_tbl[4][4] = {
{ { 0u,Right }, { 1u,Top }, { 3u,Left }, { 2u,Top } },
{ { 2u,Down }, { 3u,Right }, { 1u,Down }, { 0u,Left } },
{ { 2u,Left }, { 1u,Left }, { 3u,Top }, { 0u,Down } },
{ { 0u,Top }, { 3u,Down }, { 1u,Right }, { 2u,Right } }
};
static const boost::mpi::environment* g_env;
static const boost::mpi::communicator* g_world;
// ---------------- Служебные структуры сообщений и slab-аллокатор ----------------
__declspec(align(16)) struct CrossMsg final {
float s_x1, s_x2, e_x1, e_x2, Rtop;
template<typename Archive> __declspec(noalias) __forceinline void serialize(Archive& ar, const unsigned int) noexcept { ar& s_x1& s_x2& e_x1& e_x2& Rtop; }
};
__declspec(align(16)) struct CtrlMsg final {
uint8_t kind; CrossMsg xchg;
template<typename Archive> __declspec(noalias) __forceinline void serialize(Archive& ar, const unsigned int) noexcept { ar& kind& xchg; }
};
__declspec(align(16)) struct Slab final {
char* const base; char* current; char* const end;
__forceinline Slab(void* const memory, const size_t usable) : base((char*)memory), current(base), end(base + (usable & ~(size_t)63u)) {}
};
static tbb::enumerable_thread_specific<Slab*> tls( noexcept {
void* memory = _aligned_malloc(16777216u, 64u);
Slab* slab = (Slab*)_aligned_malloc(32u, 64u);
new (slab) Slab(memory, 16777216u);
char* p = slab->base;
while (p < slab->end) { *p = 0; p += 4096u; }
return slab;
});
// ---------------- Параметры активной 2D-карты (как было) ----------------
__declspec(align(16)) struct Peano2DMap final {
const int levels; const float a, b, c, d; const float lenx, leny; const float inv_lenx; const uint32_t scale; const uint8_t start;
__forceinline Peano2DMap(int L, float _a, float _b, float _c, float _d, uint8_t st)
: levels(L), a(_a), b(_b), c(_c), d(_d),
lenx(_b - _a), leny(_d - _c), inv_lenx(1.0f / (_b - _a)),
scale((uint32_t)1u << (L << 1)), start(st) {
}
};
static Peano2DMap gActiveMap(0, 0, 0, 0, 0, 0);
// ---------------- Интервалы и сравнения (как было) ----------------
__declspec(align(64)) struct Interval1D final {
const float x1, x2, y1, y2, delta_y, ordinate_factor, N_factor, quadratic_term, M; float R;
__forceinline void* operator new(size_t) noexcept { Slab* s = tls.local(); char* r = s->current; s->current += 64u; return r; }
__forceinline Interval1D(float _x1, float _x2, float _y1, float _y2, float _N) noexcept
: x1(_x1), x2(_x2), y1(_y1), y2(_y2), delta_y(_y2 - _y1), ordinate_factor(-(y1 + y2) * 2.0f),
N_factor(_N == 1.0f ? _x2 - _x1 : sqrtf(_x2 - _x1)),
quadratic_term((1.0f / N_factor)* delta_y* delta_y),
M((1.0f / N_factor)* fabsf(delta_y)) {
}
__forceinline void ChangeCharacteristic(float _m) noexcept { R = fmaf(1.0f / _m, quadratic_term, fmaf(_m, N_factor, ordinate_factor)); }
};
__declspec(align(64)) struct IntervalND final {
const float x1, x2, y1, y2, delta_y, ordinate_factor;
float N_factor, quadratic_term, M, R;
uint64_t i1, i2; float diam; int span_level;
__forceinline void* operator new(size_t) noexcept { Slab* s = tls.local(); char* r = s->current; s->current += 64u; return r; }
__forceinline IntervalND(float _x1, float _x2, float _y1, float _y2) noexcept
: x1(_x1), x2(_x2), y1(_y1), y2(_y2), delta_y(_y2 - _y1), ordinate_factor(-(y1 + y2) * 2.0f),
N_factor(0), quadratic_term(0), M(0), R(0), i1(0), i2(0), diam(0), span_level(0) {
}
__forceinline void compute_span_level(const struct MortonND& map) noexcept; // объявим ниже
__forceinline void set_metric(float d_alpha) noexcept { N_factor = d_alpha; quadratic_term = (1.0f / N_factor) * delta_y * delta_y; M = (1.0f / N_factor) * fabsf(delta_y); }
__forceinline void ChangeCharacteristic(float _m) noexcept { R = fmaf(1.0f / _m, quadratic_term, fmaf(_m, N_factor, ordinate_factor)); }
};
__forceinline bool ComparePtr1D(const Interval1D* a, const Interval1D* b) noexcept { return a->R < b->R; }
__forceinline bool ComparePtrND(const IntervalND* a, const IntervalND* b) noexcept { return a->R < b->R; }
// ---------------- Массовые пересчёты R (как было) ----------------
__forceinline void RecomputeR_ConstM_AVX2_1D(Interval1D* const* arr, size_t n, float m) {
const __m256 vm = _mm256_set1_ps(m);
__m256 vinvm = _mm256_rcp_ps(vm);
vinvm = _mm256_mul_ps(vinvm, _mm256_fnmadd_ps(vm, vinvm, _mm256_set1_ps(2.0f)));
size_t i = 0, limit = n & ~7ull;
alignas(32) float q[8], nf[8], od[8], out[8];
for (; i < limit; i += 8) {
for (int k = 0; k < 8; ++k) { const Interval1D* p = arr[i + k]; q[k] = p->quadratic_term; nf[k] = p->N_factor; od[k] = p->ordinate_factor; }
const __m256 vq = _mm256_load_ps(q), vnf = _mm256_load_ps(nf), vod = _mm256_load_ps(od);
const __m256 t = _mm256_fmadd_ps(vm, vnf, vod);
const __m256 res = _mm256_fmadd_ps(vq, vinvm, t);
_mm256_store_ps(out, res);
for (int k = 0; k < 8; ++k) arr[i + k]->R = out[k];
}
for (; i < n; ++i) arr[i]->ChangeCharacteristic(m);
}
__forceinline void RecomputeR_AffineM_AVX2_1D(Interval1D* const* arr, size_t n, float GF, float alpha) {
const __m256 vGF = _mm256_set1_ps(GF), va = _mm256_set1_ps(alpha);
size_t i = 0, limit = n & ~7ull;
alignas(32) float ln[8], Mv[8], q[8], nf[8], od[8], out[8];
for (; i < limit; i += 8) {
for (int k = 0; k < 8; ++k) {
const Interval1D* p = arr[i + k];
ln[k] = p->x2 - p->x1; Mv[k] = p->M; q[k] = p->quadratic_term; nf[k] = p->N_factor; od[k] = p->ordinate_factor;
}
const __m256 vln = _mm256_load_ps(ln), vM = _mm256_load_ps(Mv), vq = _mm256_load_ps(q), vnf = _mm256_load_ps(nf), vod = _mm256_load_ps(od);
const __m256 vm = _mm256_fmadd_ps(vGF, vln, _mm256_mul_ps(va, vM));
__m256 vinvm = _mm256_rcp_ps(vm);
vinvm = _mm256_mul_ps(vinvm, _mm256_fnmadd_ps(vm, vinvm, _mm256_set1_ps(2.0f)));
const __m256 t = _mm256_fmadd_ps(vm, vnf, vod);
const __m256 res = _mm256_fmadd_ps(vq, vinvm, t);
_mm256_store_ps(out, res);
for (int k = 0; k < 8; ++k) arr[i + k]->R = out[k];
}
for (; i < n; ++i) {
const Interval1D* p = arr[i];
const float mi = fmaf(GF, (p->x2 - p->x1), p->M * alpha);
arr[i]->R = fmaf(1.0f / mi, p->quadratic_term, fmaf(mi, p->N_factor, p->ordinate_factor));
}
}
__forceinline void RecomputeR_ConstM_AVX2_ND(IntervalND* const* arr, size_t n, float m) {
const __m256 vm = _mm256_set1_ps(m);
__m256 vinvm = _mm256_rcp_ps(vm);
vinvm = _mm256_mul_ps(vinvm, _mm256_fnmadd_ps(vm, vinvm, _mm256_set1_ps(2.0f)));
size_t i = 0, limit = n & ~7ull;
alignas(32) float q[8], nf[8], od[8], out[8];
for (; i < limit; i += 8) {
for (int k = 0; k < 8; ++k) { const IntervalND* p = arr[i + k]; q[k] = p->quadratic_term; nf[k] = p->N_factor; od[k] = p->ordinate_factor; }
const __m256 vq = _mm256_load_ps(q), vnf = _mm256_load_ps(nf), vod = _mm256_load_ps(od);
const __m256 t = _mm256_fmadd_ps(vm, vnf, vod);
const __m256 res = _mm256_fmadd_ps(vq, vinvm, t);
_mm256_store_ps(out, res);
for (int k = 0; k < 8; ++k) arr[i + k]->R = out[k];
}
for (; i < n; ++i) arr[i]->ChangeCharacteristic(m);
}
__forceinline void RecomputeR_AffineM_AVX2_ND(IntervalND* const* arr, size_t n, float GF, float alpha) {
const __m256 vGF = _mm256_set1_ps(GF), va = _mm256_set1_ps(alpha);
size_t i = 0, limit = n & ~7ull;
alignas(32) float ln[8], Mv[8], q[8], nf[8], od[8], out[8];
for (; i < limit; i += 8) {
for (int k = 0; k < 8; ++k) {
const IntervalND* p = arr[i + k];
ln[k] = p->x2 - p->x1; Mv[k] = p->M; q[k] = p->quadratic_term; nf[k] = p->N_factor; od[k] = p->ordinate_factor;
}
const __m256 vln = _mm256_load_ps(ln), vM = _mm256_load_ps(Mv), vq = _mm256_load_ps(q), vnf = _mm256_load_ps(nf), vod = _mm256_load_ps(od);
const __m256 vm = _mm256_fmadd_ps(vGF, vln, _mm256_mul_ps(va, vM));
__m256 vinvm = _mm256_rcp_ps(vm);
vinvm = _mm256_mul_ps(vinvm, _mm256_fnmadd_ps(vm, vinvm, _mm256_set1_ps(2.0f)));
const __m256 t = _mm256_fmadd_ps(vm, vnf, vod);
const __m256 res = _mm256_fmadd_ps(vq, vinvm, t);
_mm256_store_ps(out, res);
for (int k = 0; k < 8; ++k) arr[i + k]->R = out[k];
}
for (; i < n; ++i) {
const IntervalND* p = arr[i];
const float mi = fmaf(GF, (p->x2 - p->x1), p->M * alpha);
arr[i]->R = fmaf(1.0f / mi, p->quadratic_term, fmaf(mi, p->N_factor, p->ordinate_factor));
}
}
// ---------------- Быстрые степени (как было) ----------------
__forceinline float fast_pow_int(float v, int n) {
float r;
switch (n) {
case 3: { float v2 = v * v; r = v2 * v; } break;
case 4: { float v2 = v * v; r = v2 * v2; } break;
case 5: { float v2 = v * v; r = v2 * v2 * v; } break;
case 6: { float v2 = v * v; float v4 = v2 * v2; r = v4 * v2; } break;
case 7: { float v2 = v * v; float v4 = v2 * v2; r = v4 * v2 * v; } break;
case 8: { float v2 = v * v; float v4 = v2 * v2; r = v4 * v4; } break;
case 9: { float v3 = v * v * v; float v6 = v3 * v3; r = v6 * v3; } break;
case 10: { float v2 = v * v; float v4 = v2 * v2; float v8 = v4 * v4; r = v8 * v2; } break;
case 11: { float v2 = v * v; float v4 = v2 * v2; float v8 = v4 * v4; r = v8 * v2 * v; } break;
case 12: { float v3 = v * v * v; float v6 = v3 * v3; r = v6 * v6; } break;
case 13: { float v3 = v * v * v; float v6 = v3 * v3; r = v6 * v6 * v; } break;
case 14: { float v7 = v * v * v * v * v * v * v; r = v7 * v7; } break;
case 15: { float v7 = v * v * v * v * v * v * v; r = v7 * v7 * v; } break;
default: { float v2 = v * v; float v4 = v2 * v2; float v8 = v4 * v4; r = v8 * v8; }
}
return r;
}
// ---------------- Функция шага (как было) ----------------
__forceinline float Shag(float _m, float x1, float x2, float y1, float y2, float _N, float _r) {
const float diff = y2 - y1;
const float sign_mult = _mm_cvtss_f32(_mm_castsi128_ps(_mm_set1_epi32(0x3F800000u | (((((uint32_t)&diff)) & 0x80000000u) ^ 0x80000000u))));
if (_N == 1.0f)
return fmaf(-(1.0f / _m), diff, x1 + x2) * 0.5f;
if (_N == 2.0f)
return fmaf(sign_mult / (_m * _m), diff * diff * _r, x1 + x2) * 0.5f;
return fmaf((1.0f / fast_pow_int(_m, _N)) * sign_mult, fast_pow_int(fabsf(diff), _N) * _r, x1 + x2) * 0.5f;
}
// ---------------- Morton-кэш ----------------
struct MortonCachePerRank {
std::vector<int> permCache;
std::vector<uint64_t> invMaskCache;
uint32_t baseSeed;
};
static MortonCachePerRank g_mc;
static __forceinline uint64_t gray_encode(uint64_t x) noexcept { return x ^ (x >> 1); }
static __forceinline uint64_t gray_decode(uint64_t g) noexcept { g ^= g >> 32; g ^= g >> 16; g ^= g >> 8; g ^= g >> 4; g ^= g >> 2; g ^= g >> 1; return g; }
// ---------------- MortonND: ВАЖНО — поддержка levels*dim > 64 через «чанки» ----------------
struct MortonND final {
int dim, levels; // полные уровни
int eff_levels; // уровни в одном 64-битном чанке (<= floor(63/dim))
int extra_levels; // остаток уровней за пределами первого чанка (>=0)
int chunks; // количество чанков
std::vector<int> chunk_bits; // кол-во уровней в каждом чанке
std::vector<uint64_t> chunk_bases; // 2^(dim*chunk_bits[c]) на каждый чанк
uint64_t scale; // БАЗА для индекса старшего чанка (совместима со старыми местами использования)
std::vector<float> low, high, step, invStep, baseOff;
std::vector<int> perm;
std::vector<uint64_t> invMask; // длиной levels (используем кусочно)
std::vector<uint64_t> pextMask; // маски для СТАРШЕГО чанка (для эвристик/диаметра)
std::vector<uint64_t> pextMaskChunks; // маски (chunk,dim) для map01ToPoint
float invScaleLevel;
bool use_gray;
textstatic inline uint64_t make_mask(int dim, int Lc, int d) { uint64_t m = 0ull, bitpos = (uint64_t)d; for (int b = 0; b < Lc; ++b) { m |= 1ull << bitpos; bitpos += (uint64_t)dim; } return m; } __forceinline MortonND(int D, int L, const float* lows, const float* highs, const MortonCachePerRank& mc) : dim(D), levels(L), eff_levels((std::max)(1, (int)(63 / (D ? D : 1)))), extra_levels((L > eff_levels) ? (L - eff_levels) : 0), chunks((extra_levels > 0) ? (1 + (extra_levels + eff_levels - 1) / eff_levels) : 1), low(lows, lows + D), high(highs, highs + D), step(D, 0.0f), invStep(D, 0.0f), baseOff(D, 0.0f), perm(mc.permCache.begin(), mc.permCache.begin() + D), invMask(mc.invMaskCache.begin(), mc.invMaskCache.begin() + D), invScaleLevel(1.0f / (float)((uint64_t)1 << L)), use_gray(true) { // Размер шага и базовое смещение — по ПОЛНОЙ глубине: for (int d = 0; d < dim; ++d) { float rng = high[d] - low[d]; float st = rng * invScaleLevel; step[d] = st; invStep[d] = 1.0f / st; baseOff[d] = fmaf(0.5f, st, low[d]); // центр ячейки на полном уровне } // Инициализация чанков chunk_bits.resize(chunks); pextMaskChunks.resize((size_t)chunks * (size_t)dim); chunk_bases.resize(chunks); int remaining = levels; for (int c = 0; c < chunks; ++c) { int Lc = (c == 0) ? (std::min)(eff_levels, remaining) : (std::min)(eff_levels, remaining); chunk_bits[c] = Lc; remaining -= Lc; uint64_t baseC = (uint64_t)1 << (dim * Lc); chunk_bases[c] = baseC; for (int d = 0; d < dim; ++d) pextMaskChunks[(size_t)c * (size_t)dim + (size_t)d] = make_mask(dim, Lc, d); } // Маски для старшего чанка (для диаметра/эвристик) pextMask.resize(dim); for (int d = 0; d < dim; ++d) pextMask[d] = make_mask(dim, chunk_bits[0], d); // Маски инверсий — используем только нужные биты для каждого чанка // (оставляем как есть; сдвиги/маски применим динамически) // БАЗОВАЯ «scale» — база старшего чанка: 2^(dim*eff_levels_actual_for_chunk0) scale = (uint64_t)1 << (dim * chunk_bits[0]); } // Диаметр между двумя индексами старшего чанка, с учётом «скрытых» уровней: добавим extra_levels в nfree как константу __forceinline float block_diameter(uint64_t i1, uint64_t i2) const noexcept { if (i1 > i2) std::swap(i1, i2); float s2 = 0.0f; for (int d = 0; d < dim; ++d) { const int pd = perm[d]; const uint64_t varying = (i1 ^ i2) & pextMask[d]; const int nfree_hi = _mm_popcnt_u64(varying); const int nfree_total = nfree_hi + (levels - chunk_bits[0]); const float range = step[pd] * (ldexpf(1.0f, nfree_total) - 1.0f); s2 = fmaf(range, range, s2); } return sqrtf(s2); } // Преобразование t∈[0,1) → точка, используя все чанки (точность сохраняется) __forceinline void map01ToPoint(float t, float* __restrict out) const noexcept { if (t <= 0.0f) t = 0.0f; else if (t >= 1.0f) t = 0x1.fffffep-1f; // накапливаем биты для каждой размерности uint64_t accBits[32] = { 0ull }; int accShifted[32] = { 0 }; double u = (double)t; // используем double для дробной части for (int c = 0; c < chunks; ++c) { const int Lc = chunk_bits[c]; const uint64_t baseC = chunk_bases[c]; double scaled = u * (double)baseC; uint64_t idxc = (scaled >= (double)baseC) ? (baseC - 1ull) : (uint64_t)scaled; u = scaled - (double)idxc; // дробная часть для следующего чанка if (use_gray) idxc = gray_encode(idxc); // какой набор бит инверсии используется для этого чанка: int shift_from_top = 0; for (int k = 0; k <= c; ++k) shift_from_top += chunk_bits[k]; int inv_shift = levels - shift_from_top; // смещение начала подпоследовательности в общем инверс-маске for (int d = 0; d < dim; ++d) { int pd = perm[d]; uint64_t mask = pextMaskChunks[(size_t)c * (size_t)dim + (size_t)d]; uint64_t bits = _pext_u64(idxc, mask); if (inv_shift >= 0) { uint64_t invMaskSegment = 0ull; if (chunk_bits[c] < 63) { uint64_t take = ((uint64_t)1 << chunk_bits[c]) - 1ull; invMaskSegment = (invMask[pd] >> inv_shift) & take; } bits ^= invMaskSegment; } accBits[pd] = (accBits[pd] << Lc) | bits; accShifted[pd] += Lc; } } for (int d = 0; d < dim; ++d) out[d] = fmaf(step[d], (float)accBits[d], baseOff[d]); } // t для лучшего-приблизительного позиционирования (используем только старший чанк — быстро ради эвристик) __forceinline float pointToT(const float* __restrict q) const noexcept { const int bitsFull = levels; const int bitsCoarse = chunk_bits[0]; uint64_t idx0 = 0ull; for (int d = 0; d < dim; ++d) { int pd = perm[d]; float v = (q[pd] - baseOff[pd]) * invStep[pd]; int64_t cell = (int64_t)_mm_cvt_ss2si(_mm_round_ss(_mm_setzero_ps(), _mm_set_ss(v), _MM_FROUND_TO_NEAREST_INT | _MM_FROUND_NO_EXC)); // отбрасываем младшие (bitsFull - bitsCoarse) разрядов: if (cell < 0) cell = 0; int64_t maxv = ((int64_t)1 << bitsFull) - 1; if (cell > maxv) cell = maxv; uint64_t b = (uint64_t)cell >> (bitsFull - bitsCoarse); // инверсия для старшего чанка: uint64_t invMask0 = 0ull; if (bitsCoarse < 63) { uint64_t take = ((uint64_t)1 << bitsCoarse) - 1ull; invMask0 = (invMask[pd] >> (bitsFull - bitsCoarse)) & take; } b ^= invMask0; idx0 |= _pdep_u64(b, pextMask[d]); } if (use_gray) idx0 = gray_decode(idx0); return ((float)idx0 + 0.5f) / (float)scale; }
};
__forceinline void IntervalND::compute_span_level(const MortonND& map) noexcept {
span_level = 0;
for (int d = 0; d < map.dim; ++d) {
uint64_t varying = (i1 ^ i2) & map.pextMask[d];
span_level += _mm_popcnt_u64(varying);
}
// учтём, что ниже старшего чанка есть ещё уровни:
span_level += (map.levels - map.chunk_bits[0]) * map.dim;
span_level = (std::min)(span_level, 11);
}
// ---------------- Стоимость для манипулятора (добавлены приоритет положительных углов и экспоненциальный штраф) ----------------
struct ManipCost final {
int n; bool variableLen; float targetX, targetY; float minTheta;
// коэффициенты для «арки» и экспоненциального штрафа
float archBiasW, archBiasK; // слабый приоритет положительных углов
float sharpW; // вес штрафа за |θ| < minTheta
ManipCost(int _n, bool _variableLen, float _targetX, float _targetY, float _minTheta)
: n(_n), variableLen(_variableLen), targetX(_targetX), targetY(_targetY), minTheta(_minTheta),
archBiasW(0.02f), archBiasK(3.0f), sharpW(0.05f) {
}
__forceinline float operator()(const float* __restrict q, float& out_x, float& out_y) const noexcept {
const float* th = q;
const float* L = variableLen ? (q + n) : nullptr;
__declspec(align(64)) float phi[32], s_arr[32], c_arr[32];
float x = 0.0f, y = 0.0f, phi_acc = 0.0f, penC = 0.0f, archPen = 0.0f;
text// накопление фаз, затем sincos for (int i = 0; i < n; ++i) { phi_acc += th[i]; phi[i] = phi_acc; } FABE13_SINCOS(phi, s_arr, c_arr, n); const float Lc = 1.0f; if (variableLen) { for (int i = 0; i < n; ++i) { float Li = L[i]; x = fmaf(Li, c_arr[i], x); y = fmaf(Li, s_arr[i], y); } } else { for (int i = 0; i < n; ++i) { x = fmaf(Lc, c_arr[i], x); y = fmaf(Lc, s_arr[i], y); } } // штрафы for (int i = 0; i < n; ++i) { float ai = fabsf(th[i]); float v = minTheta - ai; // «слишком острый» (малый) угол if (v > 0.0f) { // экспоненциально нарастающий штраф (через exp2 для скорости) float scale = 2.0f / (minTheta + 1e-6f); penC += sharpW * (exp2f(scale * v) - 1.0f); } // лёгкий приоритет положительных углов (арка): softplus(k * (-theta)) // = ln(1 + exp(k*(-θ))) — минимальный на θ>=0, растёт на θ<0 float t = -th[i] * archBiasK; // быстрая и стабильная softplus: float sp = (t > 10.f) ? t : log1pf(expf(t)); archPen += archBiasW * sp; } float dx = x - targetX, dy = y - targetY; float dist = sqrtf(fmaf(dx, dx, dy * dy)); out_x = x; out_y = y; return dist + penC + archPen; }
};
// ---------------- Попадание по карте (как было) ----------------
__forceinline void HitTest2D_analytic(float x_param, float& out_x1, float& out_x2) {
const float a = gActiveMap.a, inv_lenx = gActiveMap.inv_lenx;
const uint32_t scale = gActiveMap.scale, scale_minus_1 = scale - 1u;
const float lenx = gActiveMap.lenx, leny = gActiveMap.leny, c = gActiveMap.c;
const uint8_t start = gActiveMap.start; const int levels = gActiveMap.levels;
textfloat norm = (x_param - a) * inv_lenx; norm = fminf(fmaxf(norm, 0.0f), 0x1.fffffep-1f); uint32_t idx = (uint32_t)(norm * (float)scale); idx = idx > scale_minus_1 ? scale_minus_1 : idx; float sx = lenx, sy = leny; float x1 = a, x2 = c; uint8_t type = start; int l = levels - 1; while (l >= 0) { const uint32_t q = (idx >> (l * 2)) & 3u; const Step s = g_step_tbl[type][q]; type = s.next; sx *= 0.5f; sy *= 0.5f; x1 += s.dx ? sx : 0.0f; x2 += s.dy ? sy : 0.0f; --l; } out_x1 = x1 + sx * 0.5f; out_x2 = x2 + sy * 0.5f;
}
__forceinline float FindX2D_analytic(float px, float py) {
const float a = gActiveMap.a, b = gActiveMap.b, c = gActiveMap.c, d = gActiveMap.d;
const float lenx = gActiveMap.lenx, leny = gActiveMap.leny; const uint32_t scale = gActiveMap.scale;
const uint8_t start = gActiveMap.start; const int levels = gActiveMap.levels;
const float clamped_px = fminf(fmaxf(px, a), b), clamped_py = fminf(fmaxf(py, c), d);
float sx = lenx, sy = leny; float x0 = a, y0 = c; uint32_t idx = 0u; uint8_t type = start; int l = 0;
while (l < levels) {
sx *= 0.5f; sy *= 0.5f; const float mx = x0 + sx, my = y0 + sy;
const uint32_t tr = (uint32_t)((clamped_px > mx) & (clamped_py > my));
const uint32_t tl = (uint32_t)((clamped_px < mx) & (clamped_py > my));
const uint32_t dl = (uint32_t)((clamped_px < mx) & (clamped_py < my));
const uint32_t none = (uint32_t)(1u ^ (tr | tl | dl));
const uint32_t dd = (tr << 1) | tr | tl | (none << 1);
const InvStep inv = g_inv_tbl[type][dd];
type = inv.next; idx = (idx << 2) | inv.q;
const uint32_t dx = dd >> 1, dy = dd & 1u; x0 += dx ? sx : 0.0f; y0 += dy ? sy : 0.0f; ++l;
}
const float scale_recip = 1.0f / (float)scale;
return fmaf((float)idx * scale_recip, lenx, a);
}
// ---------------- Сообщения ND (как было) ----------------
__declspec(align(16)) struct MultiCrossMsg final {
float intervals[15]; uint8_t count;
template<typename Archive> __declspec(noalias) __forceinline void serialize(Archive& ar, const unsigned int) noexcept { ar& intervals& count; }
};
__declspec(align(16)) struct BestSolutionMsg final {
float bestF, bestX, bestY, bestQ[32]; uint8_t dim;
template<typename Archive> __declspec(noalias) __forceinline void serialize(Archive& ar, const unsigned int) noexcept { ar& bestF& bestX& bestY& bestQ& dim; }
};
__declspec(align(16)) struct CtrlMsgND final {
uint8_t kind; CrossMsg xchg; MultiCrossMsg multiXchg; BestSolutionMsg bestSol;
template<typename Archive> __declspec(noalias) __forceinline void serialize(Archive& ar, const unsigned int) noexcept {
ar& kind; if (kind == 1) ar& xchg; else if (kind == 2) ar& multiXchg; else if (kind == 3) ar& bestSol;
}
};
// ---------------- Генерация сидов (как было) ----------------
static __forceinline int generate_lhs_seeds_lite(const MortonND& map, const int dim, float* __restrict S, int stride, uint32_t seed) {
int temp_dim = dim; const int ns = --temp_dim * temp_dim; uint32_t st = seed;
alignas(32) int permutations[32][256];
for (int d = 0; d < dim; ++d) {
for (int s = 0; s < ns; ++s) permutations[d][s] = s;
for (int s = ns - 1; s > 0; --s) { st ^= st << 13; st ^= st >> 17; st ^= st << 5; int j = st % (s + 1); std::swap(permutations[d][s], permutations[d][j]); }
}
for (int s = 0; s < ns; ++s) {
for (int d = 0; d < dim; ++d) {
st ^= st << 13; st ^= st >> 17; st ^= st << 5;
float u = (st & 0xFFFFFF) * 5.9604645e-8f;
int stratum = permutations[d][s];
float pos = ((float)stratum + u) / (float)ns;
int pd = map.perm[d];
float lo = map.low[pd], hi = map.high[pd];
S[s * stride + d] = fmaf(pos, (hi - lo), lo);
}
}
return ns;
}
static __forceinline int generate_heuristic_seeds(const ManipCost& cost, const MortonND& map, int dim, float* __restrict S, int stride, uint32_t seed) {
const int n = cost.n; const bool VL = cost.variableLen;
const float tx = cost.targetX, ty = cost.targetY;
int total_seeds = 0;
// 1) направленное в сторону цели
{
float* s0 = S + total_seeds * stride;
float phi = atan2f(ty, tx);
float rho = sqrtf(fmaf(tx, tx, ty * ty));
float len = fminf(fmaxf(rho / (float)n, 0.5f), 2.0f); // новый минимум 0.5
for (int i = 0; i < n; ++i) s0[i] = phi / (float)n;
if (VL) for (int i = 0; i < n; ++i) s0[n + i] = len;
total_seeds++;
}
// 2) попеременно
{
float* s1 = S + total_seeds * stride;
float phi = atan2f(ty, tx);
for (int i = 0; i < n; ++i) s1[i] = 0.5f * phi * ((i & 1) ? -1.0f : 1.0f);
if (VL) for (int i = 0; i < n; ++i) s1[n + i] = 1.0f * (0.8f + 0.4f * (float)i / (float)n);
total_seeds++;
}
// 3) затухающая линейка
{
float* s2 = S + total_seeds * stride;
const float inv = (n > 1) ? 1.0f / (float)(n - 1) : 0.0f;
float phi = atan2f(ty, tx);
for (int i = 0; i < n; ++i) { float pr = (float)i * inv; s2[i] = phi * (1.0f - 0.3f * pr); }
if (VL) { for (int i = 0; i < n; ++i) { float si = sinf(1.5f * (float)i); s2[n + i] = (1.0f + 0.2f * si); } }
total_seeds++;
}
// 4) LHS лайт
int lhs_count = generate_lhs_seeds_lite(map, dim, S + total_seeds * stride, stride, seed);
total_seeds += lhs_count;
return total_seeds;
}
// ---------------- Основной ND-ветвитель с улучшениями ----------------
static __forceinline void agp_run_branch_mpi(
const MortonND& map, const ManipCost& cost, int maxIter, float r, bool adaptive, float eps, unsigned seed,
std::vector<IntervalND*>& H, std::vector<float>& bestQ, float& bestF, float& bestX, float& bestY, float M_prior = 1e-3f)
{
const int n = cost.n;
const int dim = n + (cost.variableLen ? n : 0);
textalignas(32) float M_by_span[12]; for (int i = 0; i < 12; ++i) M_by_span[i] = M_prior; float Mmax = M_prior; alignas(64) float q_local[32], phi[32], s_arr[32], c_arr[32], sum_s[32], sum_c[32], q_try[32]; bestQ.reserve(dim); float x = 0.0f, y = 0.0f; int no_improve = 0; auto t_to_idx = [&](float t) -> uint64_t { float tt = (t <= 0.0f) ? 0.0f : (t >= 1.0f ? std::nextafterf(1.0f, 0.0f) : t); uint64_t idx = (uint64_t)((double)tt * (double)map.scale); // индекс старшего чанка if (idx >= map.scale) idx = map.scale - 1ull; return idx; }; auto update_pockets_and_Mmax = [&](IntervalND* I) { const int k = I->span_level; if (I->M > M_by_span[k]) M_by_span[k] = I->M; if (M_by_span[k] > Mmax) Mmax = M_by_span[k]; }; float a = 0.0f, b = 1.0f; // --- локальная оценка + небольшой градиентный шаг + ДОТОЧКА ПО ПОСЛЕДНЕМУ ЗВЕНУ --- auto evalAt = [&](float t) -> float { map.map01ToPoint(t, q_local); float f = cost(q_local, x, y); if (f < bestF * 1.25f) { // подготовим суммы для градиента float acc = 0.0f; for (int i = 0; i < n; ++i) { acc += q_local[i]; phi[i] = acc; } FABE13_SINCOS(phi, s_arr, c_arr, n); float as = 0.0f, ac = 0.0f; for (int k = n - 1; k >= 0; --k) { const float Lk = cost.variableLen ? q_local[n + k] : 1.0f; as += Lk * s_arr[k]; ac += Lk * c_arr[k]; sum_s[k] = as; sum_c[k] = ac; } const float dx = x - cost.targetX, dy = y - cost.targetY; float dist = sqrtf(fmaf(dx, dx, dy * dy)) + 1e-8f; // градиент с учётом новых штрафов: float eta = 0.125f; for (int stepI = 0; stepI < 3; ++stepI) { for (int i = 0; i < n; ++i) { float gpen = 0.0f; // экспоненциальный штраф за малые |θ| { float ai = fabsf(q_local[i]); float v = cost.minTheta - ai; if (v > 0.0f) { float scale = 2.0f / (cost.minTheta + 1e-6f); float e = exp2f(scale * v); float dpen_dtheta = cost.sharpW * (e * 0.69314718055994530941723212145818f * scale) * (-copysignf(1.0f, q_local[i])); gpen += dpen_dtheta; } } // «арка»: производная softplus(k*(-θ)) = -k * sigmoid(k*(-θ)) { float tsg = -q_local[i] * cost.archBiasK; float sig = 1.0f / (1.0f + expf(-tsg)); gpen += -(cost.archBiasW * cost.archBiasK) * sig; } float g = (dx * (-sum_s[i]) + dy * (sum_c[i])) / dist + gpen; q_try[i] = q_local[i] - eta * g; // ГРАНИЦЫ УГЛОВ (рад): для i==0 [-60°, +150°], для остальных [-150°, +150°] const float deg2rad = 3.14159265358979323846f / 180.0f; const float lo0 = -60.0f * deg2rad, hi0 = 150.0f * deg2rad; const float lo = -150.0f * deg2rad, hi = 150.0f * deg2rad; const float Lb = (i == 0) ? lo0 : lo; const float Hb = (i == 0) ? hi0 : hi; if (q_try[i] < Lb) q_try[i] = Lb; else if (q_try[i] > Hb) q_try[i] = Hb; } if (cost.variableLen) { for (int i = 0; i < n; ++i) { float g = (dx * c_arr[i] + dy * s_arr[i]) / dist; float v = q_local[n + i] - eta * g; // НОВЫЕ ГРАНИЦЫ ДЛИН: [0.5; 2.0] if (v < 0.5f) v = 0.5f; else if (v > 2.0f) v = 2.0f; q_try[n + i] = v; } } float x2, y2; float f2 = cost(q_try, x2, y2); if (f2 < f) { memcpy(q_local, q_try, dim * sizeof(float)); f = f2; x = x2; y = y2; break; } eta *= 0.5f; } // ДОТОЧКА ПО ПОСЛЕДНЕМУ УГЛУ — чтобы «довернуть» последний сегмент const int last = n - 1; const float deg2rad = 3.14159265358979323846f / 180.0f; const float lo = (last == 0) ? (-60.0f * deg2rad) : (-150.0f * deg2rad); const float hi = 150.0f * deg2rad; float bestLocF = f; float saved = q_local[last]; for (float delta = 0.05f; delta >= 0.00625f; delta *= 0.5f) { for (int sgn = -1; sgn <= 1; sgn += 2) { float cand = saved + sgn * delta; if (cand < lo) cand = lo; else if (cand > hi) cand = hi; float backup = q_local[last]; q_local[last] = cand; float x2, y2; float f2 = cost(q_local, x2, y2); if (f2 < bestLocF) { bestLocF = f2; x = x2; y = y2; saved = cand; } q_local[last] = backup; } } if (bestLocF < f) { q_local[last] = saved; f = bestLocF; } } if (f < bestF) { bestF = f; bestQ.assign(q_local, q_local + dim); bestX = x; bestY = y; no_improve = 0; } else { ++no_improve; } return f; }; float f_a = evalAt(a), f_b = evalAt(b); const int K = (std::min)((std::max)(2 * dim, 8), 128); H.reserve((size_t)maxIter + K + 16); const int rank = g_world->rank(); const int world = g_world->size(); alignas(64) float seeds[256 * 32]; const int seedCnt = generate_heuristic_seeds(cost, map, dim, seeds, 32, seed + rank * 7919u); for (int i = 0; i < seedCnt; ++i) { const float* s = seeds + i * 32; float t_seed = map.pointToT(s); float interval_size = (i < 3) ? (0.0004f * (float)dim) : (0.00031f * (float)dim) * exp2f((1.0f / (float)(seedCnt - 4)) * log2f(0.00025f / 0.00031f) * (float)(i - 3)); float t1 = fmaxf(a, t_seed - interval_size), t2 = fminf(b, t_seed + interval_size); if (t2 <= t1) continue; alignas(64) float q1[32], q2[32]; float x1, y1, x2, y2; map.map01ToPoint(t1, q1); float f1 = cost(q1, x1, y1); map.map01ToPoint(t2, q2); float f2 = cost(q2, x2, y2); IntervalND* I = new IntervalND(t1, t2, f1, f2); I->i1 = t_to_idx(t1); I->i2 = t_to_idx(t2); I->diam = map.block_diameter(I->i1, I->i2); I->compute_span_level(map); I->set_metric(I->diam); update_pockets_and_Mmax(I); I->ChangeCharacteristic(r * Mmax); if (i < 3) I->R *= fmaf(0.01f, (float)dim, 0.85f); else { float start_mult = 0.214f * (float)dim; float end_mult = 0.174f * (float)dim; float mult = start_mult * exp2f((1.0f / (float)(seedCnt - 4)) * log2f(end_mult / start_mult) * (float)(i - 3)); I->R *= mult; } H.emplace_back(I); std::push_heap(H.begin(), H.end(), ComparePtrND); if (f1 < bestF) { bestF = f1; bestQ.assign(q1, q1 + dim); bestX = x1; bestY = y1; } if (f2 < bestF) { bestF = f2; bestQ.assign(q2, q2 + dim); bestX = x2; bestY = y2; } } float prev_t = a, prev_f = f_a; for (int k = 1; k <= K; ++k) { float t = a + (b - a) * ((float)k / (K + 1)) + (float)rank / (float)(world * (K + 1)); float f = evalAt(t); IntervalND* I = new IntervalND(prev_t, t, prev_f, f); I->i1 = t_to_idx(prev_t); I->i2 = t_to_idx(t); I->diam = map.block_diameter(I->i1, I->i2); I->compute_span_level(map); I->set_metric(I->diam); update_pockets_and_Mmax(I); I->ChangeCharacteristic(r * Mmax); H.emplace_back(I); std::push_heap(H.begin(), H.end(), ComparePtrND); prev_t = t; prev_f = f; } IntervalND* tail = new IntervalND(prev_t, b, prev_f, f_b); tail->i1 = t_to_idx(prev_t); tail->i2 = t_to_idx(b); tail->diam = map.block_diameter(tail->i1, tail->i2); tail->compute_span_level(map); tail->set_metric(tail->diam); update_pockets_and_Mmax(tail); tail->ChangeCharacteristic(r * Mmax); H.emplace_back(tail); std::push_heap(H.begin(), H.end(), ComparePtrND); float dmax = b - a, initial_len = dmax, thr03 = 0.3f * initial_len, inv_thr03 = 1.0f / thr03; int it = 0; // новые «мягко-зависящие-от-размерности» константы auto kickEveryByDim = [&](int dim) -> int { // 60..120, убывает экспоненциально с размерностью float z = 120.0f * exp2f(-0.05f * (float)dim); if (z < 60.0f) z = 60.0f; return (int)z; }; auto noImproveThrByDim = [&](int dim) -> int { float z = 80.0f * exp2f(-0.08f * (float)dim); if (z < 30.0f) z = 30.0f; return (int)z; }; while (it < maxIter) { // подталкивание около лучшего — зависящее от размерности if ((it % kickEveryByDim(dim)) == 0 && no_improve > noImproveThrByDim(dim)) { float t_best = map.pointToT(bestQ.data()); for (int i = 0; i < 2; ++i) { float off = (i == 0) ? 0.01f : -0.01f; float t_seed = fminf(b, fmaxf(a, t_best + off)); float f_seed = evalAt(t_seed); IntervalND* J = new IntervalND(t_seed - 0.005f, t_seed + 0.005f, f_seed, f_seed); J->i1 = t_to_idx(t_seed - 0.005f); J->i2 = t_to_idx(t_seed + 0.005f); J->diam = map.block_diameter(J->i1, J->i2); J->compute_span_level(map); J->set_metric(J->diam); update_pockets_and_Mmax(J); J->ChangeCharacteristic(r * Mmax); J->R *= 0.9f; H.emplace_back(J); std::push_heap(H.begin(), H.end(), ComparePtrND); } no_improve = 0; } const float p = fmaf(-1.0f / initial_len, dmax, 1.0f); bool stagnation = (no_improve > 100) && (it > 270); // Плавная адаптация «T» к размерности (вместо констант 264/277) float A = 200.0f + 64.0f * exp2f(-0.06f * (float)dim); float B = 210.0f + 67.0f * exp2f(-0.06f * (float)dim); const int T = (int)fmaf(-expm1f(p), A, B); float r_eff = fmaxf(1.0f, r * (0.7f + 0.3f * (1.0f - p))); std::pop_heap(H.begin(), H.end(), ComparePtrND); IntervalND* cur = H.back(); H.pop_back(); const float x1 = cur->x1, x2 = cur->x2, y1 = cur->y1, y2 = cur->y2; float m = r_eff * Mmax; float tNew = Shag(m, x1, x2, y1, y2, (float)dim, r); tNew = fminf(fmaxf(tNew, a), b); float fNew = evalAt(tNew); IntervalND* L = new IntervalND(x1, tNew, y1, fNew); IntervalND* Rv = new IntervalND(tNew, x2, fNew, y2); L->i1 = t_to_idx(x1); L->i2 = t_to_idx(tNew); Rv->i1 = t_to_idx(tNew); Rv->i2 = t_to_idx(x2); L->diam = map.block_diameter(L->i1, L->i2); Rv->diam = map.block_diameter(Rv->i1, Rv->i2); L->compute_span_level(map); Rv->compute_span_level(map); L->set_metric(L->diam); Rv->set_metric(Rv->diam); float Mloc = (std::max)(L->M, Rv->M); update_pockets_and_Mmax(L); update_pockets_and_Mmax(Rv); const float prevMmax = Mmax; if (Mloc > Mmax) Mmax = Mloc; m = r_eff * Mmax; if (adaptive) { float len1 = tNew - x1, len2 = x2 - tNew; if (len1 + len2 == dmax) { dmax = (std::max)(len1, len2); for (auto pI : H) { float Ls = pI->x2 - pI->x1; if (Ls > dmax) dmax = Ls; } } // НЕ ТРОГАЕМ усл. с thr03 — строго как просили if ((thr03 > dmax && !(it % 3)) || (10.0f * dmax < initial_len)) { const float progress = fmaf(-dmax, inv_thr03, 1.0f); const float alpha = progress * progress; const float beta = fmaf(-alpha, 1.0f, 2.0f); const float MULT = (1.0f / dmax) * Mmax; const float global_coeff = fmaf(MULT, r_eff, -MULT); const float GF = fmaf(beta, global_coeff, 0.0f); L->ChangeCharacteristic(fmaf(GF, len1, L->M * alpha)); Rv->ChangeCharacteristic(fmaf(GF, len2, Rv->M * alpha)); size_t sz = H.size(); RecomputeR_AffineM_AVX2_ND(H.data(), sz, GF, alpha); std::make_heap(H.begin(), H.end(), ComparePtrND); } else { if (Mloc > prevMmax) { L->ChangeCharacteristic(m); Rv->ChangeCharacteristic(m); if (Mloc > 1.15f * prevMmax) { size_t sz = H.size(); RecomputeR_ConstM_AVX2_ND(H.data(), sz, m); std::make_heap(H.begin(), H.end(), ComparePtrND); } } else { L->ChangeCharacteristic(m); Rv->ChangeCharacteristic(m); } } } else { if (Mloc > prevMmax) { L->ChangeCharacteristic(m); Rv->ChangeCharacteristic(m); if (Mloc > 1.15f * prevMmax) { size_t sz = H.size(); RecomputeR_ConstM_AVX2_ND(H.data(), sz, m); std::make_heap(H.begin(), H.end(), ComparePtrND); } } else { L->ChangeCharacteristic(m); Rv->ChangeCharacteristic(m); } } H.push_back(L); std::push_heap(H.begin(), H.end(), ComparePtrND); H.push_back(Rv); std::push_heap(H.begin(), H.end(), ComparePtrND); if (!H.empty()) { _mm_prefetch((const char*)H[0], _MM_HINT_T0); if (H.size() > 1) _mm_prefetch((const char*)H[1], _MM_HINT_T0); } IntervalND* top = H.front(); float interval_len = top->x2 - top->x1; // обмен/терминация по расписанию: bool want_term = (exp2f(log2f(interval_len) / (float)dim) < eps) || (it == maxIter - 1); if (!(it % T) || want_term) { CtrlMsgND out; out.kind = want_term ? 0 : 2; if (!want_term) { uint8_t cnt = (uint8_t)((H.size() >= 3) ? 3 : H.size()); out.multiXchg.count = cnt; float* dest = out.multiXchg.intervals; IntervalND* t1 = H[0]; IntervalND* t2 = (H.size() > 1 ? H[1] : H[0]); IntervalND* t3 = (H.size() > 2 ? H[2] : H[H.size() - 1]); IntervalND* tops[3] = { t1, t2, t3 }; for (uint8_t i2 = 0; i2 < cnt; ++i2) { IntervalND* Tt = tops[i2]; dest[0] = Tt->x1; dest[1] = 0.0f; dest[2] = Tt->x2; dest[3] = 0.0f; dest[4] = Tt->R; dest += 5; } } for (int i2 = 0; i2 < world; ++i2) if (i2 != rank) g_world->isend(i2, 0, out); if (want_term) return; } if (!(it % 500)) { CtrlMsgND out; out.kind = 3; out.bestSol.bestF = bestF; out.bestSol.bestX = bestX; out.bestSol.bestY = bestY; out.bestSol.dim = (uint8_t)bestQ.size(); memcpy(out.bestSol.bestQ, bestQ.data(), bestQ.size() * sizeof(float)); for (int i2 = 0; i2 < world; ++i2) if (i2 != rank) g_world->isend(i2, 0, out); } while (g_world->iprobe(boost::mpi::any_source, 0)) { CtrlMsgND in; g_world->recv(boost::mpi::any_source, 0, in); if (in.kind == 0) return; else if (in.kind == 1) { float sx = in.xchg.s_x1, ex = in.xchg.e_x1; if (sx < 0.0f) sx = 0.0f; if (ex > 1.0f) ex = 1.0f; if (ex > sx) { alignas(64) float tmp[32]; float tx, ty; map.map01ToPoint(sx, tmp); float y1i = cost(tmp, tx, ty); map.map01ToPoint(ex, tmp); float y2i = cost(tmp, tx, ty); IntervalND* inj = new IntervalND(sx, ex, y1i, y2i); inj->i1 = t_to_idx(sx); inj->i2 = t_to_idx(ex); inj->diam = map.block_diameter(inj->i1, inj->i2); inj->compute_span_level(map); inj->set_metric(inj->diam); update_pockets_and_Mmax(inj); inj->ChangeCharacteristic(r * Mmax); if (!H.empty()) { _mm_prefetch((const char*)H[0], _MM_HINT_T0); if (H.size() > 1) _mm_prefetch((const char*)H[1], _MM_HINT_T0); } IntervalND* topH = H.front(); if (inj->R > 1.15f * topH->R) { float p2 = fmaf(-1.0f / initial_len, dmax, 1.0f); float k = (no_improve > 100 && it > 270) ? fmaf(0.5819767068693265f, expm1f(p2), 0.3f) : fmaf(0.3491860241215959f, expm1f(p2), 0.6f); inj->R = in.xchg.Rtop * k; H.emplace_back(inj); std::push_heap(H.begin(), H.end(), ComparePtrND); } } } else if (in.kind == 2) { const MultiCrossMsg& mX = in.multiXchg; for (uint8_t ii = 0; ii < mX.count; ++ii) { const float* d = &mX.intervals[ii * 5]; float sx = d[0], ex = d[2]; if (sx < 0.0f) sx = 0.0f; if (ex > 1.0f) ex = 1.0f; if (ex > sx) { alignas(64) float tmp[32]; float tx, ty; map.map01ToPoint(sx, tmp); float y1i = cost(tmp, tx, ty); map.map01ToPoint(ex, tmp); float y2i = cost(tmp, tx, ty); IntervalND* inj = new IntervalND(sx, ex, y1i, y2i); inj->i1 = t_to_idx(sx); inj->i2 = t_to_idx(ex); inj->diam = map.block_diameter(inj->i1, inj->i2); inj->compute_span_level(map); inj->set_metric(inj->diam); update_pockets_and_Mmax(inj); inj->ChangeCharacteristic(r * Mmax); if (!H.empty()) { _mm_prefetch((const char*)H[0], _MM_HINT_T0); if (H.size() > 1) _mm_prefetch((const char*)H[1], _MM_HINT_T0); } IntervalND* topH = H.front(); if (inj->R > 1.15f * topH->R) { float p2 = fmaf(-1.0f / initial_len, dmax, 1.0f); float k = (no_improve > 100 && it > 270) ? fmaf(0.5819767068693265f, expm1f(p2), 0.3f) : fmaf(0.3491860241215959f, expm1f(p2), 0.6f); inj->R = d[4] * k; H.emplace_back(inj); std::push_heap(H.begin(), H.end(), ComparePtrND); } } } } else if (in.kind == 3) { const BestSolutionMsg& bm = in.bestSol; if (bm.bestF < bestF * 1.15f) { alignas(64) float tmp_q[32]; memcpy(tmp_q, bm.bestQ, bm.dim * sizeof(float)); float t_best = map.pointToT(tmp_q); float t1 = fmaxf(a, t_best - 0.001f), t2 = fminf(b, t_best + 0.001f); if (t2 > t1) { alignas(64) float tq1[32], tq2[32]; float xx1, yy1, xx2, yy2; map.map01ToPoint(t1, tq1); float f1 = cost(tq1, xx1, yy1); map.map01ToPoint(t2, tq2); float f2 = cost(tq2, xx2, yy2); IntervalND* I = new IntervalND(t1, t2, f1, f2); I->i1 = t_to_idx(t1); I->i2 = t_to_idx(t2); I->diam = map.block_diameter(I->i1, I->i2); I->compute_span_level(map); I->set_metric(I->diam); update_pockets_and_Mmax(I); I->ChangeCharacteristic(r * Mmax); I->R *= 0.90f; H.emplace_back(I); std::push_heap(H.begin(), H.end(), ComparePtrND); } if (bm.bestF < bestF) { bestF = bm.bestF; bestX = bm.bestX; bestY = bm.bestY; bestQ.assign(bm.bestQ, bm.bestQ + bm.dim); } } } } ++it; }
}
// ---------------- ВНЕШНИЕ ФУНКЦИИ DLL ----------------
// Главная ND оптимизация с новыми границами углов/длин
extern "C" __declspec(dllexport) __declspec(noalias)
void AGP_Manip2D(int nSegments, bool variableLengths, float minTheta, float targetX, float targetY,
int peanoLevels, int maxIterPerBranch, float r, bool adaptiveMode, float epsilon,
unsigned int seed, float** out_bestQ, size_t* out_bestQLen, float* out_bestX,
float* out_bestY, float* out_bestF)
{
Slab* const __restrict slab = tls.local(); slab->current = slab->base;
while (g_world->iprobe(boost::mpi::any_source, 0)) { CtrlMsg dummy; g_world->recv(boost::mpi::any_source, 0, dummy); }
const int dim = nSegments + (variableLengths ? nSegments : 0);
text// Перестановка/инверсии (как было) g_mc.permCache.resize(dim); for (int i = 0; i < dim; ++i) g_mc.permCache[i] = i; uint32_t s = g_mc.baseSeed; for (int i = dim - 1; i > 0; --i) { s ^= s << 13; s ^= s >> 17; s ^= s << 5; uint32_t j = s % (uint32_t)(i + 1); std::swap(g_mc.permCache[i], g_mc.permCache[j]); } g_mc.invMaskCache.resize(dim); for (int k = 0; k < dim; ++k) { s ^= s << 13; s ^= s >> 17; s ^= s << 5; g_mc.invMaskCache[k] = (uint64_t)s; } // НОВЫЕ ГРАНИЦЫ: // углы: θ0 ∈ [-60°, +150°], θi>0 ∈ [-150°, +150°]; длины ∈ [0.5, 2.0] const float deg2rad = 3.14159265358979323846f / 180.0f; const float theta0Min = -60.0f * deg2rad, theta0Max = 150.0f * deg2rad; const float thetaMin = -150.0f * deg2rad, thetaMax = 150.0f * deg2rad; const float lenMin = 0.5f, lenMax = 2.0f; std::vector<float> low; low.reserve(dim); std::vector<float> high; high.reserve(dim); for (int i = 0; i < nSegments; ++i) { low.push_back(i == 0 ? theta0Min : thetaMin); high.push_back(i == 0 ? theta0Max : thetaMax); } if (variableLengths) { for (int i = 0; i < nSegments; ++i) { low.push_back(lenMin); high.push_back(lenMax); } } ManipCost cost(nSegments, variableLengths, targetX, targetY, minTheta); const int rank = g_world->rank(), world = g_world->size(); std::vector<float> bestQ; float bestF = FLT_MAX, bx = 0.0f, by = 0.0f; const int levels0 = (std::min)(peanoLevels, 8); const int maxIter0 = (int)(maxIterPerBranch * 0.2f); MortonND map0(dim, levels0, low.data(), high.data(), g_mc); std::vector<IntervalND*> H_coarse; std::vector<float> bestQ_coarse; float bestF_coarse = FLT_MAX, bx_coarse = 0.0f, by_coarse = 0.0f; float M_prior = (variableLengths ? 2.0f * nSegments : 2.0f * nSegments) * (1.0f / (float)(1u << levels0)); if (variableLengths) M_prior += 1.41421356237f * (1.0f / (float)(1u << levels0)); agp_run_branch_mpi(map0, cost, maxIter0, r, adaptiveMode, epsilon, seed, H_coarse, bestQ_coarse, bestF_coarse, bx_coarse, by_coarse, M_prior); if (bestF_coarse < bestF) { bestF = bestF_coarse; bestQ = bestQ_coarse; bx = bx_coarse; by = by_coarse; } if (levels0 < peanoLevels) { while (g_world->iprobe(boost::mpi::any_source, 0)) { CtrlMsg dummy; g_world->recv(boost::mpi::any_source, 0, dummy); } MortonND map1(dim, peanoLevels, low.data(), high.data(), g_mc); std::vector<IntervalND*> H_fine; std::vector<float> bestQ_fine = bestQ; float bestF_fine = bestF, bx_fine = bx, by_fine = by; float M_prior_fine = (variableLengths ? 2.0f * nSegments : 2.0f * nSegments) * (1.0f / (float)(1u << peanoLevels)); if (variableLengths) M_prior_fine += 1.41421356237f * (1.0f / (float)(1u << peanoLevels)); std::sort(H_coarse.begin(), H_coarse.end(), [](const IntervalND* a, const IntervalND* b) { return a->R < b->R; }); const size_t topCount = (size_t)(H_coarse.size() * 0.3f); auto t_to_idx_fine = [&](float t) -> uint64_t { float tt = (t <= 0.0f) ? 0.0f : (t >= 1.0f ? std::nextafterf(1.0f, 0.0f) : t); uint64_t idx = (uint64_t)((double)tt * (double)map1.scale); if (idx >= map1.scale) idx = map1.scale - 1ull; return idx; }; for (size_t i = 0; i < topCount && i < H_coarse.size(); ++i) { const IntervalND* C = H_coarse[i]; alignas(64) float q1[32], q2[32]; float x1, y1, x2, y2; map1.map01ToPoint(C->x1, q1); float f1 = cost(q1, x1, y1); map1.map01ToPoint(C->x2, q2); float f2 = cost(q2, x2, y2); IntervalND* I = new IntervalND(C->x1, C->x2, f1, f2); I->i1 = t_to_idx_fine(C->x1); I->i2 = t_to_idx_fine(C->x2); I->diam = map1.block_diameter(I->i1, I->i2); I->set_metric(I->diam); H_fine.push_back(I); if (f1 < bestF_fine) { bestF_fine = f1; bestQ_fine.assign(q1, q1 + dim); bx_fine = x1; by_fine = y1; } if (f2 < bestF_fine) { bestF_fine = f2; bestQ_fine.assign(q2, q2 + dim); bx_fine = x2; by_fine = y2; } } std::make_heap(H_fine.begin(), H_fine.end(), ComparePtrND); agp_run_branch_mpi(map1, cost, maxIterPerBranch - maxIter0, r, adaptiveMode, epsilon, seed, H_fine, bestQ_fine, bestF_fine, bx_fine, by_fine, M_prior_fine); if (bestF_fine < bestF) { bestF = bestF_fine; bestQ = bestQ_fine; bx = bx_fine; by = by_fine; } } BestSolutionMsg bp; bp.bestF = bestF; bp.bestX = bx; bp.bestY = by; bp.dim = (uint8_t)bestQ.size(); // Копируем данные через memcpy как в BestSolutionMsg memcpy(bp.bestQ, bestQ.data(), bp.dim * sizeof(float)); if (rank) { g_world->send(0, 2, bp); } else { int received_count = 0; BestSolutionMsg in; BestSolutionMsg best = bp; while (received_count < world - 1) { g_world->recv(boost::mpi::any_source, 2, in); if (in.bestF < best.bestF) best = in; ++received_count; } *out_bestQLen = (size_t)best.dim; *out_bestQ = (float*)CoTaskMemAlloc(sizeof(float) * (*out_bestQLen)); // Используем memcpy как при работе с BestSolutionMsg memcpy(*out_bestQ, best.bestQ, sizeof(float) * (*out_bestQLen)); *out_bestX = best.bestX; *out_bestY = best.bestY; *out_bestF = best.bestF; }
}
extern "C" __declspec(dllexport) __declspec(noalias) __forceinline int AgpInit(int peanoLevel, float a, float b, float c, float d) {
g_env = new boost::mpi::environment();
g_world = new boost::mpi::communicator();
_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
_MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);
const int rank = g_world->rank();
const int world_size = g_world->size();
if (world_size == 4) { new (&gActiveMap) Peano2DMap(peanoLevel, a, b, c, d, rank & 3); }
g_mc.baseSeed = fmaf(0x9E3779B9u, rank, 0x9E3779B9u);
return rank;
}
// 1D/2D тестовые функции/обёртки — без функциональных изменений (оставлены как были)
__forceinline float ShekelFunc(float x, float seed) {
int i = 0; float st = seed, r1, r2, res = 0.0f;
while (i < 10) {
XOR_RAND(st, r1); float xp = fmaf(-r1, 10.0f, x); XOR_RAND(st, r1); XOR_RAND(st, r2);
float d = fmaf(fmaf(r1, 20.0f, 5.0f), xp * xp, fmaf(r2, 0.2f, 1.0f)); d = copysignf(fmaxf(fabsf(d), FLT_MIN), d); res -= 1.0f / d; ++i;
}
return res;
}
__forceinline float RastriginFunc(float x1, float x2) {
const float t = fmaf(x1, x1, x2 * x2); float c1, c2;
FABE13_COS(6.28318530717958647692f * x1, c1); FABE13_COS(6.28318530717958647692f * x2, c2);
return (t - fmaf(c1 + c2, 10.0f, -14.6f)) * fmaf(-t, 0.25f, 18.42f);
}
__forceinline float HillFunc(float x, float seed) {
int j = 0; __declspec(align(32)) float ang[14u];
float st = 6.28318530717958647692f * x; while (j < 14) { ang[j] = st * (float)(j + 1); ++j; }
__declspec(align(32)) float sv[14u], cv[14u]; FABE13_SINCOS(ang, sv, cv, 14u);
float state = seed, r1, r2; XOR_RAND(state, r1); float res = fmaf(r1, 2.0f, -1.1f); --j;
while (j >= 0) { XOR_RAND(state, r1); XOR_RAND(state, r2); res += fmaf(fmaf(r1, 2.0f, -1.1f), sv[j], fmaf(r2, 2.0f, -1.1f) * cv[j]); --j; }
return res;
}
__forceinline float GrishaginFunc(float x1, float x2, float seed) {
int j = 0; __declspec(align(32)) float aj[8u], ak[8u];
while (j < 8) { float pj = 3.14159265358979323846f * (float)(j + 1); aj[j] = pj * x1; ak[j] = pj * x2; ++j; }
__declspec(align(32)) float sj[8u], cj[8u], sk[8u], ck[8u]; FABE13_SINCOS(aj, sj, cj, 8u); FABE13_SINCOS(ak, sk, ck, 8u);
--j; float p1 = 0.0f, p2 = 0.0f; float st = seed, r1, r2;
while (j >= 0) {
size_t k = 0u; while (k < 8u) {
float s = sj[j] * sj[j]; float c = ck[k] * ck[k];
XOR_RAND_GRSH(st, r1); XOR_RAND_GRSH(st, r2); p1 = fmaf(r1, s, fmaf(r2, c, p1));
XOR_RAND_GRSH(st, r1); XOR_RAND_GRSH(st, r2); p2 = fmaf(-r1, c, fmaf(r2, s, p2)); ++k;
}
--j;
}
return -sqrtf(fmaf(p1, p1, p2 * p2));
}
// Тестовые 1D/2D API оставлены без поведенческих изменений
extern "C" __declspec(dllexport) __declspec(noalias)
void AGP_1D(float global_iterations, float a, float b, float r, bool mode, float epsilon, float seed,
float** out_data, size_t* out_len) {
Slab* const __restrict slab = tls.local(); slab->current = slab->base; int schetchick = 0;
const float initial_length = b - a; float dmax = initial_length;
const float threshold_03 = 0.3f * initial_length, inv_threshold_03 = 1.0f / threshold_03;
const float start_val = ShekelFunc(a, seed); float best_f = ShekelFunc(b, seed);
float x_Rmax_1 = a, x_Rmax_2 = b; float y_Rmax_1 = start_val, y_Rmax_2 = best_f;
std::vector<float, boost::alignment::aligned_allocator<float, 16u>> Extr;
std::vector<Interval1D*, boost::alignment::aligned_allocator<Interval1D*, 64u>> R;
Extr.reserve((size_t)global_iterations << 2u); R.reserve((size_t)global_iterations << 1u);
R.emplace_back(new Interval1D(a, b, start_val, best_f, 1.0f)); float Mmax = R.front()->M; float m = r * Mmax;
textwhile (true) { const float new_point = Shag(m, x_Rmax_1, x_Rmax_2, y_Rmax_1, y_Rmax_2, 1.0f, r); const float new_value = ShekelFunc(new_point, seed); if (new_value < best_f) { best_f = new_value; Extr.emplace_back(best_f); Extr.emplace_back(new_point); } std::pop_heap(R.begin(), R.end(), ComparePtr1D); const Interval1D* pro = R.back(); const float new_x1 = pro->x1, new_x2 = pro->x2; const float len2 = new_x2 - new_point, len1 = new_point - new_x1; const float interval_len = (len1 < len2 ? len1 : len2); if (++schetchick == (int)global_iterations || interval_len < epsilon) { Extr.emplace_back((float)schetchick); Extr.emplace_back(interval_len); *out_len = Extr.size(); *out_data = (float*)CoTaskMemAlloc(sizeof(float) * (*out_len)); memcpy(*out_data, Extr.data(), sizeof(float) * (*out_len)); return; } Interval1D* curr = new Interval1D(new_x1, new_point, pro->y1, new_value, 1.0f); Interval1D* curr1 = new Interval1D(new_point, new_x2, new_value, pro->y2, 1.0f); const float currM = curr->M > curr1->M ? curr->M : curr1->M; const size_t r_size = R.size(); if (mode) { if (len2 + len1 == dmax) { dmax = len2 > len1 ? len2 : len1; for (auto p : R) { float L = p->x2 - p->x1; if (L > dmax) dmax = L; } } if (threshold_03 > dmax && !(schetchick % 3) || 10.0f * dmax < initial_length) { if (currM > Mmax) { Mmax = currM; m = r * Mmax; } const float progress = fmaf(-inv_threshold_03, dmax, 1.0f); const float alpha = progress * progress; const float betta = 2.0f - alpha; const float MULT = (1.0f / dmax) * Mmax; const float global_coeff = fmaf(MULT, r, -MULT); const float GF = betta * global_coeff; curr->ChangeCharacteristic(fmaf(GF, len1, curr->M * alpha)); curr1->ChangeCharacteristic(fmaf(GF, len2, curr1->M * alpha)); RecomputeR_AffineM_AVX2_1D(R.data(), r_size, GF, alpha); std::make_heap(R.begin(), R.end(), ComparePtr1D); } else { if (currM > Mmax) { if (currM < 1.15f * Mmax) { Mmax = currM; m = r * Mmax; curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); } else { Mmax = currM; m = r * Mmax; curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); RecomputeR_ConstM_AVX2_1D(R.data(), r_size, m); std::make_heap(R.begin(), R.end(), ComparePtr1D); } } else { curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); } } } else { if (currM > Mmax) { if (currM < 1.15f * Mmax) { Mmax = currM; m = r * Mmax; curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); } else { Mmax = currM; m = r * Mmax; curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); RecomputeR_ConstM_AVX2_1D(R.data(), r_size, m); std::make_heap(R.begin(), R.end(), ComparePtr1D); } } else { curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); } } R.back() = curr; std::push_heap(R.begin(), R.end(), ComparePtr1D); R.emplace_back(curr1); std::push_heap(R.begin(), R.end(), ComparePtr1D); const Interval1D* top = R.front(); x_Rmax_1 = top->x1; x_Rmax_2 = top->x2; y_Rmax_1 = top->y1; y_Rmax_2 = top->y2; }
}
extern "C" __declspec(dllexport) __declspec(noalias)
void AGP_2D(const float N, const float global_iterations, const float a, const float b, const float c,
const float d, const float r, const bool mode, const float epsilon, const float seed,
float** const __restrict out_data, size_t* const __restrict out_len) noexcept
{
Slab* const __restrict slab = tls.local(); slab->current = slab->base;
int schetchick = 0, no_improve = 0; const int rank = g_world->rank(); const int world_size = g_world->size();
while (g_world->iprobe(boost::mpi::any_source, 0)) { CtrlMsg dummy; g_world->recv(boost::mpi::any_source, 0, dummy); }
const float inv_divider = ldexpf(1.0f, -((gActiveMap.levels << 1) + 1));
const float x_addition = (b - a) * inv_divider, y_addition = (d - c) * inv_divider;
const float true_start = a + x_addition, true_end = b - x_addition;
float x_Rmax_1 = true_start, x_Rmax_2 = true_end;
const float initial_length = x_Rmax_2 - x_Rmax_1; float dmax = initial_length;
const float threshold_03 = 0.3f * initial_length, inv_threshold_03 = 1.0f / threshold_03;
const float start_val = rank % 3 ? RastriginFunc(true_end, d - y_addition) : RastriginFunc(true_start, c + y_addition);
float best_f = rank % 2 ? RastriginFunc(true_start, d - y_addition) : RastriginFunc(true_end, c + y_addition);
float y_Rmax_1 = start_val, y_Rmax_2 = best_f;
std::vector<float, boost::alignment::aligned_allocator<float, 16u>> Extr;
std::vector<Interval1D* __restrict, boost::alignment::aligned_allocator<Interval1D* __restrict, 64u>> R;
Extr.clear(); Extr.reserve(static_cast<size_t>(global_iterations) << 2u);
R.clear(); R.reserve(static_cast<size_t>(global_iterations) << 1u);
R.emplace_back(new Interval1D(true_start, true_end, start_val, best_f, 2.0f));
const Interval1D* __restrict top_ptr;
float Mmax = R.front()->M, m = r * Mmax;
while (true) {
const float interval_len = x_Rmax_2 - x_Rmax_1;
const bool stagnation = no_improve > 100 && schetchick > 270;
const float p = fmaf(-1.0f / initial_length, dmax, 1.0f);
while (g_world->iprobe(boost::mpi::any_source, 0)) {
CtrlMsg in; g_world->recv(boost::mpi::any_source, 0, in);
if (in.kind) {
if (!rank) {
Extr.emplace_back((float)schetchick); Extr.emplace_back(interval_len); *out_len = Extr.size();
out_data = reinterpret_cast<float __restrict>(CoTaskMemAlloc(sizeof(float) * (*out_len)));
memcpy(*out_data, Extr.data(), sizeof(float) * (out_len));
}
return;
}
const float sx = FindX2D_analytic(in.xchg.s_x1, in.xchg.s_x2);
const float ex = FindX2D_analytic(in.xchg.e_x1, in.xchg.e_x2);
Interval1D const __restrict injected = new Interval1D(sx, ex,
RastriginFunc(in.xchg.s_x1, in.xchg.s_x2), RastriginFunc(in.xchg.e_x1, in.xchg.e_x2), 2.0f);
injected->ChangeCharacteristic(m);
if (injected->R > 1.15f * top_ptr->R) {
const float k = stagnation ? fmaf(0.5819767068693265f, expm1f(p), 0.3f) : fmaf(0.3491860241215959f, expm1f(p), 0.6f);
injected->R = in.xchg.Rtop * k;
R.emplace_back(injected); std::push_heap(R.begin(), R.end(), ComparePtr1D);
}
}
// исходная формула T оставлена (это тестовая 2D), менять здесь не обязательно
const int T = (int)fmaf(-expm1f(p), 264.0f, 277.0f);
const bool want_term = interval_len < epsilon || schetchick == static_cast<int>(global_iterations);
if (!(++schetchick % T) || stagnation || want_term) {
CtrlMsg out; out.kind = want_term;
if (!out.kind) {
float s_x1, s_x2, e_x1, e_x2; HitTest2D_analytic(top_ptr->x1, s_x1, s_x2); HitTest2D_analytic(top_ptr->x2, e_x1, e_x2);
out.xchg = CrossMsg{ s_x1, s_x2, e_x1, e_x2, top_ptr->R };
}
for (int i = 0; i < world_size; ++i) if (i != rank) g_world->isend(i, 0, out);
if (out.kind) {
if (!rank) {
Extr.emplace_back((float)schetchick); Extr.emplace_back(interval_len); *out_len = Extr.size();
out_data = reinterpret_cast<float __restrict>(CoTaskMemAlloc(sizeof(float) * (out_len)));
memcpy(out_data, Extr.data(), sizeof(float) * (out_len));
}
return;
}
}
const float new_point = Shag(m, x_Rmax_1, x_Rmax_2, y_Rmax_1, y_Rmax_2, 2.0f, r);
float new_x1_val, new_x2_val; HitTest2D_analytic(new_point, new_x1_val, new_x2_val);
const float new_value = RastriginFunc(new_x1_val, new_x2_val);
if (new_value < best_f) { best_f = new_value; Extr.emplace_back(best_f); Extr.emplace_back(new_x1_val); Extr.emplace_back(new_x2_val); no_improve = 0; }
else { ++no_improve; }
std::pop_heap(R.begin(), R.end(), ComparePtr1D);
Interval1D const __restrict promej = R.back();
const float segment_x1 = promej->x1, segment_x2 = promej->x2;
const float len2 = segment_x2 - new_point, len1 = new_point - segment_x1;
Interval1D const __restrict curr = new Interval1D(segment_x1, new_point, promej->y1, new_value, 2.0f);
Interval1D const __restrict curr1 = new Interval1D(new_point, segment_x2, new_value, promej->y2, 2.0f);
const float currM = (std::max)(curr->M, curr1->M); const size_t r_size = R.size();
if (mode) {
if (len2 + len1 == dmax) { dmax = (std::max)(len1, len2); for (auto pI : R) { float L = pI->x2 - pI->x1; if (L > dmax) dmax = L; } }
if (threshold_03 > dmax && !(schetchick % 3) || 10.0f * dmax < initial_length) {
if (currM > Mmax) { Mmax = currM; m = r * Mmax; }
const float progress = fmaf(-inv_threshold_03, dmax, 1.0f);
const float alpha = progress * progress; const float betta = 2.0f - alpha;
const float MULTIPLIER = (1.0f / dmax) * Mmax;
const float global_coeff = fmaf(MULTIPLIER, r, -MULTIPLIER);
const float GLOBAL_FACTOR = betta * global_coeff;
curr->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, len1, curr->M * alpha));
curr1->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, len2, curr1->M * alpha));
RecomputeR_AffineM_AVX2_1D(R.data(), r_size, GLOBAL_FACTOR, alpha);
std::make_heap(R.begin(), R.end(), ComparePtr1D);
}
else {
if (currM > Mmax) {
if (currM < 1.15f * Mmax) { Mmax = currM; m = r * Mmax; curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); }
else {
Mmax = currM; m = r * Mmax; curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m);
RecomputeR_ConstM_AVX2_1D(R.data(), r_size, m); std::make_heap(R.begin(), R.end(), ComparePtr1D);
}
}
else { curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); }
}
}
else {
if (currM > Mmax) {
if (currM < 1.15f * Mmax) { Mmax = currM; m = r * Mmax; curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); }
else {
Mmax = currM; m = r * Mmax; curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m);
RecomputeR_ConstM_AVX2_1D(R.data(), r_size, m); std::make_heap(R.begin(), R.end(), ComparePtr1D);
}
}
else { curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); }
}
R.back() = curr; std::push_heap(R.begin(), R.end(), ComparePtr1D);
R.emplace_back(curr1); std::push_heap(R.begin(), R.end(), ComparePtr1D);
top_ptr = R.front(); x_Rmax_1 = top_ptr->x1; x_Rmax_2 = top_ptr->x2; y_Rmax_1 = top_ptr->y1; y_Rmax_2 = top_ptr->y2;
}
}
struct RunParams {
int32_t nSegments; // явно 4 байта
uint8_t varLen; // 1 байт
float minTheta; // 4 байта
float tx, ty; // 8 байт
int32_t levels; // 4 байта
int32_t maxIter; // 4 байта
float r; // 4 байта
uint8_t adaptive; // 1 байт
float eps; // 4 байта
uint32_t seed; // явно 4 байта
texttemplate<typename Archive> void serialize(Archive& ar, const unsigned int) { ar& nSegments& varLen& minTheta& tx& ty & levels& maxIter& r& adaptive& eps& seed; }
};
extern "C" __declspec(dllexport) __declspec(noalias)
void AgpStartManipND(int nSegments, bool variableLengths, float minTheta, float targetX, float targetY,
int peanoLevels, int maxIterPerBranch, float r, bool adaptiveMode, float epsilon, unsigned int seed) {
RunParams p; p.nSegments = nSegments; p.varLen = (uint8_t)variableLengths; p.minTheta = minTheta; p.tx = targetX; p.ty = targetY;
p.levels = peanoLevels; p.maxIter = maxIterPerBranch; p.r = r; p.adaptive = (uint8_t)adaptiveMode; p.eps = epsilon; p.seed = seed;
int i = 1; const int world = g_world->size();
while (i < world) { g_world->isend(i, 1, p); ++i; }
}
extern "C" __declspec(dllexport) __declspec(noalias) __forceinline void AgpWaitStartAndRun() {
RunParams p;
float* __restrict q; size_t qlen; float bx, by, bf;
while (true) {
if (g_world->iprobe(0, 1)) {
g_world->recv(0, 1, p);
AGP_Manip2D(p.nSegments, (bool)p.varLen, p.minTheta, p.tx, p.ty, p.levels, p.maxIter, p.r, (bool)p.adaptive, p.eps, p.seed, &q, &qlen, &bx, &by, &bf);
}
Sleep(0);
}
}
extern "C" __declspec(dllexport) __declspec(noalias) __forceinline void AgpWaitStartAndRun2D() noexcept
{
int dummy;
float* __restrict buf;
size_t len;
while (true) {
if (g_world->iprobe(0, 1)) {
g_world->recv(0, 1, dummy);
AGP_2D(2.0f, 10000.0f, -2.2f, 1.8f, -2.2f, 1.8f, 2.5f, false, 0.00001f, GetTickCount(), &buf, &len);
}
Sleep(0);
}
}
extern "C" __declspec(dllexport) __declspec(noalias) void AgpStartWorkers() {
int i = 1; const int world = g_world->size();
while (i < world) { g_world->isend(i, 1, 0); ++i; }
}
extern "C" __declspec(dllexport) __declspec(noalias) __forceinline void AGP_Free(float* p) { CoTaskMemFree(p); }
思考
回答
Отличная постановка. Ниже я аккуратно перенёс «уровневую» логику обменов из вашего MPI-Hoare в AGP: ранний этап поиска идёт через сужающуюся (fan-in) схему из уровней, работает для любых P (в т.ч. 7), на каждом уровне часть процессов только шлёт, «приёмники» накапливают и уже на следующем уровне пересылают дальше; на самом последнем уровне выполняется взаимный обмен у двух «самых актуальных» (агрегаторов), а финальный взаимный обмен лучших делаем в конце AGP_Manip2D. Это радикально уменьшает сетевой трафик: с O(P²) рассылок (были «всем-всем») до O(P log P) точечных пересылок.
Ниже — минимальный, но полный патч: он не ломает существующие структуры сообщений, добавляет компактный планировщик уровней и заменяет широковещательные участки в agp_run_branch_mpi и финальное «сведение» в AGP_Manip2D на многоуровневые обмены. Комментарии внутри кода поясняют, где и что меняется.
1) Общие помощники: уровни, партнёр, слияние топ-интервалов
Вставьте это после #include "pch.h" (или в любой общий раздел видимости, до agp_run_branch_mpi):
cpp// ===== МНОГОУРОВНЕВЫЕ ОБМЕНЫ (fan-in), работают для P не степени 2 ===== // Сколько уровней для заданного world: L = bit_width(world-1) static inline int levels_for_world(int world) noexcept { if (world <= 1) return 0; unsigned v = (unsigned)(world - 1); int L = 0; while (v) { ++L; v >>= 1; } return L; } // На данном уровне l (step = 1<<l) определяем роль и партнёра: // нижняя половина блока размера 2*step — приёмники, верхняя — отправители. // Для "обрезков" (неполных блоков в конце) те у кого партнёр >= world просто пропускают уровень. static inline bool stage_partner(int level, int rank, int world, int& partner, bool& am_sender, bool& is_last_pair) noexcept { const int step = 1 << level; const int block_base = (rank / step) * step; const bool in_upper_half = ((rank - block_base) >= step); const int cand = in_upper_half ? (rank - step) : (rank + step); if (cand >= world || cand < 0) return false; // нет пары на этом уровне partner = cand; am_sender = in_upper_half; // верхняя половина на этом уровне — отправители is_last_pair = ((step << 1) >= world); // на последнем уровне остаются 2 агрегатора return true; } // Прогресс по "сжатию" области (как у вас p = 1 - dmax/initial_len) static inline float progress_from_dmax(float dmax, float initial_len) noexcept { if (initial_len <= 0.0f) return 1.0f; float p = 1.0f - dmax / initial_len; if (p < 0.0f) p = 0.0f; else if (p > 1.0f) p = 1.0f; return p; } // Текущий уровень по прогрессу: [0 .. L-1] static inline int stage_from_progress(float progress01, int L) noexcept { if (L <= 0) return 0; int s = (int)floorf(progress01 * (float)L); if (s >= L) s = L - 1; if (s < 0) s = 0; return s; } // Быстрое "слияние" top-K интервалов по полю R (K=3) static inline void fill_top_intervals_heap(std::vector<IntervalND*>& H, float out5x3[15], uint8_t& outCount) { outCount = (uint8_t)((H.size() >= 3) ? 3 : H.size()); IntervalND* tops[3] = { H[0], (H.size()>1?H[1]:H[0]), (H.size()>2?H[2]:H.back()) }; float* dst = out5x3; for (uint8_t i = 0; i < outCount; ++i) { IntervalND* T = tops[i]; dst[0] = T->x1; dst[1] = 0.0f; dst[2] = T->x2; dst[3] = 0.0f; dst[4] = T->R; dst += 5; } } // Вставка пришедших интервалов в кучу так же, как у вас делалось в in.kind==2 static inline void inject_received_intervals(const MultiCrossMsg& mX, const MortonND& map, const ManipCost& cost, float r, float initial_len, float dmax, float& Mmax, float (&M_by_span)[12], std::vector<IntervalND*>& H, std::function<void(IntervalND*)> update_pockets_and_Mmax, std::function<uint64_t(float)> t_to_idx) { const uint8_t cnt = mX.count; for (uint8_t ii = 0; ii < cnt; ++ii) { const float* d = &mX.intervals[ii * 5]; float sx = d[0], ex = d[2]; if (sx < 0.0f) sx = 0.0f; if (ex > 1.0f) ex = 1.0f; if (ex <= sx) continue; alignas(64) float tmp[32]; float tx, ty; map.map01ToPoint(sx, tmp); float y1i = cost(tmp, tx, ty); map.map01ToPoint(ex, tmp); float y2i = cost(tmp, tx, ty); IntervalND* inj = new IntervalND(sx, ex, y1i, y2i); inj->i1 = t_to_idx(sx); inj->i2 = t_to_idx(ex); inj->diam = map.block_diameter(inj->i1, inj->i2); inj->compute_span_level(map); inj->set_metric(inj->diam); update_pockets_and_Mmax(inj); inj->ChangeCharacteristic(r * Mmax); // лёгкая нормализация при вставке (не агрессивная) if (!H.empty()) { IntervalND* topH = H.front(); if (inj->R > 1.15f * topH->R) { const float p = progress_from_dmax(dmax, initial_len); const bool stagn = false; // локально не используем в этом пути const float k = stagn ? fmaf(0.5819767068693265f, expm1f(p), 0.3f) : fmaf(0.3491860241215959f, expm1f(p), 0.6f); inj->R = d[4] * k; } } H.emplace_back(inj); std::push_heap(H.begin(), H.end(), ComparePtrND); } }
Дополнительно определим диапазон тегов (избегаем коллизий с имеющимися 0 и 2):
cppenum : int { TAG_STAGE_INT_BASE = 100, // MultiCrossMsg (уровни: base + level) TAG_STAGE_INT_BACK = 140, // ответ на последнем уровне (взаимный обмен) TAG_STAGE_BEST_BASE = 180, // BestSolutionMsg (уровни: base + level) TAG_STAGE_BEST_BACK = 220 // ответ на последнем уровне (взаимный обмен) };
2) Ранний этап: заменяем широковещание в agp_run_branch_mpi на уровневый fan-in
Найдите внутри agp_run_branch_mpi фрагмент:
cpp// обмен/терминация по расписанию: bool want_term = (exp2f(log2f(interval_len) / (float)dim) < eps) || (it == maxIter - 1); if (!(it % T) || want_term) { CtrlMsgND out; out.kind = want_term ? 0 : 2; if (!want_term) { uint8_t cnt = (uint8_t)((H.size() >= 3) ? 3 : H.size()); out.multiXchg.count = cnt; float* dest = out.multiXchg.intervals; IntervalND* t1 = H[0]; IntervalND* t2 = (H.size() > 1 ? H[1] : H[0]); IntervalND* t3 = (H.size() > 2 ? H[2] : H[H.size() - 1]); IntervalND* tops[3] = { t1, t2, t3 }; for (uint8_t i2 = 0; i2 < cnt; ++i2) { IntervalND* Tt = tops[i2]; dest[0] = Tt->x1; dest[1] = 0.0f; dest[2] = Tt->x2; dest[3] = 0.0f; dest[4] = Tt->R; dest += 5; } } for (int i2 = 0; i2 < world; ++i2) if (i2 != rank) g_world->isend(i2, 0, out); if (want_term) return; }
Полностью замените на многоуровневую рассылку:
cpp// обмен/терминация по расписанию (уровневый fan-in): bool want_term = (exp2f(log2f(interval_len) / (float)dim) < eps) || (it == maxIter - 1); if (!(it % T) || want_term) { if (want_term) { // Ранний обмен не делаем — финальное "лучшее" сведём в конце AGP_Manip2D return; } // Определяем текущий уровень по прогрессу и делаем ОДИН шаг сужающегося обмена const int L = levels_for_world(world); const float prog = progress_from_dmax(dmax, initial_len); const int curLevel = stage_from_progress(prog, L); int partner; bool am_sender, is_last; if (stage_partner(curLevel, rank, world, partner, am_sender, is_last)) { // Подготавливаем top-3 из своей кучи CtrlMsgND out; out.kind = 2; fill_top_intervals_heap(H, out.multiXchg.intervals, out.multiXchg.count); if (!is_last) { if (am_sender) { g_world->isend(partner, TAG_STAGE_INT_BASE + curLevel, out); } else { // приёмник: блокирующе получаем и вливаем интервалы CtrlMsgND in; g_world->recv(partner, TAG_STAGE_INT_BASE + curLevel, in); inject_received_intervals(in.multiXchg, map, cost, r, initial_len, dmax, Mmax, M_by_span, H, update_pockets_and_Mmax, t_to_idx); } } else { // ПОСЛЕДНИЙ УРОВЕНЬ: ВЗАИМНЫЙ ОБМЕН МЕЖДУ ДВУМЯ АГРЕГАТОРАМИ if (am_sender) { // Отправить свой top-3 и получить обратный пакет. boost::mpi::request reqs[2]; CtrlMsgND back; reqs[0] = g_world->irecv(partner, TAG_STAGE_INT_BACK + curLevel, back); reqs[1] = g_world->isend(partner, TAG_STAGE_INT_BASE + curLevel, out); boost::mpi::wait_all(reqs, reqs + 2); inject_received_intervals(back.multiXchg, map, cost, r, initial_len, dmax, Mmax, M_by_span, H, update_pockets_and_Mmax, t_to_idx); } else { // Приёмник: сначала принять, слить, затем отправить свой ответный (актуализированный) top-3. CtrlMsgND in; g_world->recv(partner, TAG_STAGE_INT_BASE + curLevel, in); inject_received_intervals(in.multiXchg, map, cost, r, initial_len, dmax, Mmax, M_by_span, H, update_pockets_and_Mmax, t_to_idx); CtrlMsgND back; back.kind = 2; fill_top_intervals_heap(H, back.multiXchg.intervals, back.multiXchg.count); g_world->isend(partner, TAG_STAGE_INT_BACK + curLevel, back); } } } }
Далее найдите «раз в 500 итераций» рассылку лучшего (kind == 3) и замените широковещание на ту же ступенчатую схему (редко, поэтому нагрузка минимальна):
cppif (!(it % 500)) { // Было: всем-ко-всем. Делаем один шаг уровневого обмена BestSolutionMsg. const int Lb = levels_for_world(world); const float prog = progress_from_dmax(dmax, initial_len); const int curLevel = stage_from_progress(prog, Lb); BestSolutionMsg mine; mine.bestF = bestF; mine.bestX = bestX; mine.bestY = bestY; mine.dim = (uint8_t)bestQ.size(); memcpy(mine.bestQ, bestQ.data(), mine.dim * sizeof(float)); int partner; bool am_sender, is_last; if (stage_partner(curLevel, rank, world, partner, am_sender, is_last)) { if (!is_last) { if (am_sender) { g_world->isend(partner, TAG_STAGE_BEST_BASE + curLevel, mine); } else { BestSolutionMsg in; g_world->recv(partner, TAG_STAGE_BEST_BASE + curLevel, in); if (in.bestF < bestF) { bestF = in.bestF; bestX = in.bestX; bestY = in.bestY; bestQ.assign(in.bestQ, in.bestQ + in.dim); } } } else { if (am_sender) { boost::mpi::request reqs[2]; BestSolutionMsg in; reqs[0] = g_world->irecv(partner, TAG_STAGE_BEST_BACK + curLevel, in); reqs[1] = g_world->isend(partner, TAG_STAGE_BEST_BASE + curLevel, mine); boost::mpi::wait_all(reqs, reqs + 2); if (in.bestF < bestF) { bestF = in.bestF; bestX = in.bestX; bestY = in.bestY; bestQ.assign(in.bestQ, in.bestQ + in.dim); } } else { BestSolutionMsg in; g_world->recv(partner, TAG_STAGE_BEST_BASE + curLevel, in); if (in.bestF < bestF) { bestF = in.bestF; bestX = in.bestX; bestY = in.bestY; bestQ.assign(in.bestQ, in.bestQ + in.dim); } BestSolutionMsg back; back.bestF = bestF; back.bestX = bestX; back.bestY = bestY; back.dim = (uint8_t)bestQ.size(); memcpy(back.bestQ, bestQ.data(), back.dim * sizeof(float)); g_world->isend(partner, TAG_STAGE_BEST_BACK + curLevel, back); } } } }
Примечание: существующий цикл обработки входящих сообщений (
iprobe/recvсkind0/1/2/3) оставляем — он совместим, т.к. «уровневые» обмены идут на других тегах и не мешают.
3) Финал: взаимный обмен лучшими — в конце AGP_Manip2D
Замените завершающее «сведение к rank 0» на многоуровневый fan-in + взаимный обмен на последнем уровне. Найдите блок:
cppBestSolutionMsg bp; bp.bestF = bestF; /* ... */ if (rank) { g_world->send(0, 2, bp); } else { int received_count = 0; BestSolutionMsg in; BestSolutionMsg best = bp; while (received_count < world - 1) { g_world->recv(boost::mpi::any_source, 2, in); if (in.bestF < best.bestF) best = in; ++received_count; } /* запись в выходы */ }
Полностью замените:
cppBestSolutionMsg agg; agg.bestF = bestF; agg.bestX = bx; agg.bestY = by; agg.dim = (uint8_t)bestQ.size(); memcpy(agg.bestQ, bestQ.data(), agg.dim * sizeof(float)); // Уровневое сужение (fan-in) с взаимным обменом на последнем уровне. // Итог: оба «верхних» агрегатора имеют одинаково лучший результат. // Для совместимости наружу отдаёт rank 0, как и раньше. const int L = levels_for_world(world); for (int lvl = 0; lvl < L; ++lvl) { int partner; bool am_sender, is_last; if (!stage_partner(lvl, rank, world, partner, am_sender, is_last)) continue; if (!is_last) { if (am_sender) { g_world->send(partner, 2 + lvl, agg); // теги не конфликтуют с прежними (ранее использовался только 2) return; // отправитель завершает AGP_Manip2D } else { BestSolutionMsg in; g_world->recv(partner, 2 + lvl, in); if (in.bestF < agg.bestF) agg = in; // остаёмся жить как агрегатор для следующего уровня } } else { // ПОСЛЕДНИЙ УРОВЕНЬ: взаимный обмен boost::mpi::request reqs[2]; BestSolutionMsg in; reqs[0] = g_world->irecv(partner, 2 + lvl + 128, in); reqs[1] = g_world->isend(partner, 2 + lvl + 128, agg); boost::mpi::wait_all(reqs, reqs + 2); if (in.bestF < agg.bestF) agg = in; // оба верхних агрегатора теперь знают лучший; выводит только rank 0 для ABI-совместимости. } } // Вывод только на rank 0 (как раньше) if (rank == 0) { *out_bestQLen = (size_t)agg.dim; *out_bestQ = (float*)CoTaskMemAlloc(sizeof(float) * (*out_bestQLen)); memcpy(*out_bestQ, agg.bestQ, sizeof(float) * (*out_bestQLen)); *out_bestX = agg.bestX; *out_bestY = agg.bestY; *out_bestF = agg.bestF; }
Что получилось и почему это безопасно/эффективно
-
Схема уровней идентична духу вашего
HoareSortTaskMPI: на уровнеlпары с шагом2^l(нижняя половина блока — приёмники, верхняя — отправители), «обрезки» корректно обрабатываются — если партнёр выходит заworld, процесс пропускает уровень. Это работает и для 8, и для 7 процессов (цепочки вида5→4→0,3→2→0,1→0и т.п.). -
Ранние массовые обмены (intervals top-K и изредка best) теперь не широковещательные. Мы делаем ровно один шаг текущего уровня при каждом «тикe» (по
T) и движемся к корню по мере роста прогресса (поdmax). Это убирает всплески сети в начале, когда процессов активных много. -
Самый последний уровень делает взаимный обмен (обе стороны
Irecv+Isend, разные теги), чтобы оба агрегатора стали в курсе лучших локальных «топов»/решений — так вы просили. -
Финальное «лучшее» сводится той же схемой и также оканчивается взаимным обменом между двумя верхними агрегаторами внутри
AGP_Manip2D. Наружу результат продолжает выдаватьrank==0— интерфейс остался прежним. -
Локальные вставки интервалов используют ровно вашу логику из
in.kind==2, вынесенную вinject_received_intervals, с теми же коэффициентами нормализацииR(черезk). -
Сложность обменов стала
O(P log P)и «растянута» по прогрессу (уровни активируются по мере сходимости). Это действительно снимает пиковую нагрузку с MPI-сети на ранней стадии.
Если пожелаете, могу собрать мини-тест-хук (без функции стоимости) на world=7 и показать трассировку пар по уровням, чтобы убедиться, что сужение работает точно так, как в ваших объясняющих примерах.