|| || || || || ||

---

Chi2Distance.h

Go to the documentation of this file.

#ifndef Impala_Core_Vector_Chi2Distance_h
#define Impala_Core_Vector_Chi2Distance_h

#include "Core/Vector/VectorTem.h"
#ifdef SSE_USED
#include <xmmintrin.h>
#include <stdint.h>
#ifndef POINTER_ALIGNED
#define POINTER_ALIGNED(x) (!((( intptr_t)x) & 0xF))
#endif
#endif // SSE_USED

namespace Impala
{
namespace Core
{
namespace Vector
{


inline double
Chi2Distance(const Vector::VectorTem<double>& v1,
             const Vector::VectorTem<double>& v2)
{
    int length = v1.Size();
    double accumulator = 0;
    for (int i=0 ; i<length ; i++)
    {
        double sum = v1[i] + v2[i];
        if (sum)
            accumulator += (v1[i] - v2[i]) * (v1[i] - v2[i]) / sum;
    }
    // make it between 0..1, instead of 0..2 (this statement is only true if
    // the input vectors sum to one)
    return accumulator / 2.0;
}


#ifdef SSE_USED

// assumption: vectorSize is a multiple of 4 (or at least the memory is big
// enough/zero-padded), pointers are 16-byte aligned
inline void
DoSlabChi2DistanceSSE(float* C, const float* A, const float* B,
                      unsigned int slabSizeA, unsigned int slabSizeB,
                      unsigned int vectorSize)
{
    const unsigned int alignFactor = 4;
    const int SSELength = IntAlignUp(vectorSize, alignFactor) / alignFactor;
    const __m128 zero = _mm_setzero_ps();
    const __m128 two = _mm_set1_ps(2.0f);
    #pragma omp parallel for
    for (unsigned int a = 0; a < slabSizeA; a++)
    {
        const __m128* baseA = (__m128*) (A + (a * vectorSize));
        for (unsigned int b = 0; b < slabSizeB; b++)
        {
            const __m128* baseB = (__m128*) (B + (b * vectorSize));
            __m128 sum = _mm_setzero_ps();
            for (unsigned int k = 0; k < SSELength; k++)
            {
                __m128 As = baseA[k];
                __m128 Bs = baseB[k];
                __m128 tmp = _mm_sub_ps(As, Bs);
                // sum += (As - Bs) * (As - Bs) / (As + Bs);
                tmp =_mm_div_ps(_mm_mul_ps(tmp, tmp), _mm_add_ps(As, Bs));
                // fix the NaN for 0/0
                tmp = _mm_max_ps(tmp, zero);
                sum = _mm_add_ps(sum, tmp);
            }
            sum = _mm_div_ps(sum, two);
            // add the components of sum together
            // shuffle = [1 0 3 2]
            // sum     = [3+1 2+0 1+3 0+2]
            // shuffle = [2+0 3+1 0+2 1+3]
            // res     = [3+1+2+0 3+1+2+0 1+3+0+2 0+2+1+3]
            __m128 shuffle;
            shuffle = _mm_shuffle_ps(sum, sum, _MM_SHUFFLE(1, 0, 3, 2));
            sum     = _mm_add_ps(sum, shuffle) ;
            shuffle = _mm_shuffle_ps(sum, sum, _MM_SHUFFLE(2, 3, 0, 1));
            __m128 result  = _mm_add_ps(sum, shuffle) ;
            C[b * slabSizeA + a] = _mm_cvtss_f32(result);
        }
    }
}

// assumption: vectorSize is a multiple of 4 (or at least the memory is big
// enough/zero-padded), pointers are 16-byte aligned
inline void
DoSlabChi2DistanceSSE(double* C, const double* A, const double* B,
                      unsigned int slabSizeA, unsigned int slabSizeB,
                      unsigned int vectorSize)
{
    const unsigned int alignFactor = 2;
    const int SSELength = IntAlignUp(vectorSize, alignFactor) / alignFactor;
    const __m128d zero = _mm_setzero_pd();
    const __m128d two = _mm_set1_pd(2.0);
    #pragma omp parallel for
    for (unsigned int a = 0; a < slabSizeA; a++)
    {
        const __m128d* baseA = (__m128d*) (A + (a * vectorSize));
        for (unsigned int b = 0; b < slabSizeB; b++)
        {
            const __m128d* baseB = (__m128d*) (B + (b * vectorSize));
            __m128d sum = _mm_setzero_pd();
            for (unsigned int k = 0; k < SSELength; k++)
            {
                __m128d As = baseA[k];
                __m128d Bs = baseB[k];
                __m128d tmp = _mm_sub_pd(As, Bs);
                // sum += (As - Bs) * (As - Bs) / (As + Bs);
                tmp =_mm_div_pd(_mm_mul_pd(tmp, tmp), _mm_add_pd(As, Bs));
                // fix the NaN for 0/0
                tmp = _mm_max_pd(tmp, zero);
                sum = _mm_add_pd(sum, tmp);
            }
            sum = _mm_div_pd(sum, two);
            // add the components of sum together
            __m128d shuffle = _mm_shuffle_pd(sum, sum, _MM_SHUFFLE2(0, 1));
            __m128d result  = _mm_add_pd(sum, shuffle) ;
            C[b * slabSizeA + a] = _mm_cvtsd_f64(result);
        }
    }
}

#endif // SSE_USED


template <class FType>
void
DoSlabChi2DistanceStd(FType* C, const FType* A, const FType* B,
                      unsigned int slabSizeA, unsigned int slabSizeB,
                      unsigned int vectorSize)
{
    #pragma omp parallel for
    for (unsigned int a = 0; a < slabSizeA; a++)
    {
        for (unsigned int b = 0; b < slabSizeB; b++)
        {
            const FType* baseA = A + (a * vectorSize);
            const FType* baseB = B + (b * vectorSize);
            FType sum = 0;
            for (unsigned int k = 0; k < vectorSize; k++)
            {
                FType As = baseA[k];
                FType Bs = baseB[k];
                if(As == 0)
                {
                    sum += Bs;
                }
                else
                {
                    sum += (As - Bs) * (As - Bs) / (As + Bs);
                }
            }
            C[b * slabSizeA + a] = (FType)(sum / 2);
        }
    }
}

template <class FType>
inline void
DoSlabChi2Distance(FType* C, const FType* A, const FType* B,
                   unsigned int slabSizeA, unsigned int slabSizeB,
                   unsigned int vectorSize)
{
#ifdef SSE_USED
    return DoSlabChi2DistanceSSE(C, A, B, slabSizeA, slabSizeB, vectorSize);
#endif
    DoSlabChi2DistanceStd(C, A, B, slabSizeA, slabSizeB, vectorSize);
}

} // namespace Vector
} // namespace Core
} // namespace Impala

#endif

---