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_SSE2_MATMUL_H
#define _EXMAT_SIMD_SEE2_MATMUL_H


#if HAVE_PARTIAL_SPECIALIZATION

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

namespace exmat {
namespace SIMD {
namespace PNS {


#define _mm_load_pd_1(input, count) (_mm_shuffle_pd(input, input, _MM_SHUFFLE2(count, count)))
#define _mm_load_sd_1(input, count) (_mm_load_sd(((double*)&input) + count))


template<int iRow1, int iCol, int iRow2, int i = iRow1> struct SSE2MatMulRowEvaluator;
template<int iRow1, int iCol, int iRow2, int i = iRow1> struct SSE2MatMulRowEvaluatorS;
template<int iRow1, int iCol, int iRow2, int i = iCol>  struct SSE2MatMulColEvaluator;
template<int iRow, int iCol> struct SSE2MatMulRowElmEvaluatorS;
template<int iRow, int iCol> struct SSE2MatMulRowElmEvaluator;

template<int iRow1, int iCol, int iRow2, int i> // column wise
struct SSE2MatMulColEvaluator
{
    EXMAT_INLINE2 static void
    eval(const __m128d* pmatRowStart, const __m128d* pmat, __m128d* pmatDest)
    {
        SSE2MatMulRowEvaluator<iRow1, (iCol + 1) / 2, iRow2>::eval(pmatRowStart, pmat, pmatDest);
        SSE2MatMulColEvaluator<iRow1, iCol, iRow2, i - 2>::eval(pmatRowStart, pmat + 1, pmatDest + 1);
    };
};

template<int iRow1, int iCol, int iRow2> // column wise
struct SSE2MatMulColEvaluator<iRow1, iCol, iRow2, 2>
{
    EXMAT_INLINE2 static void
    eval(const __m128d* pmatRowStart, const __m128d* pmat, __m128d* pmatDest)
    {
        SSE2MatMulRowEvaluator<iRow1, (iCol + 1) / 2, iRow2>::eval(pmatRowStart, pmat, pmatDest);
    };
};

template<int iRow1, int iCol, int iRow2> // column wise
struct SSE2MatMulColEvaluator<iRow1, iCol, iRow2, 1>
{
    EXMAT_INLINE2 static void
    eval(const __m128d* pmatRowStart, const __m128d* pmat, __m128d* pmatDest)
    {
        SSE2MatMulRowEvaluatorS<iRow1, (iCol + 1) / 2, iRow2>::eval(pmatRowStart, pmat, pmatDest);
    };
};

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

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

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

template<int iRow1, int iCol, int iRow2> // row wise
struct SSE2MatMulRowEvaluator<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 __m128d* pmatRowStart, const __m128d* pmat, __m128d* pmatDest)
    {
        pmatDest[0] = SSE2MatMulRowElmEvaluator<iRow2, iCol>::eval(pmatRowStart, pmat);
    };
};

template<int iRow, int iCol>
struct SSE2MatMulRowElmEvaluatorS //perform Row * Col
{
    EXMAT_INLINE2 static __m128d eval(const __m128d* pmatRow, const __m128d* pmat)
    {
        return _mm_add_ps(_mm_add_sd(_mm_mul_sd(_mm_load_sd_1(*pmatRow, 0), pmat[0]), _mm_mul_sd(_mm_load_sd_1(*pmatRow, 1), pmat[1 * iCol])),
            SSE2MatMulRowElmEvaluatorS<iRow - 2, iCol>::eval(pmatRow + 1, pmat + 2 * iCol));
    };
};

template<int iCol>
struct SSE2MatMulRowElmEvaluatorS<2, iCol> //perform Row * Col
{
    EXMAT_INLINE2 static __m128d eval(const __m128d* pmatRow, const __m128d* pmat)
    {
        return _mm_add_sd(_mm_mul_sd(_mm_load_sd_1(*pmatRow, 0), pmat[0]), _mm_mul_sd(_mm_load_sd_1(*pmatRow, 1), pmat[1 * iCol]));
    };
};

template<int iCol>
struct SSE2MatMulRowElmEvaluatorS<1, iCol> //perform Row * Col
{
    EXMAT_INLINE2 static __m128d eval(const __m128d* pmatRow, const __m128d* pmat)
    {
        return _mm_mul_sd(_mm_load_sd_1(*pmatRow, 0), pmat[0]);
    };
};

template<int iRow, int iCol>
struct SSE2MatMulRowElmEvaluator //perform Row * Col
{
    EXMAT_INLINE2 static __m128d eval(const __m128d* pmatRow, const __m128d* pmat)
    {
        return _mm_add_pd(_mm_add_pd(_mm_mul_pd(_mm_load_pd_1(*pmatRow, 0), pmat[0]), _mm_mul_pd(_mm_load_pd_1(*pmatRow, 1), pmat[1 * iCol])),
            SSE2MatMulRowElmEvaluator<iRow - 2, iCol>::eval(pmatRow + 1, pmat + 2 * iCol));
    };
};

template<int iCol>
struct SSE2MatMulRowElmEvaluator<2, iCol> //perform Row * Col
{
    EXMAT_INLINE2 static __m128d eval(const __m128d* pmatRow, const __m128d* pmat)
    {
        return _mm_add_pd(_mm_mul_pd(_mm_load_pd_1(*pmatRow, 0), pmat[0]), _mm_mul_pd(_mm_load_pd_1(*pmatRow, 1), pmat[1 * iCol]));
    };
};

template<int iCol>
struct SSE2MatMulRowElmEvaluator<1, iCol> //perform Row * Col
{
    EXMAT_INLINE2 static __m128d eval(const __m128d* pmatRow, const __m128d* pmat)
    {
        return _mm_mul_pd(_mm_load_pd_1(*pmatRow, 0), pmat[0]);
    };
};

#undef _mm_load_pd_1
#undef _mm_load_sd_1

// For RMat is a column vector, we pretend the RHS col vec as a LHS row vec
template<class DMat, class LMat, class RMat> EXMAT_INLINE2 static void
SIMDMatMul_inner_double(DMat& d, const LMat& l, const RMat& r, Int2Type<true> isRHSCVec) {
    SSE2MatMulColEvaluator<1,RMat::ROWS,LMat::COLS>::
        eval(
            (const __m128d*)(&r.unitData[0]),
            (const __m128d*)(&l.unitData[0]),
            (__m128d*)(&d.unitData[0])
        );
}
template<class DMat, class LMat, class RMat> EXMAT_INLINE2 static void
SIMDMatMul_inner_double(DMat& d, const LMat& l, const RMat& r, Int2Type<false> isRHSCVec) {
    SSE2MatMulColEvaluator<LMat::ROWS,LMat::COLS,RMat::COLS>::
        eval(
            (const __m128d*)(&l.unitData[0]),
            (const __m128d*)(&r.unitData[0]),
            (__m128d*)(&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<double>) {
    using namespace PNS;
    enum {
        isRHSCVec = (RMat::COLS == 1)
    };
    SIMDMatMul_inner_double(d, l, r, Int2Type<isRHSCVec>());
}


};
};

#endif  // HAVE_PARTIAL_SPECIALIZATION


#endif  // _EXMAT_SIMD_SSE22_MATMUL_H

---