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

---

MatNorm2DistSSE.h

Go to the documentation of this file.

#ifndef Impala_Core_Matrix_MatNorm2DistSSE_h
#define Impala_Core_Matrix_MatNorm2DistSSE_h

#include "Core/Matrix/MatFunc.h"
#include "Core/Matrix/MatSum.h"
#include "Core/Matrix/MatMul.h"
#include "Core/Array/Abs.h"
#include "Core/Array/Add.h"
#include "Core/Array/Mul.h"
#include "Core/Array/MulVal.h"
#include "Core/Array/Sqrt.h"
#include "Basis/Timer.h"

#include <xmmintrin.h>
#include <stdint.h>
#define POINTER_ALIGNED(x) (!((( intptr_t)x) & 0xF))

namespace Impala
{
namespace Core
{
namespace Matrix
{

inline Mat*
SIMD_MatSquareAndSumAxis1(Mat* m)
{
    int n = MatNrRow(m);
    Mat* res = MatCreate<Mat>(n, 1);
    for (int i=0 ; i<n ; i++)
    {
        __m128d tmp = _mm_setzero_pd();
        __m128d* base = (__m128d*) MatE(m, i, 0);
        for (int j=0 ; j<IntAlignDown(MatNrCol(m), 2) / 2 ; j++)
            tmp = _mm_add_pd(tmp, _mm_mul_pd(base[j], base[j]));

        __m128d shuffle = _mm_shuffle_pd(tmp, tmp, _MM_SHUFFLE2(0, 1));
        Real64 result = _mm_cvtsd_f64(_mm_add_pd(tmp, shuffle));
        for (int j=IntAlignDown(MatNrCol(m), 2); j<MatNrCol(m) ; j++)
            result += result + *MatE(m, i, j);
        *MatE(res, i, 0) = result;
    }
    return res;
}

inline Mat32*
SIMD_MatSquareAndSumAxis1(Mat32* m)
{
    int n = MatNrRow(m);
    Mat32* res = MatCreate<Mat32>(n, 1);
    for (int i=0 ; i<n ; i++)
    {
        __m128 tmp = _mm_setzero_ps();
        __m128* base = (__m128*) MatE(m, i, 0);
        for (int j=0 ; j<IntAlignDown(MatNrCol(m), 4) / 4 ; j++)
            tmp = _mm_add_ps(tmp, _mm_mul_ps(base[j], base[j]));

        // 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, sum;
        shuffle = _mm_shuffle_ps(tmp, tmp, _MM_SHUFFLE(1, 0, 3, 2));
        sum     = _mm_add_ps(tmp, shuffle) ;
        shuffle = _mm_shuffle_ps(sum, sum, _MM_SHUFFLE(2, 3, 0, 1));
        Real32 result  = _mm_cvtss_f32(_mm_add_ps(sum, shuffle));
        for (int j=IntAlignDown(MatNrCol(m), 4); j<MatNrCol(m) ; j++)
            result += result + *MatE(m, i, j);
        *MatE(res, i, 0) = result;
    }
    return res;
}

/* sum over the columns (result = size(1,C))
   As a special SIMD feature, the vector is enlarged to be a multiple of
   alignCount and the first few values of the vector are replicated. */
template<class ArrayT>
inline ArrayT*
SIMD_MatSquareAndSumAxis0(ArrayT* m, int alignCount)
{
    int cn = MatNrCol(m);
    ArrayT* res = MatCreate<ArrayT>(1, IntAlignUp(cn, alignCount));
    SetVal(res, res, 0.0);
    for (int i=0 ; i<MatNrRow(m) ; i++)
        for (int j=0 ; j<cn ; j++)
            *MatE(res, 0, j) = *MatE(res, 0, j) + (*MatE(m, i, j)) * (*MatE(m, i, j));
    for(int j = 0; j < IntAlignUp(cn, alignCount) - cn; j++)
    {
        *MatE(res, 0, cn + j) = *MatE(res, 0, j);
    }
    return res;
}

template<class ArrayT>
ArrayT*
MatPadZeros(ArrayT* data, int rowCount, int columnCount)
{
    ArrayT* res = MatCreate<ArrayT>(rowCount, columnCount);
    for(int i = 0; i < MatNrRow(data); i++)
    {
        int j = 0;
        for(; j < MatNrCol(data); j++)
        {
            *MatE(res, i, j) = *MatE(data, i, j);
        }
        for(; j < columnCount; j++);
        {
            *MatE(res, i, j) = 0.0;
        }
    }
    for(int i = MatNrRow(data); i < rowCount; i++)
        for(int j = 0; j < columnCount; j++)
            *MatE(res, i, j) = 0.0;
    return res;
}


// WARNING: function returns a matrix of size (M,N+i) where i can be 0/1/2/3
template<class ArrayT>
inline ArrayT*
MatNorm2DistSSE(ArrayT* aT, ArrayT* b)
{
    ILOG_VAR(Core.Matrix.MatNorm2DistSSE);
    if (MatNrCol(aT) != MatNrRow(b)) {
        ILOG_ERROR("MatNorm2DistSSE operands: dimensionality problem");
    }

    if(!POINTER_ALIGNED(MatE(aT, 0, 0)))
    {
        ILOG_ERROR("Pointer alignment error aT");
    }
    if(!POINTER_ALIGNED(MatE(b, 0, 0)))
    {
        ILOG_ERROR("Pointer alignment error b");
    }

    Timer timer;
    int alignCount = 4;
    if(sizeof(typename ArrayT::StorType) == 8)
        alignCount = 2;
    

    // aa = sum(multiply(a,a), axis=0)      # aa = sum(a.*a,1)
    ArrayT* aaT = SIMD_MatSquareAndSumAxis1(aT);

    ILOG_DEBUG("aaT: " << timer.SplitTimeStr());

    // bb = sum(multiply(b,b), axis=0)      # bb = sum(b.*b,1)
    // this operation ensures that later on the number of columns is always
    // a multiple of 2 inside bb
    ArrayT* bb = SIMD_MatSquareAndSumAxis0(b, alignCount);

    ILOG_DEBUG("bb: " << timer.SplitTimeStr());

    // ab = dot(transpose(a),b)             # ab = a'*b
    //ArrayT* aT = MatTranspose(a);
    //ILOG_DEBUG("transpose: " << timer.SplitTimeStr());

    ILOG_DEBUG("aT " << MatNrRow(aT) << " " << MatNrCol(aT));
    ILOG_DEBUG("b  " << MatNrRow(b) << " " << MatNrCol(b));

    ArrayT* paddedB = MatPadZeros(b, MatNrRow(b), IntAlignUp(MatNrCol(b), alignCount));
    ArrayT* ab = MatMul(aT, paddedB);
    delete paddedB;

    ILOG_DEBUG("matmul: " << timer.SplitTimeStr());

    // return sqrt(abs(transpose(repmat(aa,bb.shape[0],1)) + repmat(bb,aa.shape[0],1) - 2*ab))
    // #d = sqrt(abs(repmat(aa',[1 size(bb,2)]) + repmat(bb,[size(aa,2) 1]) - 2*ab));

    //ILOG_DEBUG("aaT " << MatNrRow(aaT) << " " << MatNrCol(aaT));
    //ILOG_DEBUG("bb  " << MatNrRow(bb) << " " << MatNrCol(bb));

    if(sizeof(typename ArrayT::StorType) == 8)
    {
        #pragma omp parallel for
        for (int i=0 ; i<MatNrRow(aaT) ; i++)
        {
            const __m128d tmp = _mm_set1_pd(*MatE(aaT, i, 0));
            const __m128d minusTwo = _mm_set1_pd(-2.0);
            const UInt64 mask = 0x7FFFFFFFFFFFFFFFL;
            const __m128d absMask = _mm_set1_pd(*(double*)(&mask));
            __m128d* baseBB = (__m128d*)(MatE(bb, 0, 0));
            __m128d* baseAB = (__m128d*)(MatE(ab, i, 0));
            const int SSELength = MatNrCol(bb) / 2;
            for(int j = 0; j < SSELength; j++)
            {
                __m128d intermediate = _mm_add_pd(_mm_add_pd(baseBB[j], tmp), _mm_mul_pd(minusTwo, baseAB[j]));
                baseAB[j] = _mm_sqrt_pd(_mm_and_pd(intermediate, absMask));
            }
        }
    }
    else
    {
        #pragma omp parallel for
        for (int i=0 ; i<MatNrRow(aaT) ; i++)
        {
            const __m128 tmp = _mm_set1_ps(*MatE(aaT, i, 0));
            const __m128 minusTwo = _mm_set1_ps(-2.0);
            const UInt32 mask = 0x7FFFFFFF;
            const __m128 absMask = _mm_set1_ps(*(Real32*)(&mask));
            __m128* baseBB = (__m128*)(MatE(bb, 0, 0));
            __m128* baseAB = (__m128*)(MatE(ab, i, 0));
            const int SSELength = MatNrCol(bb) / 4;
            for(int j = 0; j < SSELength; j++)
            {
                __m128 intermediate = _mm_add_ps(_mm_add_ps(baseBB[j], tmp), _mm_mul_ps(minusTwo, baseAB[j]));
                baseAB[j] = _mm_sqrt_ps(_mm_and_ps(intermediate, absMask));
            }
        }
    }
    ILOG_DEBUG("sqrt(abs(aa + bb - 2ab)): " << timer.SplitTimeStr());
    delete bb;
    delete aaT;

    ILOG_DEBUG("delete: " << timer.SplitTimeStr());
//WriteRaw(a, "matrix_a.raw", &Util::Database::GetInstance(), 1);
//WriteRaw(b, "matrix_b.raw", &Util::Database::GetInstance(), 1);
//WriteRaw(repAAT, "matrix_c.raw", &Util::Database::GetInstance(), 1);
    //ILOG_INFO("cpu-matnorm2dist-total: " << timer.SplitTime());
    ILOG_DEBUG(timer.SplitTime() << " (cpu-matnorm2dist-total)");
    return ab;
}

} // namespace Matrix
} // namespace Core
} // namespace Impala

#endif

---