пишу крайне высокопроизводительное приложение став...

Question

пишу крайне высокопроизводительное приложение ставлю рекорд гинеса по скорости АГП вот мой код из dll библиотеки который будет использоваться в управляемом коде c++/clr дай конкретные советы по увеличению скорости в частности по классу который нужен для предварительного выделения памяти и выравнивания - читаемость и надёжность абсолютно не важны - только максимальная скорость #include "pch.h"
#define LCG_RAND(state)
do {
float reduced = fmodf((state), 1.618033988749895f);
reduced = fmaf(reduced, 1.618033988749895f, reduced * 0.5772156649f);
reduced -= (float)(int)(reduced);
reduced = fmaf(reduced, reduced + 1.0f, 0.0f);
reduced -= (float)(int)(reduced);
(state) = reduced;
} while(0)

#define LCG_RAND_GRSH(state)
do {
float reduced = fmodf((state), 1.618033988749895f);
reduced = fmaf(reduced, 1.618033988749895f, reduced * 0.5772156649f);
reduced -= (float)(int)(reduced);
reduced = fmaf(reduced, reduced + 1.0f, 0.0f);
reduced -= (float)(int)(reduced);
(state) = fmaf(reduced, 2.0f, -1.0f);
} while(0)

#define FABE13_COS(x, result_var)
do {
float abs_val = fabsf((x));
float reduced = fmodf(abs_val, 6.28318530718f);
if (reduced > 3.14159265359f) {
reduced = 6.28318530718f -reduced;
}
if (reduced < 1.57079632679f) {
float val2 = reduced * reduced;
result_var = fmaf(val2, fmaf(val2, 0.0416666667f, -0.5f), 1.0f);
} else {
reduced = 3.14159265359f -reduced;
float val2 = reduced * reduced;
result_var = -fmaf(val2, fmaf(val2, 0.0416666667f, -0.5f), 1.0f);
}
} while(0)

#define FABE13_SINCOS(in, sin_out, cos_out, n)
do {
int i = 0;
int limit = (n) & ~7;
if ((n) >= 8) {
const __m256 VEC_TWOPI = _mm256_set1_ps(6.28318530718f);
const __m256 VEC_PI = _mm256_set1_ps(3.14159265359f);
const __m256 VEC_PI_2 = _mm256_set1_ps(1.57079632679f);
const __m256 INV_TWOPI = _mm256_set1_ps(0.15915494309189535f);
const __m256 BIAS = _mm256_set1_ps(12582912.0f);
const __m256 VEC_COS_P3 = _mm256_set1_ps(0.0416666667f);
const __m256 VEC_COS_P1 = _mm256_set1_ps(-0.5f);
const __m256 VEC_COS_P0 = _mm256_set1_ps(1.0f);
const __m256 VEC_SIN_P1 = _mm256_set1_ps(-0.16666666f);
const __m256 VEC_SIN_P0 = _mm256_set1_ps(1.0f);
const __m256 VEC_ZERO = _mm256_setzero_ps();
while (i < limit) {
__m256 vx = _mm256_loadu_ps(&(in)[i]);
__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);
__m256 r = _mm256_fnmadd_ps(VEC_TWOPI, q, vax);
__m256 r1 = _mm256_min_ps(r, _mm256_sub_ps(VEC_TWOPI, r));
__m256 r2 = _mm256_min_ps(r1, _mm256_sub_ps(VEC_PI, r1));
__m256 t2 = _mm256_mul_ps(r2, r2);
__m256 cosv = _mm256_fmadd_ps(t2, _mm256_fmadd_ps(t2, VEC_COS_P3, VEC_COS_P1), VEC_COS_P0);
__m256 sinv = _mm256_mul_ps(_mm256_fmadd_ps(t2, VEC_SIN_P1, VEC_SIN_P0), r2);
__m256 cflip = _mm256_cmp_ps(r1, VEC_PI_2, _CMP_GT_OQ);
__m256 sflip = _mm256_xor_ps(_mm256_cmp_ps(vx, VEC_ZERO, _CMP_LT_OQ),
_mm256_cmp_ps(r, VEC_PI, _CMP_GT_OQ));
cosv = _mm256_blendv_ps(cosv, _mm256_sub_ps(VEC_ZERO, cosv), cflip);
sinv = _mm256_blendv_ps(sinv, _mm256_sub_ps(VEC_ZERO, sinv), sflip);
_mm256_storeu_ps(&(cos_out)[i], cosv);
_mm256_storeu_ps(&(sin_out)[i], sinv);
i += 8;
}
}
while (i < (n)) {
float x = (in)[i];
float ax = fabsf(x);
float q = fmaf(ax, 0.15915494309189535f, 12582912.0f);
q -= 12582912.0f;
float r = fmaf(-6.28318530718f, q, ax);
bool sflip = (r > 3.14159265359f);
if (sflip) r = 6.28318530718f - r;
bool cflip = (r > 1.57079632679f);
if (cflip) r = 3.14159265359f - r;
float t2 = r * r;
float c = fmaf(t2, fmaf(t2, 0.0416666667f, -0.5f), 1.0f);
float s = fmaf(t2, -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, res = 0.0f;
while (i < 10) {
LCG_RAND(current_state);
const float x_part = fmaf(-current_state, 10.0f, x);
LCG_RAND(current_state);
current_res = current_state;
LCG_RAND(current_state);
float delimiter = fmaf(fmaf(current_res, 20.0f, 5.0f),
x_part * x_part,
fmaf(current_state, 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(16)) float angles[14u];
const float start_angle = 6.28318530717958647692f * x;
#pragma loop(ivdep)
while (j < 14) {
angles[j] = start_angle * (j + 1);
++j;
}
__declspec(align(16)) float sin_vals[14u];
__declspec(align(16)) float cos_vals[14u];
FABE13_SINCOS(angles, sin_vals, cos_vals, 14u);
float current_state = seed;
LCG_RAND(current_state);
float res = fmaf(current_state, 2.0f, -1.1f);
--j;
while (j >= 0) {
LCG_RAND(current_state);
float tmp_state = current_state;
LCG_RAND(current_state);
res += fmaf(fmaf(tmp_state, 2.0f, -1.1f), sin_vals[j],
fmaf(current_state, 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(16)) float angles_j[8u];
__declspec(align(16)) float angles_k[8u];
#pragma loop(ivdep)
while (j < 8) {
const float pj_mult = 3.14159265358979323846f * (j + 1);
angles_j[j] = pj_mult * x1;
angles_k[j] = pj_mult * x2;
++j;
}
__declspec(align(16)) float sin_j[8u], cos_j[8u];
__declspec(align(16)) 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;
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];
LCG_RAND_GRSH(current_state);
float tmp_state = current_state;
LCG_RAND_GRSH(current_state);
part1 = fmaf(tmp_state, sin_term,
fmaf(current_state, cos_term, part1));
LCG_RAND_GRSH(tmp_state);
LCG_RAND_GRSH(current_state);
part2 = fmaf(-tmp_state, cos_term,
fmaf(current_state, 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(
0x3F800000 | ((((const int)&diff) & 0x80000000) ^ 0x80000000)
)));
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 Interval final {
const float x1;
const float x2;
const float y1;
const float y2;
const float delta_y;
const float N_factor;
const float M;
float R;

text
__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),
	N_factor(_N == 1.0f ? _x2 - _x1 :
		_N == 2.0f ? sqrtf(_x2 - _x1) :
		powf(_x2 - _x1, 1.0f / _N)),
	M(fabsf(delta_y)* (1.0f / N_factor)) {
}

__declspec(noalias) __forceinline void ChangeCharacteristic(const float _m) noexcept
{
	R = fmaf(
		-(y2 + y1),
		2.0f,
		fmaf(_m, N_factor, (delta_y * delta_y) * (1.0f / (_m * N_factor))));
}

};

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

__declspec(align(16)) struct PreallocatedPeanoMemoryManager final {
__declspec(align(16)) struct MemoryBlock final {
char* __restrict data;
char* __restrict current;
const char* __restrict end;

text
	__declspec(noalias) __forceinline void initialize(const size_t required_size) noexcept {
		data = static_cast<char*>(_aligned_malloc(required_size, 16u));
		current = data;
		end = data + required_size;
	}

	__declspec(noalias) __forceinline void* allocate_from_block(const size_t alloc_size) noexcept {
		const size_t aligned_size = (alloc_size + 15u) & ~(15u);
		char* result = current;
		current += aligned_size;
		return result;
	}
};

tbb::enumerable_thread_specific<MemoryBlock> thread_local_blocks;

__declspec(noalias) __forceinline size_t calculate_tree_size(const int max_depth) const noexcept {
	static const size_t precomputed_sizes[] = {
			1398080u,
			5592384u,
			22369600u,
			89478464u,
			357913920u,
			1431655744u,
			5726623040u
	};
	return precomputed_sizes[max_depth - 7];
}

__declspec(noalias) __forceinline void preallocate_for_depth(const int max_depth) noexcept {
	const size_t required_size = calculate_tree_size(max_depth);
	thread_local_blocks.clear();
	thread_local_blocks = tbb::enumerable_thread_specific<MemoryBlock>(
		[required_size]() {
			MemoryBlock block;
			block.initialize(required_size);
			return block;
		}
	);
}

__declspec(noalias) __forceinline void* allocate(const size_t size) noexcept {
	return thread_local_blocks.local().allocate_from_block(size);
}

};

static PreallocatedPeanoMemoryManager peano_memory_manager;

const enum List : uint8_t {
Top = 0b00,
Down = 0b01,
Left = 0b10,
Right = 0b11
};

__declspec(align(16)) struct PeanoCurve_2D final {
const List Type;
const int razvertka;
const float a;
const float b;
const float c;
const float d;
const float x1;
const float x2;
PeanoCurve_2D* __restrict DownLeft;
PeanoCurve_2D* __restrict TopLeft;
PeanoCurve_2D* __restrict TopRight;
PeanoCurve_2D* __restrict DownRight;

text
__declspec(noalias) __forceinline void* operator new(size_t size) {
	return peano_memory_manager.allocate(size);
}

__declspec(noalias) __forceinline PeanoCurve_2D(
	const List _Type,
	int _razvertka,
	const float _a,
	const float _b,
	const float _c,
	const float _d) noexcept
	: Type(_Type),
	razvertka(_razvertka),
	a(_a),
	b(_b),
	c(_c),
	d(_d),
	x1(fmaf(0.5f, _a, 0.5f * _b)),
	x2(fmaf(0.5f, _c, 0.5f * _d)),
	DownLeft(nullptr),
	TopLeft(nullptr),
	TopRight(nullptr),
	DownRight(nullptr) {

	if (_razvertka-- != 0) {
		if (_razvertka > 8) {
			switch (_Type) {
			case Top:
				tbb::parallel_invoke(
					[&]() { DownLeft = new PeanoCurve_2D(Right, _razvertka, _a, x1, _c, x2); },
					[&]() { TopLeft = new PeanoCurve_2D(Top, _razvertka, _a, x1, x2, _d); },
					[&]() { TopRight = new PeanoCurve_2D(Top, _razvertka, x1, _b, x2, _d); },
					[&]() { DownRight = new PeanoCurve_2D(Left, _razvertka, x1, _b, _c, x2); }
				);
				break;
			case Down:
				tbb::parallel_invoke(
					[&]() { TopRight = new PeanoCurve_2D(Left, _razvertka, x1, _b, x2, _d); },
					[&]() { DownRight = new PeanoCurve_2D(Down, _razvertka, x1, _b, _c, x2); },
					[&]() { DownLeft = new PeanoCurve_2D(Down, _razvertka, _a, x1, _c, x2); },
					[&]() { TopLeft = new PeanoCurve_2D(Right, _razvertka, _a, x1, x2, _d); }
				);
				break;
			case Right:
				tbb::parallel_invoke(
					[&]() { DownLeft = new PeanoCurve_2D(Top, _razvertka, _a, x1, _c, x2); },
					[&]() { DownRight = new PeanoCurve_2D(Right, _razvertka, x1, _b, _c, x2); },
					[&]() { TopRight = new PeanoCurve_2D(Right, _razvertka, x1, _b, x2, _d); },
					[&]() { TopLeft = new PeanoCurve_2D(Down, _razvertka, _a, x1, x2, _d); }
				);
				break;
			case Left:
				tbb::parallel_invoke(
					[&]() { TopRight = new PeanoCurve_2D(Down, _razvertka, x1, _b, x2, _d); },
					[&]() { TopLeft = new PeanoCurve_2D(Left, _razvertka, _a, x1, x2, _d); },
					[&]() { DownLeft = new PeanoCurve_2D(Left, _razvertka, _a, x1, _c, x2); },
					[&]() { DownRight = new PeanoCurve_2D(Top, _razvertka, x1, _b, _c, x2); }
				);
			}
		}
		else {
			switch (_Type) {
			case Top:
				DownLeft = new PeanoCurve_2D(Right, _razvertka, _a, x1, _c, x2);
				TopLeft = new PeanoCurve_2D(Top, _razvertka, _a, x1, x2, _d);
				TopRight = new PeanoCurve_2D(Top, _razvertka, x1, _b, x2, _d);
				DownRight = new PeanoCurve_2D(Left, _razvertka, x1, _b, _c, x2);
				break;
			case Down:
				TopRight = new PeanoCurve_2D(Left, _razvertka, x1, _b, x2, _d);
				DownRight = new PeanoCurve_2D(Down, _razvertka, x1, _b, _c, x2);
				DownLeft = new PeanoCurve_2D(Down, _razvertka, _a, x1, _c, x2);
				TopLeft = new PeanoCurve_2D(Right, _razvertka, _a, x1, x2, _d);
				break;
			case Right:
				DownLeft = new PeanoCurve_2D(Top, _razvertka, _a, x1, _c, x2);
				DownRight = new PeanoCurve_2D(Right, _razvertka, x1, _b, _c, x2);
				TopRight = new PeanoCurve_2D(Right, _razvertka, x1, _b, x2, _d);
				TopLeft = new PeanoCurve_2D(Down, _razvertka, _a, x1, x2, _d);
				break;
			case Left:
				TopRight = new PeanoCurve_2D(Down, _razvertka, x1, _b, x2, _d);
				TopLeft = new PeanoCurve_2D(Left, _razvertka, _a, x1, x2, _d);
				DownLeft = new PeanoCurve_2D(Left, _razvertka, _a, x1, _c, x2);
				DownRight = new PeanoCurve_2D(Top, _razvertka, x1, _b, _c, x2);
			}
		}
	}
}

const __declspec(noalias) __forceinline PeanoCurve_2D* __restrict HitTest_2D(float x) const noexcept
{
	int i = 0;
	const int _razvertka = this->razvertka;
	int num;
	const float this_a = this->a;
	x -= this_a;
	const float b_minus_a = this->b - this_a;
	const float inv_b_minus_a = 1.0f / b_minus_a;
	const PeanoCurve_2D* __restrict Curr = this;
	while (i != _razvertka) {
		const int shift = 1 << ++i + i;
		num = shift * x * inv_b_minus_a;
		x = fmaf(-ldexp(1.0f, -(i << 1)) * num, b_minus_a, x);
		const List currType = Curr->Type;
		switch (num) {
		case 0:
			Curr = (currType == Top || currType == Right)
				? Curr->DownLeft
				: Curr->TopRight;
			break;
		case 1:
			Curr = (currType == Top || currType == Left)
				? Curr->TopLeft
				: Curr->DownRight;
			break;
		case 2:
			Curr = (currType == Top || currType == Right)
				? Curr->TopRight
				: Curr->DownLeft;
			break;
		default:
			Curr = (currType == Top || currType == Left)
				? Curr->DownRight
				: Curr->TopLeft;
		}
	}
	return const_cast<PeanoCurve_2D * __restrict>(Curr);
}

const __declspec(noalias) __forceinline float FindX_2D(const float target_x1, const float target_x2) const noexcept
{
	int _razvertka = this->razvertka;
	int _razvertka1 = _razvertka;
	float x1, x2, x = this->a;
	const float b_minus_a = this->b - x;
	const PeanoCurve_2D* __restrict Curr = this;
	while (_razvertka != 0) {
		const int exponent = _razvertka1 - _razvertka-- << 1;
		x1 = Curr->x1;
		x2 = Curr->x2;
		const List currType = Curr->Type;
		if (target_x1 > x1 && target_x2 > x2) {
			Curr = Curr->TopRight;
			if (currType == Top || currType == Right) {
				x = fmaf(ldexpf(1.0f, -exponent) * 0.5f, b_minus_a, x);
			}
		}
		else if (target_x1 < x1 && target_x2 > x2) {
			Curr = Curr->TopLeft;
			if (currType == Top || currType == Left) {
				x = fmaf(ldexpf(1.0f, -exponent) * 0.25f, b_minus_a, x);
			}
			else {
				x = fmaf(ldexpf(1.0f, -exponent) * 0.75f, b_minus_a, x);
			}
		}
		else if (target_x1 < x1 && target_x2 < x2) {
			Curr = Curr->DownLeft;
			if (currType == Down || currType == Left) {
				x = fmaf(ldexpf(1.0f, -exponent) * 0.5f, b_minus_a, x);
			}
		}
		else {
			Curr = Curr->DownRight;
			if (currType == Top || currType == Left) {
				x = fmaf(ldexpf(1.0f, -exponent) * 0.75f, b_minus_a, x);
			}
			else {
				x = fmaf(ldexpf(1.0f, -exponent) * 0.25f, b_minus_a, x);
			}
		}
	}
	return x;
}

};

static thread_local std::vector<float, tbb::scalable_allocator<float>> Extr;
static thread_local std::vector<Interval*, tbb::scalable_allocator<Interval*>> R;
static boost::mpi::environment* g_env;
static boost::mpi::communicator* g_world;
static std::aligned_storage_t<sizeof(PeanoCurve_2D), alignof(PeanoCurve_2D)> curveStorage;
static PeanoCurve_2D& Curve = reinterpret_cast<PeanoCurve_2D&>(curveStorage);
static boost::optionalboost::mpi::status probe_status;

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 {
peano_memory_manager.preallocate_for_depth(peanoLevel);
g_env = ::new boost::mpi::environment();
g_world = ::new boost::mpi::communicator();
int rank = g_world->rank();
if (rank) {
::new (&Curve) PeanoCurve_2D(List::Down, peanoLevel, a, b, c, d);
}
else {
::new (&Curve) PeanoCurve_2D(List::Top, peanoLevel, a, b, c, d);
}
return rank;
}

extern "C" __declspec(dllexport) __declspec(noalias) __forceinline void AgpFinalize() noexcept {
delete g_world;
delete g_env;
}

__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& __restrict const ar, const int) noexcept { ar& s_x1& s_x2& e_x1& e_x2& Rtop; }
};

extern "C" __declspec(dllexport) __declspec(noalias)
void Base_LNA_1_2_Mer_AGP(
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** __restrict out_data, size_t* __restrict out_len) noexcept
{
const int rank = g_world->rank();
const int partner = rank ^ 1;
float divider, inv_divider, x_addition, y_addition, true_start, true_end;
float x_Rmax_1, x_Rmax_2;
float initial_length, dmax, threshold_03, inv_threshold_03;
float start_val, best_f, y_Rmax_1, y_Rmax_2;
float Mmax, m;
int schetchick = 0, mcQueenSpeed = 1;
float new_point, new_value;
const PeanoCurve_2D* __restrict pc;
float new_x1, new_x2;
Interval* __restrict promejutochny_otrezok;
Interval* __restrict curr;
Interval* __restrict curr1;
float currM, len1, len2, len_item;
size_t r_size;
float progress, alpha, betta, MULTIPLIER, global_coeff, GLOBAL_FACTOR;
Interval* __restrict top_ptr;
float interval_len;
int dummy;
const PeanoCurve_2D* __restrict p1L;
const PeanoCurve_2D* __restrict p2L;
CrossMsg outbound, inbound{};
Interval* __restrict injected;
float newParam1_remote1, newParam2_remote1, newParam1_remote2, newParam2_remote2;
float cooling;
int T;
float k;
int tag;

text
if (N == 1.0f)
{
	initial_length = b - a;
	dmax = initial_length;
	threshold_03 = 0.3f * initial_length;
	inv_threshold_03 = 1.0f / threshold_03;
	start_val = ShekelFunc(a, seed);
	best_f = ShekelFunc(b, seed);
	x_Rmax_1 = a;
	x_Rmax_2 = b;
	y_Rmax_1 = start_val;
	y_Rmax_2 = best_f;
	Extr.reserve(static_cast<size_t>(global_iterations) * 4u);
	Extr.clear();
	R.reserve(static_cast<size_t>(global_iterations) * 2u);
	R.clear();
	R.emplace_back(new Interval(a, b, start_val, best_f, N));
	std::push_heap(R.begin(), R.end(), ComparePtr);
	Mmax = R.front()->M;
	m = r * Mmax;
	float x1, y1, x2, y2;
	while (true)
	{
		new_point = Shag(m, x_Rmax_1, x_Rmax_2, y_Rmax_1, y_Rmax_2, N, r);
		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);
		promejutochny_otrezok = R.back();
		R.pop_back();
		x1 = promejutochny_otrezok->x1;
		y1 = promejutochny_otrezok->y1;
		x2 = promejutochny_otrezok->x2;
		y2 = promejutochny_otrezok->y2;
		len2 = x2 - new_point;
		len1 = new_point - x1;
		interval_len = (len1 < len2) ? len1 : len2;
		if (++schetchick == static_cast<int>(global_iterations) || interval_len < epsilon)
		{
			if (rank == 0) {
				Extr.emplace_back(static_cast<float>(schetchick));
				Extr.push_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;
		}
		curr = new Interval(x1, new_point, y1, new_value, N);
		curr1 = new Interval(new_point, x2, new_value, y2, N);
		currM = (curr->M > curr1->M) ? curr->M : curr1->M;
		r_size = R.size();
		if ((len2 + len1) == dmax) {
			dmax = (len2 > len1) ? len2 : len1;
			size_t i = 0u;

#pragma loop(ivdep)
while (i < r_size) {
len_item = R[i]->x2 - R[i]->x1;
if (len_item > dmax) dmax = len_item;
++i;
}
}
if (mode) {
if ((threshold_03 > dmax && fmodf(static_cast<float>(schetchick), 3.0f) == 0.0f) || 10.0f * dmax < initial_length) {
if (currM > Mmax) { Mmax = currM; m = r * Mmax; }
progress = fmaf(-inv_threshold_03, dmax, 1.0f);
alpha = fmaf(progress, progress, 1.0f);
betta = 2.0f - alpha;
MULTIPLIER = (1.0f / dmax) * Mmax;
global_coeff = fmaf(MULTIPLIER, r, -MULTIPLIER);
GLOBAL_FACTOR = betta * global_coeff;
curr->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, len1, curr->M * alpha));
curr1->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, len2, curr1->M * alpha));
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;
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
else {
if (currM > Mmax) {
Mmax = currM; m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) { R[i]->ChangeCharacteristic(m); ++i; }
std::make_heap(R.begin(), R.end(), ComparePtr);
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}
}
else {
if (currM > Mmax) {
Mmax = currM; m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) { R[i]->ChangeCharacteristic(m); ++i; }
std::make_heap(R.begin(), R.end(), ComparePtr);
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}
R.emplace_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;
}
}
else
{
divider = ldexpf(1.0f, (Curve.razvertka << 1) + 1);
inv_divider = 1.0f / divider;
x_addition = (b - a) * inv_divider;
y_addition = (d - c) * inv_divider;
true_start = a + x_addition;
true_end = b - x_addition;
x_Rmax_1 = true_start;
x_Rmax_2 = true_end;
initial_length = true_end - true_start;
dmax = initial_length;
threshold_03 = 0.3f * initial_length;
inv_threshold_03 = 1.0f / threshold_03;
start_val = rank ? RastriginFunc(true_end, d - y_addition)
: RastriginFunc(true_start, c + y_addition);
y_Rmax_1 = start_val;
best_f = rank ? RastriginFunc(true_start, d - y_addition)
: RastriginFunc(true_end, c + y_addition);
y_Rmax_2 = best_f;
Extr.reserve(static_cast<size_t>(global_iterations) * 4u);
Extr.clear();
R.reserve(static_cast<size_t>(global_iterations) * 2u);
R.clear();
R.emplace_back(new Interval(true_start, true_end, start_val, best_f, N));
Mmax = R.front()->M;
m = r * Mmax;
while (true)
{
if (++schetchick == static_cast<int>(global_iterations)) {
g_world->isend(partner, 2, 0);
if (partner) {
top_ptr = R.front();
Extr.emplace_back(static_cast<float>(schetchick));
Extr.push_back(top_ptr->x2 - top_ptr->x1);
*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;
}
cooling = ldexpf(1.0f, -(1.0f / 138.63f) * schetchick);
T = static_cast<int>(fmaf(20.0f, cooling, 10.0f));
k = fmaf(0.2f, cooling, 0.7f);
if (schetchick % (T - 1) == 0) {
probe_status = g_world->iprobe(partner, boost::mpi::any_tag);
if (probe_status) {
tag = probe_status->tag();
if (tag == 0) {
g_world->irecv(partner, 0, mcQueenSpeed);
}
else if (tag == 2) {
g_world->irecv(partner, 2, dummy);
if (partner) {
top_ptr = R.front();
Extr.emplace_back(static_cast<float>(schetchick));
Extr.push_back(top_ptr->x2 - top_ptr->x1);
*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;
}
}
}
if (!(schetchick % T) == mcQueenSpeed) {
top_ptr = R.front();
p1L = Curve.HitTest_2D(top_ptr->x1);
p2L = Curve.HitTest_2D(top_ptr->x2);
outbound = CrossMsg{ p1L->x1, p1L->x2, p2L->x1, p2L->x2, top_ptr->R };
if (mcQueenSpeed) {
g_world->isend(partner, 0, 0);
}
else {
mcQueenSpeed = 1;
}
g_world->isend(partner, 1, outbound);
g_world->irecv(partner, 1, inbound);
newParam1_remote1 = Curve.FindX_2D(inbound.s_x1, inbound.s_x2);
newParam2_remote1 = Curve.FindX_2D(inbound.e_x1, inbound.e_x2);
newParam1_remote2 = RastriginFunc(inbound.s_x1, inbound.s_x2);
newParam2_remote2 = RastriginFunc(inbound.e_x1, inbound.e_x2);
injected = new Interval(newParam1_remote1, newParam2_remote1, newParam1_remote2, newParam2_remote2, N);
injected->R = inbound.Rtop * k;
R.emplace_back(injected);
std::push_heap(R.begin(), R.end(), ComparePtr);
}
new_point = Shag(m, x_Rmax_1, x_Rmax_2, y_Rmax_1, y_Rmax_2, N, r);
pc = Curve.HitTest_2D(new_point);
new_x1 = pc->x1;
new_x2 = pc->x2;
new_value = RastriginFunc(new_x1, new_x2);
if (new_value < best_f) {
best_f = new_value;
Extr.emplace_back(best_f);
Extr.emplace_back(new_x1);
Extr.emplace_back(new_x2);
}
std::pop_heap(R.begin(), R.end(), ComparePtr);
promejutochny_otrezok = R.back();
R.pop_back();
curr = new Interval(promejutochny_otrezok->x1, new_point, promejutochny_otrezok->y1, new_value, N);
curr1 = new Interval(new_point, promejutochny_otrezok->x2, new_value, promejutochny_otrezok->y2, N);
currM = (curr->M > curr1->M) ? curr->M : curr1->M;
len2 = promejutochny_otrezok->x2 - new_point;
len1 = new_point - promejutochny_otrezok->x1;
r_size = R.size();
if ((len2 + len1) == dmax) {
dmax = (len2 > len1) ? len2 : len1;
size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) {
len_item = R[i]->x2 - R[i]->x1;
if (len_item > dmax) dmax = len_item;
++i;
}
}
if (mode) {
if ((threshold_03 > dmax && fmodf(static_cast<float>(schetchick), 3.0f) == 0.0f) || 10.0f * dmax < initial_length) {
if (currM > Mmax) { Mmax = currM; m = r * Mmax; }
progress = fmaf(-inv_threshold_03, dmax, 1.0f);
alpha = fmaf(progress, progress, 1.0f);
betta = 2.0f - alpha;
MULTIPLIER = (1.0f / dmax) * Mmax;
global_coeff = fmaf(MULTIPLIER, r, -MULTIPLIER);
GLOBAL_FACTOR = betta * global_coeff;
curr->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, len1, curr->M * alpha));
curr1->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, len2, curr1->M * alpha));
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;
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
else {
if (currM > Mmax) {
Mmax = currM; m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) { R[i]->ChangeCharacteristic(m); ++i; }
std::make_heap(R.begin(), R.end(), ComparePtr);
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}
}
else {
if (currM > Mmax) {
Mmax = currM; m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) { R[i]->ChangeCharacteristic(m); ++i; }
std::make_heap(R.begin(), R.end(), ComparePtr);
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}
R.emplace_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();
interval_len = top_ptr->x2 - top_ptr->x1;
if (interval_len < epsilon) {
g_world->isend(partner, 2, 0);
if (partner) {
Extr.emplace_back(static_cast<float>(schetchick));
Extr.push_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;
}
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 Base_LNA_1_2_Mer_AGP_Free(float* p) {
CoTaskMemFree(p);
}

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

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

extern "C" __declspec(dllexport) __declspec(noalias) __forceinline void AgpWaitStartAndRun() {
int command = 0;
while (true) {
probe_status = g_world->iprobe(0, 0);
if (probe_status) {
g_world->irecv(0, 0, command);
return;
}
if (command == 1) {
float* __restrict dummy = nullptr;
size_t dummy_len = 0u;
Base_LNA_1_2_Mer_AGP(2.0f, 1000.0f, -2.2f, 1.8f, -2.2f, 1.8f,
2.5f, false, 0.001f, GetTickCount(),
&dummy, &dummy_len);
}
else if (command == 2) {
break;
}
}
}

Answer 1

ниже — только про скорость. читаемость/надёжность игнорирую, в фокусе — память/выравнивание, кэш/TLB, стоимость аллокатора и векторизация.

Что прямо «болит»

Вы выделяете память заранее только под PeanoCurve_2D. Основная же горячая часть плодит кучу Interval через new Interval(...) внутри tight-loop — это медленно (общий аллокатор, синхронизация, фрагментация).
Ваш preallocator кладёт блоки с выравниванием 16 байт, а в коде есть AVX2/FMA. Это лишает вас _mm256_load_ps/_store_ps, cache-line friendliness и увеличивает вероятность split-load’ов.
Нет huge pages/NUMA pinning ⇒ лишние TLB-промахи на больших деревьях.
Векторные буферы на стеке помечены align(16) при том, что вы используете AVX2.
std::push_heap/std::pop_heap + «мелкие new» для интервалов съедают много тактов — часть этого уйдёт, если сделать typed bump-allocator.

Ультра-быстрый bump-arena с правильным выравниванием (дроп-ин)

выравнивание 64 B (cache line, AVX-friendly; если AVX-512 — тоже ок);
попытка huge pages (2 MB) с fallback на _aligned_malloc;
NUMA-локальность «по месту» (выделение в локальном для потока узле);
без проверок границ (вы просили скорость);
один и тот же арена-класс используется и для PeanoCurve_2D, и для Interval.

Вставьте это вместо вашего PreallocatedPeanoMemoryManager и добавьте operator new у Interval.

cpp
// @@ drop-in: быстрая арена с huge pages, NUMA и 64B alignment
#include <windows.h>
#include <processthreadsapi.h>
#include <memoryapi.h>
#include <stdint.h>

__declspec(align(64)) struct PreallocatedPeanoMemoryManager final {
    struct MemoryBlock final {
        char*  data;
        char*  current;
        char*  end;
        size_t mask;   // alignment-1 (для быстрого округления)

        __declspec(noalias) __forceinline void initialize(size_t bytes, size_t alignment = 64, bool tryHuge = true) noexcept {
            mask = alignment - 1u;
            // запас под выравнивание
            size_t need = bytes + alignment;

            // попытка huge pages (если дано право SeLockMemoryPrivilege)
            void* p = nullptr;
            if (tryHuge) {
                SIZE_T hp = GetLargePageMinimum();
                if (hp && (need % hp)) need = (need + hp - 1) & ~(hp - 1);
                if (hp) {
                    p = VirtualAlloc(nullptr, need, MEM_RESERVE | MEM_COMMIT | MEM_LARGE_PAGES, PAGE_READWRITE);
                }
            }
            if (!p) {
                // обычные страницы, но с крупным резервом
                SIZE_T gran = 64 * 1024; // крупная грануляция для TLB locality
                need = (need + gran - 1) & ~(gran - 1);
                p = _aligned_malloc(need, alignment);
            }

            data = reinterpret_cast<char*>(p);
            // вручную довыравниваем старт
            uintptr_t base = reinterpret_cast<uintptr_t>(data);
            uintptr_t aligned = (base + mask) & ~mask;
            current = reinterpret_cast<char*>(aligned);
            end     = reinterpret_cast<char*>(base + need);
        }

        __declspec(noalias) __forceinline void* allocate_from_block(size_t alloc_size) noexcept {
            size_t n = (alloc_size + mask) & ~mask;
            char* p  = current;
            current += n;           // без проверок — MAX speed
            return p;
        }

        __declspec(noalias) __forceinline void release() noexcept {
            // huge pages надо освобождать через VirtualFree, _aligned_malloc — через _aligned_free.
            // Ради скорости не держим ветку; вызываться может один раз на финализации процесса.
        }
    };

    // ETS даёт один блок на поток без CAS/lock при доступе.
    tbb::enumerable_thread_specific<MemoryBlock> tls;

    __declspec(noalias) __forceinline size_t calculate_tree_size_bytes(const int max_depth) const noexcept {
        // Ваши предрасчитанные размеры, но умножим на sizeof(PeanoCurve_2D) / 64 округлённо.
        // Если хотите сохранить точные — оставьте ваш массив. Ниже — ваш original:
        static const size_t bytes[] = {
            1398080u, 5592384u, 22369600u, 89478464u, 357913920u, 1431655744u, 5726623040u
        };
        return bytes[max_depth - 7];
    }

    __declspec(noalias) __forceinline void preallocate_for_depth(const int max_depth,
                                                                 size_t align = 64,
                                                                 bool tryHuge = true) noexcept {
        const size_t required_size = calculate_tree_size_bytes(max_depth);
        tls.clear();
        tls = tbb::enumerable_thread_specific<MemoryBlock>(
            [=]() {
                MemoryBlock b;
                b.initialize(required_size, align, tryHuge);
                return b;
            }
        );
    }

    __declspec(noalias) __forceinline void* allocate(const size_t size) noexcept {
        return tls.local().allocate_from_block(size);
    }
};

// Глобальный арена-менеджер
static PreallocatedPeanoMemoryManager peano_memory_manager;

// Подключаем арену к обоим горячим типам:
__declspec(align(16)) struct Interval final {
    // ... поля как у вас ...
    // быстрый оператор new без delete (arena-подход)
    __declspec(noalias) __forceinline void* operator new(size_t sz) { return peano_memory_manager.allocate(sz); }
    __declspec(noalias) __forceinline void  operator delete(void*) noexcept {} // no-op

    // ваш остальной код Interval без изменений
};

__declspec(align(16)) struct PeanoCurve_2D final {
    // ... ваши поля ...
    __declspec(noalias) __forceinline void* operator new(size_t sz) { return peano_memory_manager.allocate(sz); }
    __declspec(noalias) __forceinline void  operator delete(void*) noexcept {} // no-op
    // ваш остальной код PeanoCurve_2D без изменений
};

Как применять: оставьте ваш вызов peano_memory_manager.preallocate_for_depth(peanoLevel); в AgpInit — теперь он подготовит каждому потоку локальный bump-буфер с 64-байтным выравниванием и (если доступно) на huge pages.

Быстрые точечные правки вокруг выравнивания/векторизации

AVX-friendly стеки. Везде, где у вас массивы для векторной математики, делайте как минимум 32-байтное выравнивание (лучше 64, чтобы совпасть с cache line и арендой):

cpp
alignas(64) float angles[16];
alignas(64) float sin_vals[16];
alignas(64) float cos_vals[16];

(у вас сейчас align(16))

Переключитесь на aligned-loads, когда можете. После пункта (1) входные буферы в FABE13_SINCOS и локальные массивы гарантированно выровнены — можно заменить loadu/storeu на load/store внутри ветки, где вы полностью контролируете источник (ваши стеки/арена):

cpp
// вариант для гарантированно выровненных входов/выходов
#define FABE13_SINCOS_ALIGNED(in, sin_out, cos_out, n) /* как ваш макрос, но _mm256_load_ps/_store_ps */

Подменяйте вызовы на _ALIGNED там, где входные массивы формируете вы сами.

Выравнивание арены → выравнивание полей. Поскольку арена даёт 64-байтный старт и кратные 64-байтам шаги, экземпляры Interval и PeanoCurve_2D тоже будут лежать кратно 64. Это уменьшит ложные sharing-конфликты между потоками и обеспечит best-case для предвыборки.
Non-temporal store для больших стримов. Если Extr растёт крупным (вы его потом копируете в CoTaskMem), при массовой записи промежуточных массивов используйте _mm_stream_ps/_mm256_stream_ps (только для выровненных и без последующего чтения из кэша). Это снимает давление на L3/LLC в «стадионах» итераций.

Убрать «мелкие new» в критической петле

Вы уже резервируете R/Extr, но сами Interval — всё ещё heap allocations. С перегрузкой operator new у Interval (см. код выше) вы превращаете это в простой pointer-bump без синхронизации — обычно минус десятки наносекунд на каждый сплит.

Ещё 10 быстрых выигрышей (низкий риск, высокая отдача)

Компиляция: /O2 /Ot /GL /LTCG /GS- /Gw /Gy /fp:fast /arch:AVX2 (и /Qpar- — TBB сам распараллелит; /Oi для intrinsic).
RNG-макросы: fmodf и приведение к int — дорогие. Если вам не нужна крипто-стойкость, замените на int-LCG/xorshift в uint32_t и масштабируйте к [0,1) умножением на 2^-32. Это срежет десятки тактов в Shekel/Grishagin/Hill.
powf в Shag/Interval: у вас уже есть ветки для N==1/2. Для остальных используйте быстрый exp2f(log2f(x)*p) (с /fp:fast это дешевле powf) или предтабличьте распространённые N.
Heap операций меньше: вы часто делаете pop_heap + два push_heap. Дешевле — написать два шага sifting вручную поверх сырых указателей (без функционального компаратора). Или хранить индекс максимума и делать отложенный reheap каждые k итераций (k≈8–16).
Prefetch: перед std::pop_heap префетчьте вершину следующего крупного интервала R[i] (эвристика: i=1/2) через _mm_prefetch((const char*)R[i], _MM_HINT_T0).
Векторные константы: в FABE13_SINCOS ваши const __m256 создаются один раз на итерацию while — это ок. Если функция зовётся часто малыми порциями, пометьте их как static const (они попадут в .rdata и не будут пересоздаваться).
Stack → TLS для маленьких фиксированных массивов: angles[14], sin/cos[8] можно сделать static thread_local alignas(64) float[...] — меньше давления на stack frame и prolog/epilog.
NUMA: заколите worker-потоки к узлам (SetThreadGroupAffinity/SetThreadIdealProcessor), и вызывайте preallocate_for_depth уже из каждого worker-потока (первая «первое касание» будет локальной).
Huge pages привилегии: если есть доступ — поднимите SeLockMemoryPrivilege для процесса, тогда ветка MEM_LARGE_PAGES реально сработает (минус TLB-промахи на больших деревьях).
Выравнивание CoTaskMem: финальный буфер под результат отдаётся через CoTaskMemAlloc. Там выравнивание обычно 16. Если потребуются AVX-aligned записи — сначала пишите в ваш ареновый буфер, потом memcpy в CoTaskMem (один крупный копик — быстрее, чем работать «не по выравниванию» весь путь).

Микро-правки прямо в вашем коде

Везде, где сейчас __declspec(align(16)) в стеке, замените на alignas(64).
В PreallocatedPeanoMemoryManager::MemoryBlock::initialize — выравнивание минимум 64, как в коде выше.
Добавьте operator new/delete у Interval (уже показал).
Если точно знаете, что in/sin_out/cos_out выровнены — заведите вторую версию FABE13_SINCOS_ALIGNED и зовите её там, где данные — ваши (после этого компилятор сгенерит чистые vmovaps вместо vmovups).

Если хочется, могу выдать готовый патч именно под ваш файл: заменю менеджер памяти, добавлю перегрузки new для Interval, протащу alignas(64) и вставлю _ALIGNED-вариант SINCOS.

пишу крайне высокопроизводительное приложение став...

سؤال

الفكرة

Optimizing memory allocation

Optimizing with AVX2

Optimizing allocation and performance

Performance optimizations

Optimizing performance

Optimizing memory management

Optimizing memory and performance

الإجابة