MatMul.h

/*
 * Expression Template Matrix Library
 *
 * Copyright (C) 2004 - 2006 Ricky Lung <mtlung@users.sourceforge.net>
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */


#ifndef _EXMAT_SIMD_SSE_MATMUL_H
#define _EXMAT_SIMD_SSE_MATMUL_H


#if HAVE_PARTIAL_SPECIALIZATION

#include "../../PlatformSpec.h"
#include "../../Metaprogramming.h"

namespace exmat {
namespace SIMD {
namespace PNS {


#define _mm_load_ps_1(input, count) (_mm_shuffle_ps(input, input, _MM_SHUFFLE(count, count, count, count)))
#define _mm_load_ss_1(input, count) (_mm_load_ss(((float*)&input) + count))


template<int iRow1, int iCol, int iRow2, int i = iRow1> struct SSEMatMulRowEvaluator;
template<int iRow1, int iCol, int iRow2, int i = iRow1> struct SSEMatMulRowEvaluatorS;
template<int iRow1, int iCol, int iRow2, int i = iCol>  struct SSEMatMulColEvaluator;
template<int iRow, int iCol> struct SSEMatMulRowElmEvaluatorS;
template<int iRow, int iCol> struct SSEMatMulRowElmEvaluator;

template<int iRow1, int iCol, int iRow2, int i> // column wise
struct SSEMatMulColEvaluator
{
    EXMAT_INLINE2 static void
    eval(const __m128* pmatRowStart, const __m128* pmat, __m128* pmatDest)
    {
        SSEMatMulRowEvaluator<iRow1, (iCol + 3) / 4, iRow2>::eval(pmatRowStart, pmat, pmatDest);
        SSEMatMulColEvaluator<iRow1, iCol, iRow2, i - 4>::eval(pmatRowStart, pmat + 1, pmatDest + 1);
    };
};

template<int iRow1, int iCol, int iRow2> // column wise
struct SSEMatMulColEvaluator<iRow1, iCol, iRow2, 4>
{
    EXMAT_INLINE2 static void
    eval(const __m128* pmatRowStart, const __m128* pmat, __m128* pmatDest)
    {
        SSEMatMulRowEvaluator<iRow1, (iCol + 3) / 4, iRow2>::eval(pmatRowStart, pmat, pmatDest);
    };
};

template<int iRow1, int iCol, int iRow2> // column wise
struct SSEMatMulColEvaluator<iRow1, iCol, iRow2, 3>
{
    EXMAT_INLINE2 static void
    eval(const __m128* pmatRowStart, const __m128* pmat, __m128* pmatDest)
    {
        SSEMatMulRowEvaluator<iRow1, (iCol + 3) / 4, iRow2>::eval(pmatRowStart, pmat, pmatDest);
    };
};


template<int iRow1, int iCol, int iRow2> // column wise
struct SSEMatMulColEvaluator<iRow1, iCol, iRow2, 2>
{
    EXMAT_INLINE2 static void
    eval(const __m128* pmatRowStart, const __m128* pmat, __m128* pmatDest)
    {
        SSEMatMulRowEvaluator<iRow1, (iCol + 3) / 4, iRow2>::eval(pmatRowStart, pmat, pmatDest);
    };
};

template<int iRow1, int iCol, int iRow2> // column wise
struct SSEMatMulColEvaluator<iRow1, iCol, iRow2, 1>
{
    EXMAT_INLINE2 static void
    eval(const __m128* pmatRowStart, const __m128* pmat, __m128* pmatDest)
    {
        SSEMatMulRowEvaluatorS<iRow1, (iCol + 3) / 4, iRow2>::eval(pmatRowStart, pmat, pmatDest);
    };
};

template<int iRow1, int iCol, int iRow2, int i> // row wise
struct SSEMatMulRowEvaluatorS
{
    EXMAT_INLINE2 static void
    eval(const __m128* pmatRowStart, const __m128* pmat, __m128* pmatDest)
    {
        pmatDest[0] = SSEMatMulRowElmEvaluatorS<iRow2, iCol>::eval(pmatRowStart, pmat);
        SSEMatMulRowEvaluatorS<iRow1, iCol, iRow2, i - 1>::eval(pmatRowStart + 1 * iCol, pmat, pmatDest + 1 * iCol);
    };
};

template<int iRow1, int iCol, int iRow2> // row wise
struct SSEMatMulRowEvaluatorS<iRow1, iCol, iRow2, 1>
{
    EXMAT_INLINE2 static void
    eval(const __m128* pmatRowStart, const __m128* pmat, __m128* pmatDest)
    {
        pmatDest[0] = SSEMatMulRowElmEvaluatorS<iRow2, iCol>::eval(pmatRowStart, pmat);
    };
};

template<int iRow1, int iCol, int iRow2, int i> // row wise
struct SSEMatMulRowEvaluator
{
    EXMAT_INLINE2 static void
    eval(const __m128* pmatRowStart, const __m128* pmat, __m128* pmatDest)
    {
//      _mm_prefetch( (char*)(pmatRowStart+2), 1 );
        pmatDest[0] = SSEMatMulRowElmEvaluator<iRow2, iCol>::eval(pmatRowStart, pmat);
        SSEMatMulRowEvaluator<iRow1, iCol, iRow2, i - 1>::eval(pmatRowStart + 1 * iCol, pmat, pmatDest + 1 * iCol);
    };
};

template<int iRow1, int iCol, int iRow2> // row wise
struct SSEMatMulRowEvaluator<iRow1, iCol, iRow2, 1>//this iRow is the iRow of first matrix
                                        //where the iRow of ElmEva, is the iRow of second matrix. a fix is needed
{
    EXMAT_INLINE2 static void
    eval(const __m128* pmatRowStart, const __m128* pmat, __m128* pmatDest)
    {
        pmatDest[0] = SSEMatMulRowElmEvaluator<iRow2, iCol>::eval(pmatRowStart, pmat);
    };
};

template<int iRow, int iCol>
struct SSEMatMulRowElmEvaluatorS //perform Row * Col
{
    EXMAT_INLINE2 static __m128 eval(const __m128* pmatRow, const __m128* pmat)
    {
        return _mm_add_ss(_mm_add_ss(_mm_add_ss(_mm_add_ss(_mm_mul_ss(_mm_load_ss_1(*pmatRow, 0), pmat[0]), _mm_mul_ss(_mm_load_ss_1(*pmatRow, 1), pmat[1 * iCol])), _mm_mul_ss(_mm_load_ss_1(*pmatRow, 2), pmat[2 * iCol])), _mm_mul_ss(_mm_load_ss_1(*pmatRow, 3), pmat[3 * iCol])), 
            SSEMatMulRowElmEvaluatorS<iRow - 4, iCol>::eval(pmatRow + 1, pmat + 4 * iCol));
    };
};

template<int iCol>
struct SSEMatMulRowElmEvaluatorS<4, iCol> //perform Row * Col
{
    EXMAT_INLINE2 static __m128 eval(const __m128* pmatRow, const __m128* pmat)
    {
        return _mm_add_ss(_mm_add_ss(_mm_add_ss(_mm_mul_ss(_mm_load_ss_1(*pmatRow, 0), pmat[0]), _mm_mul_ss(_mm_load_ss_1(*pmatRow, 1), pmat[1 * iCol])), _mm_mul_ss(_mm_load_ss_1(*pmatRow, 2), pmat[2 * iCol])), _mm_mul_ss(_mm_load_ss_1(*pmatRow, 3), pmat[3 * iCol]));
    };
};

template<int iCol>
struct SSEMatMulRowElmEvaluatorS<3, iCol> //perform Row * Col
{
    EXMAT_INLINE2 static __m128 eval(const __m128* pmatRow, const __m128* pmat)
    {
        return _mm_add_ss(_mm_add_ss(_mm_mul_ss(_mm_load_ss_1(*pmatRow, 0), pmat[0]), _mm_mul_ss(_mm_load_ss_1(*pmatRow, 1), pmat[1 * iCol])), _mm_mul_ss(_mm_load_ss_1(*pmatRow, 2), pmat[2 * iCol]));
    };
};

template<int iCol>
struct SSEMatMulRowElmEvaluatorS<2, iCol> //perform Row * Col
{
    EXMAT_INLINE2 static __m128 eval(const __m128* pmatRow, const __m128* pmat)
    {
        return _mm_add_ss(_mm_mul_ss(_mm_load_ss_1(*pmatRow, 0), pmat[0]), _mm_mul_ss(_mm_load_ss_1(*pmatRow, 1), pmat[1 * iCol]));
    };
};

template<int iCol>
struct SSEMatMulRowElmEvaluatorS<1, iCol> //perform Row * Col
{
    EXMAT_INLINE2 static __m128 eval(const __m128* pmatRow, const __m128* pmat)
    {
        return _mm_mul_ss(_mm_load_ss_1(*pmatRow, 0), pmat[0]);
    };
};

template<int iRow, int iCol>
struct SSEMatMulRowElmEvaluator //perform Row * Col
{
    EXMAT_INLINE2 static __m128 eval(const __m128* pmatRow, const __m128* pmat)
    {
        return _mm_add_ps(_mm_add_ps(_mm_add_ps(_mm_add_ps(_mm_mul_ps(_mm_load_ps_1(*pmatRow, 0), pmat[0]), _mm_mul_ps(_mm_load_ps_1(*pmatRow, 1), pmat[1 * iCol])), _mm_mul_ps(_mm_load_ps_1(*pmatRow, 2), pmat[2 * iCol])), _mm_mul_ps(_mm_load_ps_1(*pmatRow, 3), pmat[3 * iCol])), 
            SSEMatMulRowElmEvaluator<iRow - 4, iCol>::eval(pmatRow + 1, pmat + 4 * iCol));
    };
};

template<int iCol>
struct SSEMatMulRowElmEvaluator<4, iCol> //perform Row * Col
{
    EXMAT_INLINE2 static __m128 eval(const __m128* pmatRow, const __m128* pmat)
    {
//      _mm_prefetch( (char*)(pmatRow+1), 1 );
        return _mm_add_ps(_mm_add_ps(_mm_add_ps(_mm_mul_ps(_mm_load_ps_1(*pmatRow, 0), pmat[0]), _mm_mul_ps(_mm_load_ps_1(*pmatRow, 1), pmat[1 * iCol])), _mm_mul_ps(_mm_load_ps_1(*pmatRow, 2), pmat[2 * iCol])), _mm_mul_ps(_mm_load_ps_1(*pmatRow, 3), pmat[3 * iCol]));
    };
};

template<int iCol>
struct SSEMatMulRowElmEvaluator<3, iCol> //perform Row * Col
{
    EXMAT_INLINE2 static __m128 eval(const __m128* pmatRow, const __m128* pmat)
    {
        return _mm_add_ps(_mm_add_ps(_mm_mul_ps(_mm_load_ps_1(*pmatRow, 0), pmat[0]), _mm_mul_ps(_mm_load_ps_1(*pmatRow, 1), pmat[1 * iCol])), _mm_mul_ps(_mm_load_ps_1(*pmatRow, 2), pmat[2 * iCol]));
    };
};

template<int iCol>
struct SSEMatMulRowElmEvaluator<2, iCol> //perform Row * Col
{
    EXMAT_INLINE2 static __m128 eval(const __m128* pmatRow, const __m128* pmat)
    {
        return _mm_add_ps(_mm_mul_ps(_mm_load_ps_1(*pmatRow, 0), pmat[0]), _mm_mul_ps(_mm_load_ps_1(*pmatRow, 1), pmat[1 * iCol]));
    };
};

template<int iCol>
struct SSEMatMulRowElmEvaluator<1, iCol> //perform Row * Col
{
    EXMAT_INLINE2 static __m128 eval(const __m128* pmatRow, const __m128* pmat)
    {
        return _mm_mul_ps(_mm_load_ps_1(*pmatRow, 0), pmat[0]);
    };
};

#undef _mm_load_ps_1
#undef _mm_load_ss_1

// For RMat is a column vector, we pretended the RHS col vec as a LHS row vec
template<class DMat, class LMat, class RMat> EXMAT_INLINE2 static void
SIMDMatMul_inner_float(DMat& d, const LMat& l, const RMat& r, Int2Type<true> isRHSCVec) {
    SSEMatMulColEvaluator<1,RMat::ROWS,LMat::COLS>::
        eval(
            (const __m128*)(&r.unitData[0]),
            (const __m128*)(&l.unitData[0]),
            (__m128*)(&d.unitData[0])
        );
}
template<class DMat, class LMat, class RMat> EXMAT_INLINE2 static void
SIMDMatMul_inner_float(DMat& d, const LMat& l, const RMat& r, Int2Type<false> isRHSCVec) {
    SSEMatMulColEvaluator<LMat::ROWS,RMat::COLS,RMat::ROWS>::
        eval(
            (const __m128*)(&l.unitData[0]),
            (const __m128*)(&r.unitData[0]),
            (__m128*)(&d.unitData[0])
        );
}

}   // namespace exmat


template<class DMat, class LMat, class RMat> EXMAT_INLINE0 static void
SIMDMatMul(DMat& d, const LMat& l, const RMat& r, Type2Type<float>) {
    using namespace PNS;
    enum {
        isRHSCVec = (RMat::COLS == 1)
    };
    _mm_prefetch( (char*)(&l)+8, 1 );
    SIMDMatMul_inner_float(d, l, r, Int2Type<isRHSCVec>());
}


};  // namespace SIMD
}   // namespace exmat

#endif  // HAVE_PARTIAL_SPECIALIZATION


#endif  // _EXMAT_SIMD_SSE_MATMUL_H

---