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

---

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



//
//Uncomment this to see the performance in usec
//#define ENABLE_PERFORMANCE_TEST
//
#ifdef ENABLE_PERFORMANCE_TEST
    #include "sys/time.h"
#endif


//  Uncomment this part to test the numerical differences 
//  between the ref. implementation
//#define ENABLE_PRECISION_TEST
//
#ifdef ENABLE_PRECISION_TEST
    #include "Core/Array/Sub.h"
    #include "Core/Array/PixSum.h"
    #include "Core/Array/PrintData.h"
#endif

namespace Impala
{
namespace Core
{
namespace Matrix
{

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

#ifdef ENABLE_PERFORMANCE_TIMING
    timeval start,end,dur;
    gettimeofday(&start,NULL);
#endif

    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));
#ifdef ENABLE_PRECISION_TEST
    ArrayT* res2;
    res2 = MatCreate<ArrayT>(MatNrRow(m), MatNrCol(a));
    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(res2, i, j) = sum;
        }
    }
    delete aa;
    ArrayT* r=0;
    Sub(r,res,res2);
    Impala::Real64 diff=PixSum(r);
    if(diff!=0){
        ILOG_WARN("The Arrays are different!");
        PrintData(r);
        
    }
    else{
        ILOG_INFO("Arrays are same!");
    }
#endif

#endif

#ifdef ENABLE_PERFORMANCE_TIMING
    gettimeofday(&end,NULL);
    long long int sec = end.tv_sec - start.tv_sec;
    long long int usec= end.tv_usec - start.tv_usec;

    ILOG_DEBUG("MatMul : "<<sec*1000000+usec<<" usec");
#endif

    return res;

}

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

#endif

---