OpenMP SIMD Euclidean Distance I'm using cv::norm() to compute the L2 distance between row vectors in thousands of dimensions (at least 4096 dimensions). I'm using an AVX2 machine and icpc version 2017. 

Since I'm computing this function millions of times, I need the possible most efficient version for this function and I highly doubt that cv::norm() exploits vectorization properly.

So I implemnted this function:

#include &lt;opencv2/opencv.hpp&gt;

float euclDist (const cv::Mat1f &amp;first, const cv::Mat1f &amp;second){
    float dist = 0;
    const float *firstp = first.ptr&lt;float&gt;(0);
    const float *secondp = second.ptr&lt;float&gt;(0);
    #pragma omp simd reduction(+:dist)
    for(size_t i=0; i&lt;first.cols; i++)
        dist += std::pow(firstp[i]-secondp[i],2);
    return std::sqrt(dist);
}

int main (int argc, const char *argv[]){
    int codeSize = 4096;
    cv::Mat1f mat(1, codeSize); 
    float low = 0;
    float high = 2.0;
     srand (0);
     cv::theRNG().state = cv::getTickCount() ;
    randu(mat, cv::Scalar(low), cv::Scalar(high));
    euclDist(mat, mat);
    return 0;
}


The example above is similar what the actually program does: it generates a row vector and it computes the distance with itself (i.e. 0). In the actual program generates a matrix and it randomly compute the distance between some rows. I don't think that I can get only that part of the code easily.

I already noticed an improvement w.r.t. cv::norm(), but scales badly when computing multiple distances in parallel. According to Intel VTune I fully use all the processors, but I get poor memory performance, as it shown from this general exploration analysis:



This is the result from: sudo lshw -short -C memory

H/W path       Device      Class          Description
=====================================================
/0/0                       memory         128KiB BIOS
/0/4/b                     memory         32KiB L1 cache
/0/4/c                     memory         256KiB L2 cache
/0/4/d                     memory         6MiB L3 cache
/0/a                       memory         32KiB L1 cache
/0/2a                      memory         12GiB System Memory
/0/2a/0                    memory         4GiB SODIMM DDR3 Synchronous 1600 MHz (0.6 ns)
/0/2a/1                    memory         DIMM [empty]
/0/2a/2                    memory         8GiB SODIMM DDR3 Synchronous 1600 MHz (0.6 ns)
/0/2a/3                    memory         DIMM [empty]


What I could do?
I had to adapt the test program to get meaningful results using /usr/bin/time:

#include &lt;iostream&gt;
int main()
{
    static const int codeSize = 4096;
    cv::Mat1f mat1(1, codeSize);
    cv::Mat1f mat2(1, codeSize);
    float low = 0;
    float high = 2.0;
    srand(0);
    cv::theRNG().state = cv::getTickCount() ;

    randu(mat1, cv::Scalar(low), cv::Scalar(high));
    randu(mat2, cv::Scalar(low), cv::Scalar(high));

    double d = 0;
    for (auto i = 0;  i &lt; 1000000;  ++i) {
        d += euclDist(mat1, mat2);
    }

    std::cout &lt;&lt; d &lt;&lt; std::endl;
}


(I used two different inputs to avoid a predictable result of 0, made sure I used the result to prevent optimising out, and executed a million iterations to get the run time over a second).

Here's my compilation flags (used with my standard Makefile):


167501: PKGS += opencv
167501: CXXFLAGS += -O3 -march=native
167501: CXXFLAGS += -fopenmp




Having done that, I tried a couple of variations, all compiled with GCC version 7.1.0, and executed them on an Intel i7-6700K with 128KiB/1MiB/8MiB L1d/2/3 caches.

The original clocks in at 14.6 seconds user time without the OpenMP directive and 4.5 seconds with it.

Changing the std::pow() call to a plain multiply improved those times to 4.2 and 0.6 seconds respectively.

I got another factor of three improvement in wall-clock time, at the expense of much more processor time, by distributing across 8 cores with

#pragma omp parallel for simd reduction(+:dist)


For this code, the loop is a bit small for parallelization, it seems.



float euclDist(const cv::Mat1f&amp; first, const cv::Mat1f&amp; second)
{
    float dist = 0;
    const float *const firstp = first.ptr&lt;float&gt;(0);
    const float *const secondp = second.ptr&lt;float&gt;(0);
#pragma omp simd reduction(+:dist)
    for(auto i = 0;  i &lt; first.cols;  ++i) {
        auto const d = firstp[i] - secondp[i];
        dist += d * d;
    }
    return std::sqrt(dist);
}