struct timeval tv;

gettimeofday (&tv, nullptr);

return tv.tv_sec * 1e3 + tv.tv_usec * 1e-3;

{

int pid = getpid ();

char fname[256];

snprintf (fname, 256, "/proc/%d/status", pid);

FILE * f = fopen (fname, "r");

FAISS_ASSERT (f || !"cannot open proc status file") ;

size_t sz = 0;

for (;;) {

char buf [256];

if (!fgets (buf, 256, f)) break;

if (sscanf (buf, "VmRSS: %ld kB", &sz) == 1) break;

}

fclose (f);

return sz;

{

fprintf(stderr, "WARN: get_mem_usage_kb not implemented on the mac\n");

return 0;

{

int32_t a;

random_r (&rand_data, &a);

return a;

{

int32_t a, b;

random_r (&rand_data, &a);

random_r (&rand_data, &b);

return long(a) | long(b) << 31;

{

memset (&rand_data, 0, sizeof (rand_data));

initstate_r (seed, rand_state, sizeof (rand_state), &rand_data);

{

memcpy (rand_state, other.rand_state, sizeof(rand_state));

rand_data = other.rand_data;

setstate_r (rand_state, &rand_data);

{

// RAND_MAX is 31 bits

// try to add more randomness in the lower bits

int lowbits = rand_r(&rand_state) >> 15;

return rand_r(&rand_state) ^ lowbits;

{ // this suffers form non-uniform probabilities when max is not a

// power of 2, but if RAND_MAX >> max the bias is limited.

return rand_int () % max;

{

// only try to parallelize on large enough arrays

const size_t nblock = n < 1024 ? 1 : 1024;

RandomGenerator rng0 (seed);

int a0 = rng0.rand_int (), b0 = rng0.rand_int ();

#pragma omp parallel for

for (size_t j = 0; j < nblock; j++) {

RandomGenerator rng (a0 + j * b0);

const size_t istart = j * n / nblock;

const size_t iend = (j + 1) * n / nblock;

for (size_t i = istart; i < iend; i++)

x[i] = rng.rand_float ();

}

{

// only try to parallelize on large enough arrays

const size_t nblock = n < 1024 ? 1 : 1024;

RandomGenerator rng0 (seed);

int a0 = rng0.rand_int (), b0 = rng0.rand_int ();

#pragma omp parallel for

for (size_t j = 0; j < nblock; j++) {

RandomGenerator rng (a0 + j * b0);

double a = 0, b = 0, s = 0;

int state = 0; /* generate two number per "do-while" loop */

const size_t istart = j * n / nblock;

const size_t iend = (j + 1) * n / nblock;

for (size_t i = istart; i < iend; i++) {

/* Marsaglia's method (see Knuth) */

if (state == 0) {

do {

a = 2.0 * rng.rand_double () - 1;

b = 2.0 * rng.rand_double () - 1;

s = a * a + b * b;

} while (s >= 1.0);

x[i] = a * sqrt(-2.0 * log(s) / s);

}

else

x[i] = b * sqrt(-2.0 * log(s) / s);

state = 1 - state;

}

}

{

// only try to parallelize on large enough arrays

const size_t nblock = n < 1024 ? 1 : 1024;

RandomGenerator rng0 (seed);

int a0 = rng0.rand_int (), b0 = rng0.rand_int ();

#pragma omp parallel for

for (size_t j = 0; j < nblock; j++) {

RandomGenerator rng (a0 + j * b0);

const size_t istart = j * n / nblock;

const size_t iend = (j + 1) * n / nblock;

for (size_t i = istart; i < iend; i++)

x[i] = rng.rand_long ();

}

{

for (size_t i = 0; i < n; i++) perm[i] = i;

RandomGenerator rng (seed);

for (size_t i = 0; i + 1 < n; i++) {

int i2 = i + rng.rand_int (n - i);

std::swap(perm[i], perm[i2]);

}

{

// only try to parallelize on large enough arrays

const size_t nblock = n < 1024 ? 1 : 1024;

RandomGenerator rng0 (seed);

int a0 = rng0.rand_int (), b0 = rng0.rand_int ();

#pragma omp parallel for

for (size_t j = 0; j < nblock; j++) {

RandomGenerator rng (a0 + j * b0);

const size_t istart = j * n / nblock;

const size_t iend = (j + 1) * n / nblock;

size_t i;

for (i = istart; i < iend; i++)

x[i] = rng.rand_long ();

}

float * __restrict x,

size_t n, size_t d, size_t nu)

{

size_t i, j, l;

for (i = 0; i < n; i++) {

const float * up = u;

for (l = 0; l < nu; l++) {

float ip1 = 0, ip2 = 0;

for (j = 0; j < d; j+=2) {

ip1 += up[j] * x[j];

ip2 += up[j+1] * x[j+1];

}

float ip = 2 * (ip1 + ip2);

for (j = 0; j < d; j++)

x[j] -= ip * up[j];

up += d;

}

x += d;

}

{

size_t i, j, l;

for (i = 0; i < n; i++) {

const float * up = u;

for (l = 0; l < nu; l++) {

double ip = 0;

for (j = 0; j < d; j++)

ip += up[j] * x[j];

ip *= 2;

for (j = 0; j < d; j++)

x[j] -= ip * up[j];

up += d;

}

x += d;

}

{

// check which values to clear out

__m128i d4 = _mm_set1_epi32 (d);

__m128i lim = _mm_set_epi32 (3, 2, 1, 0);

__m128i mask = _mm_cmpgt_epi32 (d4, lim);

return _mm_and_ps (x, (__m128)mask);

const float * y,

size_t d)

{

__m128 mx, my;

__m128 msum1 = _mm_setzero_ps();

while (d > 4) {

mx = _mm_loadu_ps (x); x += 4;

my = _mm_loadu_ps (y); y += 4;

const __m128 a_m_b1 = _mm_sub_ps (mx, my);

msum1 = _mm_add_ps (msum1, _mm_mul_ps (a_m_b1, a_m_b1));

d -= 4;

}

// add the last 1, 2, 3, or 4 values

mx = _mm_loadu_ps (x);

my = _mm_loadu_ps (y);

__m128 a_m_b1 = maskout (d, _mm_sub_ps (mx, my));

msum1 = _mm_add_ps (msum1, _mm_mul_ps (a_m_b1, a_m_b1));

msum1 = _mm_hadd_ps (msum1, msum1);

msum1 = _mm_hadd_ps (msum1, msum1);

return _mm_cvtss_f32 (msum1);

const float * y,

size_t d)

{

size_t i;

float res_ = 0;

for (i = 0; i < d; i++) {

const float tmp = x[i] - y[i];

res_ += tmp * tmp;

}

return res_;

const float * y,

size_t d)

{

__m128 mx, my;

__m128 msum1 = _mm_setzero_ps();

while (d > 4) {

mx = _mm_loadu_ps (x); x += 4;

my = _mm_loadu_ps (y); y += 4;

msum1 = _mm_add_ps (msum1, _mm_mul_ps (mx, my));

d -= 4;

}

// add the last 1, 2, 3, or 4 values

mx = _mm_loadu_ps (x);

my = _mm_loadu_ps (y);

__m128 prod = maskout (d, _mm_mul_ps (mx, my));

msum1 = _mm_add_ps (msum1, prod);

msum1 = _mm_hadd_ps (msum1, msum1);

msum1 = _mm_hadd_ps (msum1, msum1);

return _mm_cvtss_f32 (msum1);

const float * y,

size_t d)

{

size_t i;

float res_ = 0;

for (i = 0; i < d; i++)

res_ += x[i] * y[i];

return res_;

size_t d)

{

__m128 mx;

__m128 msum1 = _mm_setzero_ps();

while (d > 4) {

mx = _mm_loadu_ps (x); x += 4;

msum1 = _mm_add_ps (msum1, _mm_mul_ps (mx, mx));

d -= 4;

}

mx = maskout (d, _mm_loadu_ps (x));

msum1 = _mm_add_ps (msum1, _mm_mul_ps (mx, mx));

msum1 = _mm_hadd_ps (msum1, msum1);

msum1 = _mm_hadd_ps (msum1, msum1);

return _mm_cvtss_f32 (msum1);

size_t d)

{

size_t i;

double res_ = 0;

for (i = 0; i < d; i++)

res_ += x[i] * x[i];

return res_;

const float * x,

const float * y,

size_t d, size_t ny)

{

for (size_t i = 0; i < ny; i++) {

ip[i] = fvec_inner_product (x, y, d);

y += d;

}

const float * x,

const float * y,

size_t d, size_t ny)

{

for (size_t i = 0; i < ny; i++) {

dis[i] = fvec_L2sqr (x, y, d);

y += d;

}

const float * __restrict x,

size_t d, size_t nx)

{

#pragma omp parallel for

for (size_t i = 0; i < nx; i++) {

nr[i] = sqrtf (fvec_norm_L2sqr (x + i * d, d));

}

const float * __restrict x,

size_t d, size_t nx)

{

#pragma omp parallel for

for (size_t i = 0; i < nx; i++)

nr[i] = fvec_norm_L2sqr (x + i * d, d);

{

#pragma omp parallel for

for (size_t i = 0; i < nx; i++) {

float * __restrict xi = x + i * d;

float nr = fvec_norm_L2sqr (xi, d);

if (nr > 0) {

size_t j;

const float inv_nr = 1.0 / sqrtf (nr);

for (j = 0; j < d; j++)

xi[j] *= inv_nr;

}

}

const float * y,

size_t d, size_t nx, size_t ny,

float_minheap_array_t * res)

{

size_t k = res->k;

#pragma omp parallel for

for (size_t i = 0; i < nx; i++) {

const float * x_ = x + i * d;

const float * y_ = y;

float * __restrict simi = res->get_val(i);

long * __restrict idxi = res->get_ids (i);

minheap_heapify (k, simi, idxi);

for (size_t j = 0; j < ny; j++) {

float ip = fvec_inner_product (x_, y_, d);

if (ip > simi[0]) {

minheap_pop (k, simi, idxi);

minheap_push (k, simi, idxi, ip, j);

}

y_ += d;

}

minheap_reorder (k, simi, idxi);

}

const float * x,

const float * y,

size_t d, size_t nx, size_t ny,

float_maxheap_array_t * res)

{

size_t k = res->k;

#pragma omp parallel for

for (size_t i = 0; i < nx; i++) {

const float * x_ = x + i * d;

const float * y_ = y;

size_t j;

float * __restrict simi = res->get_val(i);

long * __restrict idxi = res->get_ids (i);

maxheap_heapify (k, simi, idxi);

for (j = 0; j < ny; j++) {

float disij = fvec_L2sqr (x_, y_, d);

if (disij < simi[0]) {

maxheap_pop (k, simi, idxi);

maxheap_push (k, simi, idxi, disij, j);

}

y_ += d;

}

maxheap_reorder (k, simi, idxi);

}

const float * x,

const float * y,

size_t d, size_t nx, size_t ny,

float_minheap_array_t * res)

{

res->heapify ();

// BLAS does not like empty matrices

if (nx == 0 || ny == 0) return;

/* block sizes */

const size_t bs_x = 4096, bs_y = 1024;

// const size_t bs_x = 16, bs_y = 16;

float *ip_block = new float[bs_x * bs_y];

for (size_t i0 = 0; i0 < nx; i0 += bs_x) {

size_t i1 = i0 + bs_x;

if(i1 > nx) i1 = nx;

for (size_t j0 = 0; j0 < ny; j0 += bs_y) {

size_t j1 = j0 + bs_y;

if (j1 > ny) j1 = ny;

/* compute the actual dot products */

{

float one = 1, zero = 0;

FINTEGER nyi = j1 - j0, nxi = i1 - i0, di = d;

sgemm_ ("Transpose", "Not transpose", &nyi, &nxi, &di, &one,

y + j0 * d, &di,

x + i0 * d, &di, &zero,

ip_block, &nyi);

}

/* collect maxima */

res->addn (j1 - j0, ip_block, j0, i0, i1 - i0);

}

}

delete [] ip_block;

res->reorder ();

const float * y,

size_t d, size_t nx, size_t ny,

float_maxheap_array_t * res,

const DistanceCorrection &corr)

{

res->heapify ();

// BLAS does not like empty matrices

if (nx == 0 || ny == 0) return;

size_t k = res->k;

/* block sizes */

const size_t bs_x = 4096, bs_y = 1024;

// const size_t bs_x = 16, bs_y = 16;

float *ip_block = new float[bs_x * bs_y];

float *x_norms = new float[nx];

fvec_norms_L2sqr (x_norms, x, d, nx);

float *y_norms = new float[ny];

fvec_norms_L2sqr (y_norms, y, d, ny);

for (size_t i0 = 0; i0 < nx; i0 += bs_x) {

size_t i1 = i0 + bs_x;

if(i1 > nx) i1 = nx;

for (size_t j0 = 0; j0 < ny; j0 += bs_y) {

size_t j1 = j0 + bs_y;

if (j1 > ny) j1 = ny;

/* compute the actual dot products */

{

float one = 1, zero = 0;

FINTEGER nyi = j1 - j0, nxi = i1 - i0, di = d;

sgemm_ ("Transpose", "Not transpose", &nyi, &nxi, &di, &one,

y + j0 * d, &di,

x + i0 * d, &di, &zero,

ip_block, &nyi);

}

/* collect minima */

#pragma omp parallel for

for (size_t i = i0; i < i1; i++) {

float * __restrict simi = res->get_val(i);

long * __restrict idxi = res->get_ids (i);

const float *ip_line = ip_block + (i - i0) * (j1 - j0);

for (size_t j = j0; j < j1; j++) {

float ip = *ip_line++;

float dis = x_norms[i] + y_norms[j] - 2 * ip;

dis = corr (dis, i, j);

if (dis < simi[0]) {

maxheap_pop (k, simi, idxi);

maxheap_push (k, simi, idxi, dis, j);

}

}

}

}

}

res->reorder ();

delete [] ip_block;

delete [] x_norms;

delete [] y_norms;

const float * y,

size_t d, size_t nx, size_t ny,

float_minheap_array_t * res)

{

if (d % 4 == 0 && nx < 20) {

knn_inner_product_sse (x, y, d, nx, ny, res);

} else {

knn_inner_product_blas (x, y, d, nx, ny, res);

}

float operator()(float dis, size_t qno, size_t bno) const {

return dis;

}

const float * y,

size_t d, size_t nx, size_t ny,

float_maxheap_array_t * res)

{

if (d % 4 == 0 && nx < 20) {

knn_L2sqr_sse (x, y, d, nx, ny, res);

} else {

NopDistanceCorrection nop;

knn_L2sqr_blas (x, y, d, nx, ny, res, nop);

}

const float *base_shift;

float operator()(float dis, size_t qno, size_t bno) const {

return dis - base_shift[bno];

}

const float * x,

const float * y,

size_t d, size_t nx, size_t ny,

float_maxheap_array_t * res,

const float *base_shift)

{

BaseShiftDistanceCorrection corr = {base_shift};

knn_L2sqr_blas (x, y, d, nx, ny, res, corr);

const float * x,

const float * y,

const long * __restrict ids, /* for y vecs */

size_t d, size_t nx, size_t ny)

{

#pragma omp parallel for

for (size_t j = 0; j < nx; j++) {

const long * __restrict idsj = ids + j * ny;

const float * xj = x + j * d;

float * __restrict ipj = ip + j * ny;

for (size_t i = 0; i < ny; i++) {

if (idsj[i] < 0)

continue;

ipj[i] = fvec_inner_product (xj, y + d * idsj[i], d);

}

}

const float * x,

const float * y,

const long * __restrict ids, /* ids of y vecs */

size_t d, size_t nx, size_t ny)

{

#pragma omp parallel for

for (size_t j = 0; j < nx; j++) {

const long * __restrict idsj = ids + j * ny;

const float * xj = x + j * d;

float * __restrict disj = dis + j * ny;

for (size_t i = 0; i < ny; i++) {

if (idsj[i] < 0)

continue;

disj[i] = fvec_L2sqr (xj, y + d * idsj[i], d);

}

}

const float * y,

const long * __restrict ids,

size_t d, size_t nx, size_t ny,

float_minheap_array_t * res)

{

size_t k = res->k;

#pragma omp parallel for

for (size_t i = 0; i < nx; i++) {

const float * x_ = x + i * d;

const long * idsi = ids + i * ny;

size_t j;

float * __restrict simi = res->get_val(i);

long * __restrict idxi = res->get_ids (i);

minheap_heapify (k, simi, idxi);

for (j = 0; j < ny; j++) {

if (idsi[j] < 0) break;

float ip = fvec_inner_product (x_, y + d * idsi[j], d);

if (ip > simi[0]) {

minheap_pop (k, simi, idxi);

minheap_push (k, simi, idxi, ip, idsi[j]);

}

}

minheap_reorder (k, simi, idxi);

}

const float * y,

const long * __restrict ids,

size_t d, size_t nx, size_t ny,

float_maxheap_array_t * res)

{

size_t k = res->k;

#pragma omp parallel for

for (size_t i = 0; i < nx; i++) {

const float * x_ = x + i * d;

const long * __restrict idsi = ids + i * ny;

float * __restrict simi = res->get_val(i);

long * __restrict idxi = res->get_ids (i);

maxheap_heapify (res->k, simi, idxi);

for (size_t j = 0; j < ny; j++) {

float disij = fvec_L2sqr (x_, y + d * idsi[j], d);

if (disij < simi[0]) {

maxheap_pop (k, simi, idxi);

maxheap_push (k, simi, idxi, disij, idsi[j]);

}

}

maxheap_reorder (res->k, simi, idxi);

}