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

---

MatMul.h

Go to the documentation of this file.

#ifndef Impala_Core_Matrix_MatMul_h
#define Impala_Core_Matrix_MatMul_h

#include "Core/Matrix/Mat.h"
#include "Core/Matrix/MatFunc.h"
#include "Core/Matrix/MatTranspose.h"


#ifdef GSL_CBLAS
#include <gsl/gsl_cblas.h>
#endif

#ifdef MKL_CBLAS
#include <mkl_cblas.h>
#endif

#ifdef ATLAS_CBLAS
extern "C"{
    #include <cblas.h>
}
#endif

#ifdef ACML_CBLAS
//ACML defines the signatures of CBlas functions differently
//there is no cblas_ prefix to the BLAS functions
//Also there is no option for choosing ColMajor or RowMajor
//It defaults to rowmajor
    #include <acml.h>

    #define CblasTrans 'T'
    #define CblasNoTrans 'N'

    #define SGEMM(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12,X13) \
            sgemm(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12,X13)
    #define DGEMM(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12,X13) \
            dgemm(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12,X13)

#else
//For other CBlas interfaces that suit the cblas.h definitions,
//we only need to add a cblas_ prefix, and mention it's RowMajor
//
    #define SGEMM(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12,X13) \
            cblas_sgemm(CblasRowMajor,X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12,X13)

    #define DGEMM(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12,X13) \
            cblas_dgemm(CblasRowMajor,X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12,X13)
    #endif

//If no external library is used, use the reference implementation
#ifndef GSL_CBLAS
#ifndef MKL_CBLAS
#ifndef ATLAS_CBLAS
#ifndef ACML_CBLAS
    #define BASE_CBLAS
#endif
#endif
#endif
#endif


namespace Impala
{
namespace Core
{
namespace Matrix
{

template<class ArrayT>
inline ArrayT*
MatMul(ArrayT* m, ArrayT* a)
{
    ArrayT* res;
    ILOG_VAR(Core.Matrix.MatMul);

    if (MatNrCol(m) != MatNrRow(a))
    {
        ILOG_ERROR("nonconformant MatMul operands.");
        return 0;
    }
    res = MatCreate<ArrayT>(MatNrRow(m), MatNrCol(a));

#ifdef BASE_CBLAS
    Mat* aa = MatTranspose(a);
    for (int i=0 ; i<MatNrRow(m) ; i++)
    {
        double* rowM = MatE(m, i, 0);
        for (int j=0 ; j<MatNrCol(a) ; j++)
        {
            double* rowAA = MatE(aa, j, 0);
            MatStorType sum = 0;
            for (int k=0 ; k<MatNrCol(m) ; k++)
            {
                sum += rowM[k] * rowAA[k];
            }
            *MatE(res, i, j) = sum;
        }
    }
    delete aa;
#else
    DGEMM(CblasNoTrans, CblasNoTrans,
          MatNrRow(m), MatNrCol(a), MatNrRow(a),
          1.0, m->PB(), MatNrCol(m),
          a->PB(), MatNrCol(a),
          0.0, res->PB(), MatNrCol(res));
#endif

    return res;
}

inline Mat32*
MatMul(Mat32* m, Mat32* a)
{
    Mat32* res;
    ILOG_VAR(Core.Matrix.MatMul);

    if (MatNrCol(m) != MatNrRow(a))
    {
        ILOG_ERROR("nonconformant MatMul operands.");
        return 0;
    }
    res = MatCreate<Mat32>(MatNrRow(m), MatNrCol(a));

#ifdef BASE_CBLAS
    Mat32* aa = MatTranspose(a);
    for (int i=0 ; i<MatNrRow(m) ; i++)
    {
        Real32* rowM = MatE(m, i, 0);
        for (int j=0 ; j<MatNrCol(a) ; j++)
        {
            Real32* rowAA = MatE(aa, j, 0);
            Real32 sum = 0;
            for (int k=0 ; k<MatNrCol(m) ; k++)
            {
                sum += rowM[k] * rowAA[k];
            }
            *MatE(res, i, j) = sum;
        }
    }
    delete aa;
#else
    SGEMM(CblasNoTrans, CblasNoTrans,
          MatNrRow(m), MatNrCol(a), MatNrRow(a),
          1.0, m->PB(), MatNrCol(m),
          a->PB(), MatNrCol(a),
          0.0, res->PB(), MatNrCol(res));
#endif

    return res;
}

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

#endif

---