logoChat01.ai

extern "C" __declspec(dllexport) __declspec(noalia...

Creado el: 19 de octubre de 2025

Respondido usando GPT-5 Thinking por Chat01

Pregunta

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
{
int schetchick = 0;
const int rank = g_world->rank();
const int partner = rank ^ 1;
int dummy;

text
const float inv_divider = ldexpf(1.0f, -((gActiveMap.levels << 1) + 1));
const float x_addition = (b - a) * inv_divider;
const float y_addition = (d - c) * inv_divider;
const float true_start = a + x_addition;
const float true_end = b - x_addition;

float x_Rmax_1 = true_start;
float x_Rmax_2 = true_end;
const float initial_length = true_end - true_start;
float dmax = initial_length;
const float threshold_03 = 0.3f * initial_length;
const float inv_threshold_03 = 1.0f / threshold_03;
const float start_val = rank ? RastriginFunc(true_end, d - y_addition) : RastriginFunc(true_start, c + y_addition);
float best_f = rank ? RastriginFunc(true_start, d - y_addition) : RastriginFunc(true_end, c + y_addition);
float y_Rmax_1 = start_val;
float y_Rmax_2 = best_f;

std::vector<float, boost::alignment::aligned_allocator<float, 16u>> Extr;
std::vector<Interval* __restrict, boost::alignment::aligned_allocator<Interval* __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 Interval(true_start, true_end, start_val, best_f, 2.0f));

float Mmax = R.front()->M;
float m = r * Mmax;

while (true) {
	const float cooling = ldexpf(1.0f, -(1.0f / 138.63f) * ++schetchick);
	const int T = static_cast<int>(fmaf(20.0f, cooling, 10.0f));
	const float k = fmaf(0.2f, cooling, 0.7f);

	const Interval* __restrict top_ptr = R.front();
	if (schetchick % T == 0 || g_world->iprobe(partner, 0)) {
		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);

		const CrossMsg outbound = CrossMsg{ s_x1, s_x2, e_x1, e_x2, top_ptr->R };
		CrossMsg inbound;
		g_world->sendrecv(partner, 0, outbound, partner, 0, inbound);

		const float sx = FindX2D_analytic(inbound.s_x1, inbound.s_x2);
		const float ex = FindX2D_analytic(inbound.e_x1, inbound.e_x2);

		Interval* const __restrict injected = new Interval(sx, ex, RastriginFunc(inbound.s_x1, inbound.s_x2), RastriginFunc(inbound.e_x1, inbound.e_x2), 2.0f);
		injected->R = inbound.Rtop * k;
		R.emplace_back(injected);
		std::push_heap(R.begin(), R.end(), ComparePtr);
	}

	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);
	}

	std::pop_heap(R.begin(), R.end(), ComparePtr);
	Interval* const __restrict promejutochny_otrezok = R.back();

	const float segment_x1 = promejutochny_otrezok->x1;
	const float segment_x2 = promejutochny_otrezok->x2;
	const float len2 = segment_x2 - new_point;
	const float len1 = new_point - segment_x1;
	const float interval_len = len1 < len2 ? len1 : len2;

	if (interval_len < epsilon || schetchick == static_cast<int>(global_iterations) || g_world->iprobe(partner, 2)) {
		if (partner) {
			Extr.emplace_back(static_cast<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));
		}
		g_world->sendrecv(partner, 2, 0, partner, 2, dummy);
		return;
	}

	Interval* const __restrict curr = new Interval(segment_x1, new_point, promejutochny_otrezok->y1, new_value, 2.0f);
	Interval* const __restrict curr1 = new Interval(new_point, segment_x2, new_value, promejutochny_otrezok->y2, 2.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;
			size_t i = 0u;

#pragma loop(ivdep)
while (i < r_size) {
const float len_item = R[i]->x2 - R[i]->x1;
if (len_item > dmax) dmax = len_item;
++i;
}
}

text
		if (threshold_03 > dmax && schetchick % 3 == 0 || 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 = fmaf(progress, progress, 1.0f);
			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));

			if (r_size < 64u) {
				size_t i = 0u;

#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, R[i]->x2 - R[i]->x1, R[i]->M * alpha));
++i;
}
}
else {
RecomputeR_AffineM_AVX2(R.data(), r_size, GLOBAL_FACTOR, alpha);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
else {
if (currM > Mmax) {
if (currM - Mmax < Mmax * 0.15f) {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
else {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
if (r_size < 64u) {
size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(m);
++i;
}
}
else {
RecomputeR_ConstM_AVX2(R.data(), r_size, m);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}
}
else {
if (currM > Mmax) {
if (currM - Mmax < Mmax * 0.15f) {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
else {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
if (r_size < 64u) {
size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(m);
++i;
}
}
else {
RecomputeR_ConstM_AVX2(R.data(), r_size, m);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}

text
	R.back() = curr;
	std::push_heap(R.begin(), R.end(), ComparePtr);
	R.emplace_back(curr1);
	std::push_heap(R.begin(), R.end(), ComparePtr);

	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;
}

} - рассмотри этот код на предмет гонок данных и возможных дэдлоков - я ожидаю поведение что если какой-то процесс пришёл к фазе обмена или к фазе завершения раньше - то другой процесс дойдя до условия обмена или условия завершения - даже если эти условия не выполнены - зайдёт в эти фазы так как сработает iprobe так как iprobe проверяет сообщения в очереди а sendrecv добавляет сообщения в очередь - так ли это работает? будь особенно внимателен и проверь все возможные случаи приводящие к гонкам данных или дэдлоку и скажи будет ли ожидаемое поведение - думай дольше для лучшего ответа, ответь на русском, я использую здесь sendrecv только как синхронизацию перед завершением - по идее можно использовать неблокирующую отправку если в условие выполнен вход не через iprobe с последующим ожиданием объекта request для этого случая и блокирующий приём для случаев когда вход в условия выполнен через iprobe - НО ЭТО ЛОМАЕТ ЛОГИКУ ЕСЛИ ОБА ПРОЦЕССА ОДНОВРЕМЕННО ДОСТИГНУТ УСЛОВИЯ ОБМЕНА ИЛИ ЗАВЕРШЕНИЯ И ОБА НАЧНУТ ОТПРАВКУ - ПОЭТОМУ ИСПОЛЬЗУЮ sendrecv для возможности и синхронного обмена при одновременном попадании в эти условия, проблема в том - что в теории это должно работать всегда, но на практике когда я запускаю приложение использующее эту функцию то несколько нажатий кнопки обрабатываются - но на некотором нажатии происходит дэдлок - приложение перестаёт отвечать а после закрытия неотвечающего приложения получаю: C:\Users\maxim>mpiexec -n 2 D:\TEST_AGP\x64\Release\TEST_AGP.exe

job aborted:
[ranks] message

[0] process exited without calling finalize

[1] terminated

---- error analysis -----

[0] on DESKTOP-C1VT6JF
D:\TEST_AGP\x64\Release\TEST_AGP.exe ended prematurely and may have crashed. exit code 0xcfffffff

---- error analysis ----- - также проанализируй код на предмет других сомнительных моментов - проверь например такие мелочи как чрезмерное использование const и укажи места где это не просто бессмысленно а может даже навредить, также внимательно проверь макрос для векторизованного вычисления синуса и косинуса - корректно ли он реализован - совпадает ли там реализация со скалярной использующей аппроксимацию полиномами по смыслу и проверь другие сомнительные моменты в коде и дай рекомендации, полный код dll: #include "pch.h"

#define XOR_RAND(state, result_var)
do {
int s = state;
s ^= s << 13;
s ^= s >> 17;
s ^= s << 5;
state = s;
result_var = state * 0x1.0p-32f;
} while(0)

#define XOR_RAND_GRSH(state, result_var)
do {
int s = state;
s ^= s << 13;
s ^= s >> 17;
s ^= s << 5;
state = s;
result_var = fmaf(state, 0x1.0p-31f, -1.0f);
} while(0)

#define FABE13_COS(x, result_var)
do {
const float abs_val = fabsf(x);
float reduced = fmodf(abs_val, 6.28318530718f);
if(reduced > 3.14159265359f) {
reduced = 6.28318530718f - reduced;
}
if(reduced < 1.57079632679f) {
const float val2 = reduced * reduced;
const float val4 = val2 * val2;
result_var = fmaf(val4, fmaf(val2, -0.0013888889f, 0.0416666667f), fmaf(val2, -0.5f, 1.0f));
} else {
reduced = 3.14159265359f - reduced;
const float val2 = reduced * reduced;
const float val4 = val2 * val2;
result_var = -fmaf(val4, fmaf(val2, -0.0013888889f, 0.0416666667f), fmaf(val2, -0.5f, 1.0f));
}
} 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)

static __declspec(noalias) __forceinline float ShekelFunc(const float x, const float seed) noexcept
{
int i = 0;
float current_state = seed, current_res, current_res2, res = 0.0f;
while (i < 10) {
XOR_RAND(current_state, current_res);
const float x_part = fmaf(-current_res, 10.0f, x);
XOR_RAND(current_state, current_res);
XOR_RAND(current_state, current_res2);
float delimiter = fmaf(fmaf(current_res, 20.0f, 5.0f), x_part * x_part, fmaf(current_res2, 0.2f, 1.0f));
delimiter = copysignf(fmaxf(fabsf(delimiter), FLT_MIN), delimiter);
res -= 1.0f / delimiter;
++i;
}
return res;
}

static __declspec(noalias) __forceinline float RastriginFunc(const float x1, const float x2) noexcept
{
const float term1 = fmaf(x1, x1, x2 * x2);
float cos1, cos2;
FABE13_COS(6.28318530717958647692f * x1, cos1);
FABE13_COS(6.28318530717958647692f * x2, cos2);
return (term1 - fmaf(cos1 + cos2, 10.0f, -14.6f)) * fmaf(-term1, 0.25f, 18.42f);
}

static __declspec(noalias) __forceinline float HillFunc(const float x, const float seed) noexcept
{
int j = 0;
__declspec(align(32)) float angles[14u];
const float start_angle = 6.28318530717958647692f * x;
#pragma loop(ivdep)
while (j < 14) {
angles[j] = start_angle * static_cast<float>(j + 1);
++j;
}
__declspec(align(32)) float sin_vals[14u];
__declspec(align(32)) float cos_vals[14u];
FABE13_SINCOS(angles, sin_vals, cos_vals, 14u);
float current_state = seed, current_res, current_res2;
XOR_RAND(current_state, current_res);
float res = fmaf(current_res, 2.0f, -1.1f);
--j;
while (j >= 0) {
XOR_RAND(current_state, current_res);
XOR_RAND(current_state, current_res2);
res += fmaf(fmaf(current_res, 2.0f, -1.1f), sin_vals[j], fmaf(current_res2, 2.0f, -1.1f) * cos_vals[j]);
--j;
}
return res;
}

static __declspec(noalias) __forceinline float GrishaginFunc(const float x1, const float x2, const float seed) noexcept
{
int j = 0;
__declspec(align(32)) float angles_j[8u];
__declspec(align(32)) float angles_k[8u];
#pragma loop(ivdep)
while (j < 8) {
const float pj_mult = 3.14159265358979323846f * static_cast<float>(j + 1);
angles_j[j] = pj_mult * x1;
angles_k[j] = pj_mult * x2;
++j;
}
__declspec(align(32)) float sin_j[8u], cos_j[8u];
__declspec(align(32)) float sin_k[8u], cos_k[8u];
FABE13_SINCOS(angles_j, sin_j, cos_j, 8u);
FABE13_SINCOS(angles_k, sin_k, cos_k, 8u);
--j;
float part1 = 0.0f;
float part2 = 0.0f;
float current_state = seed, current_res, current_res2;
while (j >= 0) {
size_t k = 0u;
while (k < 8u) {
const float sin_term = sin_j[j] * sin_j[j];
const float cos_term = cos_k[k] * cos_k[k];
XOR_RAND_GRSH(current_state, current_res);
XOR_RAND_GRSH(current_state, current_res2);
part1 = fmaf(current_res, sin_term, fmaf(current_res2, cos_term, part1));
XOR_RAND_GRSH(current_state, current_res);
XOR_RAND_GRSH(current_state, current_res2);
part2 = fmaf(-current_res, cos_term, fmaf(current_res2, sin_term, part2));
++k;
}
--j;
}
return -sqrtf(fmaf(part1, part1, part2 * part2));
}

static __declspec(noalias) __forceinline float Shag(const float _m, const float x1, const float x2, const float y1,
const float y2, const float _N, const float _r) noexcept
{
const float diff = y2 - y1;
const float sign_mult = _mm_cvtss_f32(_mm_castsi128_ps(_mm_set1_epi32(
0x3F800000u | ((reinterpret_cast<const uint32_t>(&diff) & 0x80000000u) ^ 0x80000000u))));
return _N == 1.0f
? fmaf(-(1.0f / _m), diff, x1 + x2) * 0.5f
: _N == 2.0f
? fmaf(sign_mult / (_m * _m), diff * diff * _r, x1 + x2) * 0.5f
: fmaf(sign_mult / powf(_m, _N), powf(diff, _N) * _r, x1 + x2) * 0.5f;
}

__declspec(align(16)) struct Slab final {
char* const __restrict base;
char* __restrict current;
char* const __restrict end;

text
__declspec(noalias) __forceinline Slab(void* const __restrict memory, const size_t usable_size) noexcept
	: base(static_cast<char* __restrict>(memory))
	, current(base)
	, end(base + (usable_size & ~static_cast<size_t>(63u)))
{
}

};

static tbb::enumerable_thread_specific<Slab*> tls( noexcept {
void* const __restrict memory = _aligned_malloc(16777216u, 16u);
Slab* const __restrict slab = static_cast<Slab*>(_aligned_malloc(32u, 16u));
new (slab) Slab(memory, 16777216u);

char* __restrict p = slab->base;

#pragma loop(ivdep)
while (p < slab->end) {
*p = 0u;
p += 4096u;
}
return slab;
}());

__declspec(align(64)) struct Interval final {
const float x1;
const float x2;
const float y1;
const float y2;
const float delta_y;
const float ordinate_factor;
const float N_factor;
const float quadratic_term;
const float M;
float R;

text
__declspec(noalias) __forceinline void* operator new(const size_t) noexcept {
	Slab* const __restrict s = tls.local();
	char* const __restrict result = s->current;
	s->current += 64u;
	return result;
}

__declspec(noalias) __forceinline Interval(const float _x1, const float _x2,
	const float _y1, const float _y2, const 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 : _N == 2.0f ? sqrtf(_x2 - _x1) : powf(_x2 - _x1, 1.0f / _N))
	, quadratic_term((1.0f / N_factor)* delta_y* delta_y)
	, M((1.0f / N_factor)* fabsf(delta_y))
{
}

__declspec(noalias) __forceinline void ChangeCharacteristic(const float _m) noexcept
{
	R = fmaf(1.0f / _m, quadratic_term, fmaf(_m, N_factor, ordinate_factor));
}

};

static __declspec(noalias) __forceinline bool ComparePtr(const Interval* const __restrict a, const Interval* const __restrict b) noexcept
{
return a->R < b->R;
}

const enum List : uint8_t { Top = 0b00u, Down = 0b01u, Left = 0b10u, Right = 0b11u };

__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;

text
__declspec(noalias) __forceinline Peano2DMap(
	const int L,
	const float _a,
	const float _b,
	const float _c,
	const float _d,
	const uint8_t startType
) noexcept
	: levels(L)
	, a(_a), b(_b), c(_c), d(_d)
	, lenx(_b - _a)
	, leny(_d - _c)
	, inv_lenx(1.0f / (_b - _a))
	, scale(static_cast<uint32_t>(1u) << (L << 1))
	, start(startType)
{
}

};

__declspec(align(4)) struct Step final {
const uint8_t next;
const uint8_t dx;
const uint8_t dy;
};

__declspec(align(4)) struct InvStep final {
const uint8_t q;
const uint8_t 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 Peano2DMap gActiveMap(0, 0.0f, 0.0f, 0.0f, 0.0f, 0b00u);

__declspec(noalias) __forceinline void HitTest2D_analytic(const float x_param, float& out_x1, float& out_x2) noexcept
{
const float a = gActiveMap.a;
const float inv_lenx = gActiveMap.inv_lenx;
const uint32_t scale = gActiveMap.scale;
const uint32_t scale_minus_1 = scale - 1u;
const float lenx = gActiveMap.lenx;
const float leny = gActiveMap.leny;
const float c = gActiveMap.c;
const uint8_t start = gActiveMap.start;
const int levels = gActiveMap.levels;

text
float norm = (x_param - a) * inv_lenx;
norm = fminf(fmaxf(norm, 0.0f), 0x1.fffffep-1f);

uint32_t idx = static_cast<uint32_t>(norm * static_cast<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;

#pragma loop(ivdep)
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;
}

__declspec(noalias) __forceinline float FindX2D_analytic(const float px, const float py) noexcept
{
const float a = gActiveMap.a;
const float b = gActiveMap.b;
const float c = gActiveMap.c;
const float d = gActiveMap.d;
const float lenx = gActiveMap.lenx;
const float leny = gActiveMap.leny;
const uint32_t scale = gActiveMap.scale;
const uint8_t start = gActiveMap.start;
const int levels = gActiveMap.levels;

text
const float clamped_px = fminf(fmaxf(px, a), b);
const float 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;

#pragma loop(ivdep)
while (l < levels) {
sx *= 0.5f; sy *= 0.5f;
const float mx = x0 + sx;
const float my = y0 + sy;

text
	const uint32_t tr = static_cast<uint32_t>((clamped_px > mx) & (clamped_py > my));
	const uint32_t tl = static_cast<uint32_t>((clamped_px < mx) & (clamped_py > my));
	const uint32_t dl = static_cast<uint32_t>((clamped_px < mx) & (clamped_py < my));
	const uint32_t none = static_cast<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;
	const uint32_t dy = dd & 1u;
	x0 += dx ? sx : 0.0f;
	y0 += dy ? sy : 0.0f;
	++l;
}

const float scale_reciprocal = 1.0f / static_cast<float>(scale);
return fmaf(static_cast<float>(idx) * scale_reciprocal, lenx, a);

}

static const boost::mpi::environment* __restrict g_env;
static const boost::mpi::communicator* __restrict g_world;

extern "C" __declspec(dllexport) __declspec(noalias) __forceinline int AgpInit(const int peanoLevel, const float a, const float b, const float c, const float d) noexcept
{
g_env = new boost::mpi::environment();
g_world = new boost::mpi::communicator();
const int rank = g_world->rank();

text
new(&gActiveMap) Peano2DMap(peanoLevel, a, b, c, d, rank ? static_cast<uint8_t>(Down) : static_cast<uint8_t>(Top));
return rank;

}

__declspec(align(16)) struct CrossMsg final {
float s_x1, s_x2;
float e_x1, e_x2;
float Rtop;
template<typename Archive>
__declspec(noalias) __forceinline void serialize(Archive& const __restrict ar, const unsigned int) noexcept { ar& s_x1& s_x2& e_x1& e_x2& Rtop; }
};

static __declspec(noalias) __forceinline void RecomputeR_ConstM_AVX2(Interval* const* const __restrict arr, const size_t n, const float m) noexcept {
const __m256 vm = _mm256_set1_ps(m);

text
__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;
const int limit = static_cast<int>(n & ~7u);
if (n >= 8u) {
	while (i < static_cast<size_t>(limit)) {
		__declspec(align(32)) float q[8], nf[8], ord[8];
		int k = 0;

#pragma loop(ivdep)
while (k < 8) {
const Interval* const __restrict p = arr[i + k];
q[k] = p->quadratic_term;
nf[k] = p->N_factor;
ord[k] = p->ordinate_factor;
++k;
}
const __m256 vq = _mm256_load_ps(q);
const __m256 vnf = _mm256_load_ps(nf);
const __m256 vod = _mm256_load_ps(ord);

text
		const __m256 t = _mm256_fmadd_ps(vm, vnf, vod);
		const __m256 res = _mm256_fmadd_ps(vq, vinvm, t);

		__declspec(align(32)) float out[8];
		_mm256_store_ps(out, res);
		k = 0;

#pragma loop(ivdep)
while (k < 8) {
arr[i + k]->R = out[k];
++k;
}
i += 8u;
}
}
while (i < n) {
arr[i]->ChangeCharacteristic(m);
++i;
}
}

static __declspec(noalias) __forceinline void RecomputeR_AffineM_AVX2(Interval* const* const __restrict arr, const size_t n, const float GLOBAL_FACTOR, const float alpha) noexcept {
const __m256 vGF = _mm256_set1_ps(GLOBAL_FACTOR);
const __m256 va = _mm256_set1_ps(alpha);

text
size_t i = 0;
const int limit = static_cast<int>(n & ~7u);
if (n >= 8u) {
	while (i < static_cast<size_t>(limit)) {
		__declspec(align(32)) float len[8], Mv[8], q[8], nf[8], ord[8];
		int k = 0;

#pragma loop(ivdep)
while (k < 8) {
const Interval* const p = arr[i + k];
len[k] = p->x2 - p->x1;
Mv[k] = p->M;
q[k] = p->quadratic_term;
nf[k] = p->N_factor;
ord[k] = p->ordinate_factor;
++k;
}
const __m256 vlen = _mm256_load_ps(len);
const __m256 vM = _mm256_load_ps(Mv);
const __m256 vq = _mm256_load_ps(q);
const __m256 vnf = _mm256_load_ps(nf);
const __m256 vod = _mm256_load_ps(ord);

text
		const __m256 vm = _mm256_fmadd_ps(vGF, vlen, _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);

		__declspec(align(32)) float out[8];
		_mm256_store_ps(out, res);
		k = 0;

#pragma loop(ivdep)
while (k < 8) {
arr[i + k]->R = out[k];
++k;
}
i += 8u;
}
}
while (i < n) {
const Interval* const __restrict p = arr[i];
const float mi = fmaf(GLOBAL_FACTOR, 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));
++i;
}
}

extern "C" __declspec(dllexport) __declspec(noalias)
void AGP_1D(const float global_iterations,
const float a, const float b, const float r,
const bool mode, const float epsilon, const float seed,
float** const __restrict out_data, size_t* const __restrict out_len) noexcept
{
int schetchick = 0;
const float initial_length = b - a;
float dmax = initial_length;
const float threshold_03 = 0.3f * initial_length;
const float 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;
float x_Rmax_2 = b;
float y_Rmax_1 = start_val;
float y_Rmax_2 = best_f;

text
std::vector<float, boost::alignment::aligned_allocator<float, 16u>> Extr;
std::vector<Interval* __restrict, boost::alignment::aligned_allocator<Interval* __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 Interval(a, b, start_val, best_f, 1.0f));

float Mmax = R.front()->M;
float m = r * Mmax;

while (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(), ComparePtr);
	const Interval* const __restrict promejutochny_otrezok = R.back();

	const float new_x1 = promejutochny_otrezok->x1;
	const float new_x2 = promejutochny_otrezok->x2;
	const float len2 = new_x2 - new_point;
	const float len1 = new_point - new_x1;
	const float interval_len = len1 < len2 ? len1 : len2;

	if (interval_len < epsilon || ++schetchick == static_cast<int>(global_iterations)) {
		Extr.emplace_back(static_cast<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;
	}

	Interval* const __restrict curr = new Interval(new_x1, new_point, promejutochny_otrezok->y1, new_value, 1.0f);
	Interval* const __restrict curr1 = new Interval(new_point, new_x2, new_value, promejutochny_otrezok->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;
			size_t i = 0u;

#pragma loop(ivdep)
while (i < r_size) {
const float len_item = R[i]->x2 - R[i]->x1;
if (len_item > dmax) dmax = len_item;
++i;
}
}

text
		if (threshold_03 > dmax && schetchick % 3 == 0 || 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 = fmaf(progress, progress, 1.0f);
			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));

			if (r_size < 64u) {
				size_t i = 0u;

#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, R[i]->x2 - R[i]->x1, R[i]->M * alpha));
++i;
}
}
else {
RecomputeR_AffineM_AVX2(R.data(), r_size, GLOBAL_FACTOR, alpha);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
else {
if (currM > Mmax) {
if (currM - Mmax < Mmax * 0.15f) {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
else {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
if (r_size < 64u) {
size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(m);
++i;
}
}
else {
RecomputeR_ConstM_AVX2(R.data(), r_size, m);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}
}
else {
if (currM > Mmax) {
if (currM - Mmax < Mmax * 0.15f) {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
else {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
if (r_size < 64u) {
size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(m);
++i;
}
}
else {
RecomputeR_ConstM_AVX2(R.data(), r_size, m);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}

text
	R.back() = curr;
	std::push_heap(R.begin(), R.end(), ComparePtr);
	R.emplace_back(curr1);
	std::push_heap(R.begin(), R.end(), ComparePtr);

	const Interval* const __restrict 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;
}

}

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
{
int schetchick = 0;
const int rank = g_world->rank();
const int partner = rank ^ 1;
int dummy;

text
const float inv_divider = ldexpf(1.0f, -((gActiveMap.levels << 1) + 1));
const float x_addition = (b - a) * inv_divider;
const float y_addition = (d - c) * inv_divider;
const float true_start = a + x_addition;
const float true_end = b - x_addition;

float x_Rmax_1 = true_start;
float x_Rmax_2 = true_end;
const float initial_length = true_end - true_start;
float dmax = initial_length;
const float threshold_03 = 0.3f * initial_length;
const float inv_threshold_03 = 1.0f / threshold_03;
const float start_val = rank ? RastriginFunc(true_end, d - y_addition) : RastriginFunc(true_start, c + y_addition);
float best_f = rank ? RastriginFunc(true_start, d - y_addition) : RastriginFunc(true_end, c + y_addition);
float y_Rmax_1 = start_val;
float y_Rmax_2 = best_f;

std::vector<float, boost::alignment::aligned_allocator<float, 16u>> Extr;
std::vector<Interval* __restrict, boost::alignment::aligned_allocator<Interval* __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 Interval(true_start, true_end, start_val, best_f, 2.0f));

float Mmax = R.front()->M;
float m = r * Mmax;

while (true) {
	const float cooling = ldexpf(1.0f, -(1.0f / 138.63f) * ++schetchick);
	const int T = static_cast<int>(fmaf(20.0f, cooling, 10.0f));
	const float k = fmaf(0.2f, cooling, 0.7f);

	const Interval* __restrict top_ptr = R.front();
	//if (false) {
	if (schetchick % T == 0 || g_world->iprobe(partner, 0)) {
		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);

		const CrossMsg outbound = CrossMsg{ s_x1, s_x2, e_x1, e_x2, top_ptr->R };
		CrossMsg inbound;
		g_world->sendrecv(partner, 0, outbound, partner, 0, inbound);

		const float sx = FindX2D_analytic(inbound.s_x1, inbound.s_x2);
		const float ex = FindX2D_analytic(inbound.e_x1, inbound.e_x2);

		Interval* const __restrict injected = new Interval(sx, ex, RastriginFunc(inbound.s_x1, inbound.s_x2), RastriginFunc(inbound.e_x1, inbound.e_x2), 2.0f);
		injected->R = inbound.Rtop * k;
		R.emplace_back(injected);
		std::push_heap(R.begin(), R.end(), ComparePtr);
	}

	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);
	}

	std::pop_heap(R.begin(), R.end(), ComparePtr);
	Interval* const __restrict promejutochny_otrezok = R.back();

	const float segment_x1 = promejutochny_otrezok->x1;
	const float segment_x2 = promejutochny_otrezok->x2;
	const float len2 = segment_x2 - new_point;
	const float len1 = new_point - segment_x1;
	const float interval_len = len1 < len2 ? len1 : len2;

	if (interval_len < epsilon || schetchick == static_cast<int>(global_iterations) || g_world->iprobe(partner, 2)) {
		if (partner) {
			Extr.emplace_back(static_cast<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));
		}
		g_world->sendrecv(partner, 2, 0, partner, 2, dummy);
		return;
	}

	Interval* const __restrict curr = new Interval(segment_x1, new_point, promejutochny_otrezok->y1, new_value, 2.0f);
	Interval* const __restrict curr1 = new Interval(new_point, segment_x2, new_value, promejutochny_otrezok->y2, 2.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;
			size_t i = 0u;

#pragma loop(ivdep)
while (i < r_size) {
const float len_item = R[i]->x2 - R[i]->x1;
if (len_item > dmax) dmax = len_item;
++i;
}
}

text
		if (threshold_03 > dmax && schetchick % 3 == 0 || 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 = fmaf(progress, progress, 1.0f);
			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));

			if (r_size < 64u) {
				size_t i = 0u;

#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, R[i]->x2 - R[i]->x1, R[i]->M * alpha));
++i;
}
}
else {
RecomputeR_AffineM_AVX2(R.data(), r_size, GLOBAL_FACTOR, alpha);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
else {
if (currM > Mmax) {
if (currM - Mmax < Mmax * 0.15f) {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
else {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
if (r_size < 64u) {
size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(m);
++i;
}
}
else {
RecomputeR_ConstM_AVX2(R.data(), r_size, m);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}
}
else {
if (currM > Mmax) {
if (currM - Mmax < Mmax * 0.15f) {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
else {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
if (r_size < 64u) {
size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(m);
++i;
}
}
else {
RecomputeR_ConstM_AVX2(R.data(), r_size, m);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}

text
	R.back() = curr;
	std::push_heap(R.begin(), R.end(), ComparePtr);
	R.emplace_back(curr1);
	std::push_heap(R.begin(), R.end(), ComparePtr);

	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;
}

}

extern "C" __declspec(dllexport) __declspec(noalias) __forceinline void AgpWaitStartAndRun() noexcept
{
int dummy;
float* __restrict buf;
size_t len;
while (true) {
if (g_world->iprobe(0, 3)) {
g_world->recv(0, 3, dummy);
AGP_2D(2.0f, 10000.0f, -2.2f, 1.8f, -2.2f, 1.8f, 2.5f, false, 0.00001f, GetTickCount(), &buf, &len);
}
}
}

extern "C" __declspec(dllexport) __declspec(noalias) __forceinline void AgpStartWorkers() noexcept
{
g_world->isend(1, 3, 0);
}

extern "C" __declspec(dllexport) __declspec(noalias) __forceinline void AGP_Free(float* const __restrict p) noexcept
{
CoTaskMemFree(p);
}

Pensamiento

Analyzing MPI-like behavior

The user wants me to analyze potential data races and deadlocks when using Boost.MPI with g_world->iprobe and sendrecv. They expect iprobe to detect messages queued by sendrecv, allowing processes to proceed with exchange/termination phases. I'll verify whether sendrecv correctly queues messages for each process with different tags.

Clarifying matching behavior

The semantics of iprobe and sendrecv are subtle. iprobe checks for incoming messages, while sendrecv is blocking, performing both send and receive actions. In Boost.MPI, sendrecv is collective, implemented with MPI_Sendrecv. If one process calls it and the other doesn't, the call will block until both are matched. When the other process later checks with iprobe, it can detect the message.

Understanding MPI_Sendrecv blocking

The first process in sendrecv might be blocked if it's still sending, since the send operation depends on the send mode. The MPI_Sendrecv call handles both sending and receiving. If the send is in standard mode, it may complete without the receiver posting a receive. However, iprobe only detects the message once the send is initiated and delivered. Generally, sendrecv posts the receive first and then the send, securing completion once both are posted.

Understanding sendrecv matching

For sendrecv to progress, both processes need to enter the call. If one does and the other doesn't, the first process could block waiting for a matching receive. Even if the receive is posted, without the partner sending, the operation won't complete. In this case, iprobe detects the partner's message, but the challenge is how the partner's sendrecv gets initiated.

Addressing sendrecv ordering and race conditions

The main risk is deadlock when two processes attempt different actions due to tag mismatches. For example, if one process is stuck in the exchange while the other is in termination, they can both block. The issue comes from tags and their associated iprobe checks, which need to be carefully ordered to avoid mismatches. Reordering logic, like checking termination first, can minimize this deadlock.

Analyzing memory and concurrency bugs

The main issues seem to be poor memory management and concurrency bugs. Memory isn't freed and slabs grow until they overwrite bounds, causing possible crashes. Also, there’s improper use of __restrict on pointers and invalidations of references to vector elements after modifications. MPI communication can miss termination signals because of mismatches in tags. Additionally, a lack of sync between sendrecv operations and iprobe could lead to backlogs.

Identifying deadlock and race conditions

There are multiple issues causing potential deadlocks and race conditions. The misuse of the partner variable leads to asymmetric output writing. The handling of memory pointers could be unsafe for the non-master rank. MPI communication may cause a mismatch and deadlock due to improper sync at termination. I’m considering fixes like checking termination first or organizing sendrecv calls to ensure tag consistency.

Reviewing trigonometric macros

The first macro, FABE13_COS, seems correct for computing cosine with reflection and reduction techniques, but boundary cases near π could involve rounding errors. For FABE13_SINCOS, the approach is mostly plausible, reflecting both sine and cosine properly across quadrants, though handling of sign flips and amplitude for larger values feels inconsistent. The shift in logic between scalar and vector paths might introduce discrepancies when near boundary values, but otherwise looks solid.

Identifying potential macro issues

Using static aligned constants in macros may lead to redefinition problems if the macro is included multiple times. The use of unaligned memory in AVX instructions could also cause crashes if input data isn't properly aligned, though array alignment is generally guaranteed. A small issue arises in sign reflection for sines, but overall, the logic seems consistent. Additionally, the missing MPI_Finalize call from an untracked environment pointer could lead to unexpected behavior, so ensuring proper memory cleanup sounds essential.

Addressing MPI deadlocks and code hazards

The worker function lacks an exit condition, causing MPI-related shutdown issues. Issues include unfreed memory and potential deadlocks due to MPI_Sendrecv and it's reliance on iprobe. Misalignments with message tags and lack of memory finalization may cause inconsistency, deadlock, or crashes. Suggestions include better synchronization with nonblocking calls, tag unification, and an improved handshake protocol for exchange operations. Further analysis on memory management is advised.

Improving MPI termination and constant misuse

To prevent deadlocks, both sides must regularly perform non-blocking probes and call send-receive with matching tags. Unifying payload sizes reduces confusion. On the misuse of const, it's generally fine, but using const incorrectly on references like Archive& const might confuse without impact. Also, __restrict on vector types may be misleading. Pointer conditions and global pointers need careful thread management.

Analyzing MPI termination and structuring fix

To prevent deadlocks, unify control tags for both exchange and termination. Ensure both sides always call sendrecv to avoid mismatches. Consider using MPI_Probe/MPI_Recv style for synchronization. Order operations carefully, such as using a barrier after posting termination requests to avoid issues with exchanges.

Analyzing MPI communication and potential bugs

To avoid deadlocks, we should first check for pending messages before termination, ensuring matching send/receive behavior. There are bugs in the XOR_RAND macro where float-to-int conversion leads to resolution loss and breaks RNG behavior. Also, some misuse of __restrict with std::vector is non-standard. These issues can affect algorithm correctness, leading to potential bugs in the RNG.

Investigating deadlock and memory issues

There are several potential issues: using ldexpf with a non-integer exponent causes unexpected behavior (piecewise constant steps), and missing MPI_Finalize() can lead to memory corruption. The code also mishandles thread-safety with static variables and message exchanges, especially with communication tags, which needs careful revision to prevent deadlock.

Thought for 5m 2s

Respuesta

Короткий вывод заранее:

Ожидаемая вами логика «один процесc зашёл в фазу обмена/завершения → второй, даже если его локальные условия ещё не выполнены, всё-равно зайдёт туда, потому что iprobe увидит сообщение, которое поставил в очередь sendrecv» не гарантируется MPI. Из-за этого в вашем коде реально возможны взаимные блокировки. Плюс у вас есть как минимум одна тяжёлая логическая/ресурсная ошибка, которая даёт “иногда работает, иногда зависает/падает”: локальный аренный аллокатор не сбрасывается между запусками, и память выходит за границы «слайда».

Ниже подробно — где именно залипает, почему, и как исправить; плюс ревью остального кода.

1) Почему iprobe + sendrecv у вас периодически клинит

Важные факты про MPI (актуальны и для Boost.MPI):

  1. MPI_Sendrecv (Boost.MPI sendrecv) — блокирующий вызов, который завершится только когда и отправка, и приём завершены. Он не «кладёт заранее сообщение в очередь» в смысле, который бы всегда гарантировал положительный iprobe на другой стороне.

  2. Стандартный режим отправки (MPI_Send/часть Sendrecv) может буферизовать маленькие сообщения (eager), но не обязан. Если буфера нет/не хватает — отправка блокируется до появления соответствующего приёма. Следовательно, другой процесс может не увидеть ничего в iprobe, пока вы висите внутри sendrecv.

  3. iprobe(src, tag) видит только фактически доставленные в «неожиданные» (unexpected) сообщения, и строго по тегу. Если партнёр сидит в sendrecv с другим тегом или его send не прошёл в eager, ваш iprobe будет возвращать false.

Конкретные дедлок-сценарии в вашем цикле AGP_2D

У вас два разных «канала» по тегам:

  • tag 0 — «обмен интервалами» (ветка в начале итерации: schetchick % T == 0 || iprobe(partner, 0) → sendrecv(..., tag=0)).
  • tag 2 — «завершение» (ветка ближе к концу итерации: if (interval_len < epsilon || ... || iprobe(partner, 2)) { ... sendrecv(..., tag=2); return; }).

Проблемные гонки:

A. Terminate vs Exchange (несовпадающие теги).
Процесс A заходит в «завершение» и вызывает sendrecv(tag=2). Процесс B ещё «работает» и в начале своей итерации попадает в ветку обмена (например, schetchick % T == 0 или iprobe(partner,0) ⇒ true). Он вызывает sendrecv(tag=0). Теперь:

  • A ждёт тэг 2, B ждёт тэг 0. Оба блокируются навсегда, потому что ни один не сможет «помочь» другому: оба уже сидят внутри блокирующих sendrecv с разными тегами. Никакие iprobe уже не сработают — их никто не выполняет (все в блокирующих вызовах).

B. Очередь хранит «старый» обмен (tag=0), когда партнёр уже ушёл в завершение (tag=2).
Даже если B не попал по schetchick % T, но в очереди остался прежний пакет tag=0 (или только что пришёл eager-пакет), iprobe(partner,0) сработает раньше, чем вы проверите iprobe(partner,2) (поскольку в коде проверка на обмен стоит раньше проверки на завершение). B снова уйдёт в sendrecv(tag=0), а A уже в sendrecv(tag=2) — тот же взаимный клин.

Итог: текущий порядок проверок и разделение по тегам 0/2 делают дедлок детерминированно возможным.

2) Как это починить (минимально, без полной переписки)

Есть несколько рабочих вариантов. Выберите один и сделайте последовательно — хватит одного.

Вариант 1 — Один тэг для всего (+ тип сообщения внутри полезной нагрузки)

Самый надёжный и простой для двух процессов.

  • Вводим общий тэг, например TAG_CTRL = 0 для любого сообщения.
  • Полезная нагрузка содержит поле kind:
    • kind = EXCHANGE (+ поля CrossMsg),
    • kind = TERMINATE (дальнейшие поля можно не использовать / держать заглушку фиксированного размера).
  • Условие входа в обмен/завершение сводим к одному sendrecv(TAG_CTRL):
    • Вызываем sendrecv(TAG_CTRL) если (a) пора обмениваться по локальному правилу, или (b) локально принято решение завершаться, или (c) iprobe(partner, TAG_CTRL) истинно.
    • Кто хочет завершаться — посылает kind=TERMINATE. Кто хочет обмениваться — посылает kind=EXCHANGE. Сообщения по одному тегу всегда совпадут.
  • На приёме смотрим kind партнёра:
    • Если он прислал TERMINATE, вы немедленно переходите к финализации (и опционально шлёте TERMINATE в ответ — но вы и так сделали sendrecv, так что обе стороны уже синхронно встретились и разошлись).

Это полностью убирает «перекрёст тэгов 0/2». А «раньше/позже» сводится к тому, что все важные синхронизации идут через один и тот же вызов sendrecv(TAG_CTRL).

Вариант 2 — если тэги оставлять разными (быстрый патч)

  1. Приоритет завершения. В начале итерации сначала проверяйте iprobe(partner, 2) и, если он истинный — сразу входите в путь завершения (без захода в обмен).
    То есть проверку tag=2 поднимите выше и отделите от вычислений.

  2. Перед тем как звать sendrecv(tag=2), аккуратно «осушите» возможные незакрытые обмены tag=0, чтобы не оставить партнёра висеть в sendrecv(tag=0). Например:

    cpp
    while (g_world->iprobe(partner, 0)) {
    CrossMsg junk;
    g_world->recv(partner, 0, junk);      // закрывает его send-часть
    g_world->send(partner, 0, junk);      // и его recv-часть sendrecv'а
}

    После этого вызывайте sendrecv(partner,2,...).
    Это избавляет от варианта «A ждёт tag=2, B — tag=0».

  3. В обменном блоке перед sendrecv(tag=0) на всякий случай проверьте iprobe(partner, 2) и, если он истинен — не заходите в обмен, а переходите к завершению.

Минус варианта 2: это всё равно менее надёжно, чем единый тег, и сложнее поддерживать.

Вариант 3 — всегда держать по одному Irecv на каждый тэг

  • С самого начала цикла постите постоянные MPI_Irecv на tag 0 и tag 2 (каждый — на один буфер).
  • Внутри итерации проверяйте MPI_Test соответствующего реквеста и реагируйте.
  • Для отправки используйте MPI_Isend и MPI_Wait* там, где требуется упорядочивание.
    Это тоже снимает необходимость гадать, будет ли iprobe что-то видеть.

3) Очень вероятная не MPI причина «иногда работает, иногда зависает/падает»

Ваш Slab/TLS-аллокатор никогда не «сбрасывается»

cpp
static tbb::enumerable_thread_specific<Slab*> tls([]{
    void* memory = _aligned_malloc(16777216u, 16u);
    Slab* slab = static_cast<Slab*>(_aligned_malloc(32u, 16u));
    new (slab) Slab(memory, 16777216u);

    // ... прогрев страниц ...
    return slab;
}());
  • Interval::operator new просто делает s->current += 64; — никакой проверки на выход за пределы slab->end.
  • Между вызовами AGP_2D вы никогда не делаете tls.local()->current = tls.local()->base; (или другой reset).
  • После нескольких запусков (вы сами пишете «несколько нажатий обрабатываются — на некотором нажатии дедлок») — вы уходите за пределы 16МБ слайда, начинаете переписывать произвольную память. Дальше — уже чистый UB: от «молчаливого» вырождения heap/падений до «зависания» в MPI из-за испорченных структур.

Обязательно добавьте сброс арены в начале AGP_1D/AGP_2D, например:

cpp
auto* s = tls.local();
s->current = s->base;

И, на время отладки, добавьте ассерт:

cpp
if (s->current + 64 > s->end) { /* лог/асерт/выделить новый слайс */ }

Это, на практике, зачастую «чинит» загадочные зависания.

4) Ещё проблемы/рекомендации по коду

4.1 Финализация MPI

Вы создаёте boost::mpi::environment и communicator через new и никогда не удаляете:

cpp
static const boost::mpi::environment* g_env;
static const boost::mpi::communicator* g_world;

g_env  = new boost::mpi::environment();
g_world= new boost::mpi::communicator();

Из-за этого MPI_Finalize (который вызывается деструктором environment) не вызывается, и mpiexec закономерно ругается: process exited without calling finalize.
Используйте обычные объекты с автоматическим временем жизни (RAII), либо в конце DLL явно делайте delete g_world; delete g_env;. В DLL чаще используют статические объекты на модуль.

4.2 Порядок проверок в главном цикле AGP_2D

Даже если останетесь на 0/2 тегах — поменяйте порядок:

  1. Сначала проверяйте «пришёл terminate» (iprobe(tag=2)) → если да, идите в завершение.
  2. Потом решайте вопрос об обмене (%T или iprobe(tag=0)).

Так вы хотя бы уменьшаете «окно» дедлока.

4.3 Вывод результата «только у одного ранга»

Вы делаете:

cpp
if (partner) { /* записать Extr в out_data */ }

При rank=0 → partner=1 (истина) → пишем.
При rank=1 → partner=0 (ложь) → не пишем.
Это работает, если вы действительно вызываете DLL-функцию в обоих процессах, но потребитель результата — только у rank=0. Всё равно лучше пишите явно if (rank == 0) — так понятнее и безопаснее (сейчас завязка на «парности» рангов и числе процессов ровно 2).

4.4 FABE13_SINCOS и векторный путь

  • Алгоритмически и по коэффициентам полиномы согласованы со скалярной веткой; знаки и отражения по квадрантам сделаны корректно.
  • Но вы используете _mm256_load_ps/_store_ps (требуют 32-байтного выравнивания). В ваших текущих вызовах массивы размечены __declspec(align(32)), так что ОК. Если макрос вызовут когда-нибудь на невыровненном буфере — будет крэш. Если хотите сделать безопаснее, поменяйте на loadu/storeu или задокументируйте жёсткое требование выравнивания.
  • Трюк с «магическим bias» (12582912.0f) для округления векторного q корректен для диапазона, где |x| * 1/(2π) < 2^24. У вас аргументы небольшие — нормально.

4.5 ldexpf(1.0f, -(1.0f / 138.63f) * ++schetchick);

Второй аргумент ldexpf — целый. Сейчас у вас плавающее выражение неявно приводится к int, что даёт ступенчатое поведение (0,0,0,...,-1,-1, ...). Судя по названию «cooling», вы хотели плавную экспоненту. Тогда используйте exp2f(-(1.0f/138.63f) * ++schetchick) или powf(2.0f, ...).

4.6 RNG-макросы и типы

Вы вызываете XOR_RAND(current_state, ...), где current_state — float. Макрос же рассчитывает на int state (int s = state; ... state = s;). Это:

  • постоянно делает float → int → float преобразования (потеря старших бит при state > 2^24),
  • потенциально ломает период/статистику генератора.
    Храните состояние как uint32_t и, если нужно, приводите результат отдельно.

4.7 «const-ов больше, чем нужно»

  • Archive& const __restrict ar в serialize — бессмысленно: const на ссылке ничего не даёт. Оставьте Archive& ar.
  • const enum List : uint8_t ... — const здесь тоже ни о чём.
    Это не «ломает» код, но вредит читабельности и может мешать IDE/анализаторам.

4.8 Микро-замечания

  • Проверьте возможность дубликатов указателей в R (в логике не должно быть, но если вдруг попадут) — ваши функции с __restrict могут агрессивно оптимизироваться при предположении «указатели не alias-ятся».
  • В AgpWaitStartAndRun — бесконечный busy-loop. Добавьте std::this_thread::sleep_for(…) хотя бы на миллисекунды, чтобы не сжигать ядро впустую.
  • В 2D-ветке переменная top_ptr в начале/конце итерации переиспользуется — это ок, просто следите, чтобы не использовать её после изменений кучи (push_heap/emplace_back) до переобновления.

5) Мини-рефактор (пример) с единым тегом

cpp
enum Kind : int { EXCHANGE=0, TERMINATE=1 };

struct CtrlMsg {
  Kind kind;
  CrossMsg xchg; // валиден только для EXCHANGE, можно фикс. размер
};

// Условие входа в обмен/завершение:
bool want_xchg = (schetchick % T == 0);
bool want_term = (interval_len < epsilon) || (schetchick == (int)global_iterations);

if ( want_xchg || want_term || g_world->iprobe(partner, TAG_CTRL) ) {
    CtrlMsg out{}, in{};
    out.kind = want_term ? TERMINATE : EXCHANGE;
    if (out.kind == EXCHANGE) {
       // заполнить out.xchg как у вас
    }
    g_world->sendrecv(partner, TAG_CTRL, out, partner, TAG_CTRL, in);

    if (in.kind == TERMINATE) {
        // аккуратно выходим, отправив при необходимости финальные данные
        // (можно положить их в следующий CtrlMsg или отдельным sendrecv)
        finalize_and_return();
    }
    if (out.kind == TERMINATE) {
        // мы просили завершиться; партнёр ещё мог прислать EXCHANGE — просто игнорируем
        finalize_and_return();
    }
    // иначе in.kind == EXCHANGE: применяем пришедший интервал
}

Этот рисунок исключает все «перекрёстные тэги», снимает зависимость от eager-буферизации и делает поведение воспроизводимым.

6) Резюме actionable

  1. Исправьте протокол: или единый тег + поле kind, или приоритет terminate + «осушка» tag=0 перед sendrecv(tag=2).
  2. В начале каждого запуска сбрасывайте TLS-слайс: tls.local()->current = tls.local()->base; + защита от переполнения.
  3. Приведите к порядку финализацию MPI (RAII/delete).
  4. Поправьте ldexpf → exp2f/powf.
  5. Уберите хранение состояния RNG в float.
  6. Проверьте выравнивание буферов для FABE13_SINCOS или используйте loadu/storeu.

Если сделать хотя бы (1) и (2), шансы увидеть «дедлок на N-ом нажатии» резко уменьшаются: вы убираете систематическую гонку тегов и убираете скрытую порчу памяти, которая могла маскироваться под «MPI завис».

Compartir este Q&A