Ausblick

#include <cstddef>
#include <cstdio>
#include "bench.h"
#include "ulmblas.h"
#include "gemm_refcolmajor.h"
#include "gemm_blocked.h"
#include "gematrix.h"

#ifndef COLMAJOR
#define COLMAJOR 1
#endif

#ifndef MAXDIM_M
#define MAXDIM_M    7000
#endif

#ifndef MAXDIM_N
#define MAXDIM_N    7000
#endif

#ifndef MAXDIM_K
#define MAXDIM_K    7000
#endif

#ifndef MIN_M
#define MIN_M   100
#endif

#ifndef MIN_N
#define MIN_N   100
#endif

#ifndef MIN_K
#define MIN_K   100
#endif

#ifndef MAX_M
#define MAX_M   7000
#endif

#ifndef MAX_N
#define MAX_N   7000
#endif

#ifndef MAX_K
#define MAX_K   7000
#endif

#ifndef INC_M
#define INC_M   100
#endif

#ifndef INC_N
#define INC_N   100
#endif

#ifndef INC_K
#define INC_K   100
#endif

#ifndef ALPHA
#define ALPHA   1.5
#endif

#ifndef BETA
#define BETA    1.5
#endif

int
main()
{
    using namespace matvec;

    typedef double      T;

    StorageOrder order = (COLMAJOR) ? StorageOrder::ColMajor
                                    : StorageOrder::RowMajor;

    GeMatrix<T> A_(MAXDIM_M, MAXDIM_K, order);
    GeMatrix<T> B_(MAXDIM_K, MAXDIM_N, order);
    GeMatrix<T> C1_(MAXDIM_M, MAXDIM_N, order);
    GeMatrix<T> C2_(MAXDIM_M, MAXDIM_N, order);

    initGeMatrix(A_);
    initGeMatrix(B_);
    initGeMatrix(C1_);
    initGeMatrix(C2_);

    gecopy(C1_, C2_);

    // Header-Zeile fuer die Ausgabe
    printf("%5s %5s %5s ", "m", "n", "k");
    printf("%20s %9s", "refColMajor: t", "MFLOPS");
    printf("%20s %9s %9s", "blocked GEMM: t", "MFLOPS", "diff");
    printf("\n");

    bench::WallTime<double> wallTime;

    for (long m = MIN_M, n = MIN_N, k = MIN_K;
         m <=MAX_M && n <= MAX_N && k <= MAX_K;
         m += INC_M, n += INC_N, k += INC_K)
    {
        double t, diff;

        GeMatrixView<T> A  = A_(0,0,m,k);
        GeMatrixView<T> B  = B_(0,0,k,n);
        GeMatrixView<T> C1 = C1_(0,0,m,n);
        GeMatrixView<T> C2 = C2_(0,0,m,n);

        printf("%5ld %5ld %5ld ", m, n, k);

        wallTime.tic();
        gemm_ref(ALPHA, A, B, BETA, C1);
        t = wallTime.toc();
        printf("%20.4lf %9.2lf", t, 2.*m/1000*n/1000*k/t);

        wallTime.tic();
        gemm_blk(ALPHA, A, B, BETA, C2);
        t = wallTime.toc();
        diff = asumDiffGeMatrix(C1, C2);
        printf("%20.4lf %9.2lf %9.1e", t, 2.*m/1000*n/1000*k/t, diff);
        printf("\n");
    }
}

#ifndef HPC_GEMATRIX_H
#define HPC_GEMATRIX_H 1

#include <cstddef>
#include <cassert>
#include "bench.h"
#include "ulmblas.h"
#include "gemm_refcolmajor.h"
#include "gemm_blocked.h"

namespace matvec {

enum class StorageOrder {
    ColMajor,
    RowMajor
};

template <typename T, typename I>
    struct GeMatrixView;

template <typename T, typename I=std::size_t>
struct GeMatrix
{
    typedef T                       ElementType;
    typedef I                       Index;
    typedef GeMatrix<T,Index>       NoView;
    typedef GeMatrixView<T,Index>   View;

    GeMatrix(Index m, Index n, StorageOrder order=StorageOrder::ColMajor)
        : m(m), n(n),
          incRow(order==StorageOrder::ColMajor ? 1: n),
          incCol(order==StorageOrder::RowMajor ? 1: m),
          data(new T[m*n])
    {
    }

    ~GeMatrix()
    {
        delete[] data;
    }

    const ElementType &
    operator()(Index i, Index j) const
    {
        assert(i<m && j<n);
        return data[i*incRow + j*incCol];
    }

    ElementType &
    operator()(Index i, Index j)
    {
        assert(i<m && j<n);
        return data[i*incRow + j*incCol];
    }

    View
    operator()(Index i, Index j, Index numRows, Index numCols)
    {
        assert(i+numRows<=m);
        assert(j+numCols<=n);
        return View(numRows, numCols, &(operator()(i,j)), incRow, incCol);
    }

    const Index     m, n, incRow, incCol;
    ElementType*    data;
};

template <typename T, typename I=std::size_t>
struct GeMatrixView
{
    typedef T                       ElementType;
    typedef I                       Index;
    typedef GeMatrix<T,Index>       NoView;
    typedef GeMatrixView<T,Index>   View;

    GeMatrixView(Index m, Index n, T *data, Index incRow, Index incCol)
        : m(m), n(n), incRow(incRow), incCol(incCol), data(data)
    {
    }

    GeMatrixView(const GeMatrixView &rhs)
        : m(rhs.m), n(rhs.n), incRow(rhs.incRow), incCol(rhs.incCol),
          data(rhs.data)
    {
    }

    const ElementType &
    operator()(Index i, Index j) const
    {
        assert(i<m && j<n);
        return data[i*incRow + j*incCol];
    }

    ElementType &
    operator()(Index i, Index j)
    {
        assert(i<m && j<n);
        return data[i*incRow + j*incCol];
    }

    View
    operator()(Index i, Index j, Index numRows, Index numCols)
    {
        assert(i+numRows<=m);
        assert(j+numCols<=n);
        return View(numRows, numCols, &(operator()(i,j)), incRow, incCol);
    }

    const Index     m, n, incRow, incCol;
    ElementType*    data;
};

//
//  Interface for bench
//
template <typename GeMatrix>
void
initGeMatrix(GeMatrix &A)
{
    bench::initGeMatrix(A.m, A.n, A.data, A.incRow, A.incCol);
}

template <typename GeMatrix>
typename GeMatrix::ElementType
asumDiffGeMatrix(const GeMatrix &A, const GeMatrix &B)
{
    return bench::asumDiffGeMatrix(A.m, A.n,
                                   A.data, A.incRow, A.incCol,
                                   B.data, B.incRow, B.incCol);
}

//
// Interface for ulmBLAS
//
template <typename GeMatrix>
void
gecopy(const GeMatrix &A, GeMatrix &B)
{
    assert(A.m==B.m && A.n==B.n);
    ulmBLAS::gecopy(A.m, A.n,
                    A.data, A.incRow, A.incCol,
                    B.data, B.incRow, B.incCol);

}

//
// Interface for refColMajor
//
template <typename GeMatrix>
void
gemm_ref(typename GeMatrix::ElementType alpha,
         const GeMatrix                 &A,
         const GeMatrix                 &B,
         typename GeMatrix::ElementType beta,
         GeMatrix                       &C)
{
    assert(A.n==B.m);
    assert(A.m==C.m);
    assert(A.n==C.n);

    refColMajor::gemm(C.m, C.n, A.n,
                      alpha,
                      A.data, A.incRow, A.incCol,
                      B.data, B.incRow, B.incCol,
                      beta,
                      C.data, C.incRow, C.incCol);
}

//
// Interface for blocked
//
template <typename GeMatrix>
void
gemm_blk(typename GeMatrix::ElementType alpha,
         const GeMatrix                 &A,
         const GeMatrix                 &B,
         typename GeMatrix::ElementType beta,
         GeMatrix                       &C)
{
    assert(A.n==B.m);
    assert(A.m==C.m);
    assert(A.n==C.n);

    blocked::gemm(C.m, C.n, A.n,
                  alpha,
                  A.data, A.incRow, A.incCol,
                  B.data, B.incRow, B.incCol,
                  beta,
                  C.data, C.incRow, C.incCol);
}

} // namespace matvec

#endif // HPC_GEMATRIX_H