test_euclid_ask.h (only need to read 2 functions: euclid_slow, euclid_fast)
#pragma once
#include "included.h"
double
euclid_slow(int n, double* data1, double* data2, int* mask1, int* mask2, const double weight[])
{
double result = 0.0;
double totalWeight = 0;
for (int i = 0; i < n; i++) {
if (mask1[i] && mask2[i]) {
double term = data1[i] - data2[i];
result += weight[i] * term * term;
totalWeight += weight[i];
}
}
if (totalWeight==0) return 0;
return result / totalWeight;
}
double
euclid_fast(int n, double* data1, double* data2, int* mask1, int* mask2, const double weight[])
{
double result = 0.0;
double totalWeight = 0;
double subResult[4] = { 0. };
double subTweight[4] = { 0. };
double subDiff[4] = { 0. };
double subWeight[4] = { 0. };
double subMask[4] = { 0. };
int nstep4 = n - n % 4;
for (int i = 0; i < nstep4; i += 4) {
subMask[0] = mask1[i] && mask2[i];
subMask[1] = mask1[i + 1] && mask2[i + 1];
subMask[2] = mask1[i + 2] && mask2[i + 2];
subMask[3] = mask1[i + 3] && mask2[i + 3];
if (!(subMask[0] || subMask[1] || subMask[2] || subMask[3])) continue;
subDiff[0] = data1[i] - data2[i];
subDiff[1] = data1[i + 1] - data2[i + 1];
subDiff[2] = data1[i + 2] - data2[i + 2];
subDiff[3] = data1[i + 3] - data2[i + 3];
subDiff[0] *= subDiff[0];
subDiff[1] *= subDiff[1];
subDiff[2] *= subDiff[2];
subDiff[3] *= subDiff[3];
subWeight[0] = weight[i] * subMask[0];
subWeight[1] = weight[i + 1] * subMask[1];
subWeight[2] = weight[i + 2] * subMask[2];
subWeight[3] = weight[i + 3] * subMask[3];
subTweight[0] += subWeight[0];
subTweight[1] += subWeight[1];
subTweight[2] += subWeight[2];
subTweight[3] += subWeight[3];
subResult[0] += subWeight[0] * subDiff[0];
subResult[1] += subWeight[1] * subDiff[1];
subResult[2] += subWeight[2] * subDiff[2];
subResult[3] += subWeight[3] * subDiff[3];
}
for (int i = nstep4; i < n; i++) {
if (mask1[i] && mask2[i]) {
double term = data1[i] - data2[i];
result += weight[i] * term * term;
totalWeight += weight[i];
}
}
result += subResult[0] + subResult[1] + subResult[2] + subResult[3];
totalWeight += subTweight[0] + subTweight[1] + subTweight[2] + subTweight[3];
//cout << "end fast\n";
if (!totalWeight) return 0;
return result / totalWeight;
}
void test_euclid_ask()
{
const int MAXN = 10000000, MINN = 1000000;
double* data1, * data2;
int* mask1, * mask2;
double* dataPro1, * dataPro2;
int* maskPro1, * maskPro2;
double *weight, * weightPro;
//***********
data1 = new double[MAXN + MINN + 1];
data2 = new double[MAXN + MINN + 1];
mask1 = new int[MAXN + MINN + 1];
mask2 = new int[MAXN + MINN + 1];
dataPro1 = new double[MAXN + MINN + 1];
dataPro2 = new double[MAXN + MINN + 1];
maskPro1 = new int[MAXN + MINN + 1];
maskPro2 = new int[MAXN + MINN + 1];
// ******
weight = new double[MAXN + MINN + 1];
weightPro = new double[MAXN + MINN + 1];
MyTimer timer;
int n;
double guess1, guess2, tmp, total1 = 0, total2 = 0, prev1 = 0, prev2 = 0;
for (int t = 5000; t < 6000; t++) {
if (t <= 5000) n = t;
else n = MINN + rand() % (MAXN - MINN);
cout << n << "\n";
int index = 0;
for (int i = 0; i < n; i++) {
weight[i] = int64(randomed()) % 100;
data1[i] = int64(randomed()) % 100;
data2[i] = int64(randomed()) % 100;
mask1[i] = rand() % 10;
mask2[i] = rand() % 10;
}
memcpy(weightPro, weight, n * sizeof(double));
memcpy(dataPro1, data1, n * sizeof(double));
memcpy(dataPro2, data2, n * sizeof(double));
memcpy(maskPro1, mask1, n * sizeof(int));
memcpy(maskPro2, mask2, n * sizeof(int));
//****
int tmp = flush_cache(); // do something to ensure the cache does not contain test data
cout << "ignore this " << tmp << "\n";
timer.startCounter();
guess1 = euclid_slow(n, data1, data2, mask1, mask2, weight);
tmp = timer.getCounterMicro();
total1 += tmp;
cout << "time slow = " << tmp << " us\n";
timer.startCounter();
guess2 = euclid_fast(n, dataPro1, dataPro2, maskPro1, maskPro2, weightPro);
tmp = timer.getCounterMicro();
total2 += tmp;
cout << "time fast = " << tmp << " us\n";
bool ok = fabs(guess1 - guess2) <= 0.1;
if (!ok) {
cout << "error at N = " << n << "\n";
exit(-1);
}
cout << "\n";
}
cout << "slow speed = " << (total1 / 1000) << " ms\n";
cout << "fast speed = " << (total2 / 1000) << " ms\n";
}
Basically, the function computes a kind-of Euclidean distance between 2 arrays:
result = sum(weight[i] * (data1[i] - data2[i])^2)
but only in positions where both values are available (mask1[i]==0
means it's ignored, same with mask2
). The normal code is in function euclid_slow.
So I tried to improve the code by processing 4 elements at once, hoping that SSE/AVX can speed this up. However, the result stays the same or slower(using g++ -O3 -march=native) or becomes 40% slower (using Visual Studio 2019 compiler, release mode (x64), -O2, AVX2 enabled). I tried both -O2
and -O3
, same result.
The compiler made better optimizations than what I wrote. But how can I make it actually faster?
Edit: code to test the programs here
-ffast-math
or an OpenMP pragma to auto-vectorize the reductions intoresult
andtotalWeight
.subMask[0] = mask1[i] && mask2[i];
are hard to vectorize, sincemask2[i]
must not be accessed ifmask1[i] == 0
. Also,subWeight[0] = weight[i] * subMask[0];
could result in a NaN ifweight
contains infinities, so the compiler can't replace it by a more efficient bit-and operation. As @Peter said: Either don't manually unroll anything but let the compiler do it (usually requires-ffast-math
, or a subset from that), or completely vectorize manually with intrinsics.