1
       2
       3
       4
       5
       6
       7
       8
       9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      20
      21
      22
      23
      24
      25
      26
      27
      28
      29
      30
      31
      32
      33
      34
      35
      36
      37
      38
      39
      40
      41
      42
      43
      44
      45
      46
      47
      48
      49
      50
      51
      52
      53
      54
      55
      56
      57
      58
      59
      60
      61
      62
      63
      64
      65
      66
      67
      68
      69
      70
      71
      72
      73
      74
      75
      76
      77
      78
      79
      80
      81
      82
      83
      84
      85
      86
      87
      88
      89
      90
      91
      92
      93
      94
      95
      96
      97
      98
      99
     100
     101
     102
     103
     104
     105
     106
     107
     108
     109
     110
     111
     112
     113
     114
     115
     116
     117
     118
     119
     120
     121
     122
     123
     124
     125
     126
     127
     128
     129
     130
     131
     132
     133
     134
     135
     136
     137
     138
     139
     140
     141
     142
     143
     144
     145
     146
     147
     148
     149
     150
     151
     152
     153
     154
     155
     156
     157
     158
     159
     160
     161
     162
     163
     164
     165
#include <cassert>
#include <cstddef>
#include <cstdio>
#include "bench.h"
#include "ulmblas.h"
#include "gemm_refcolmajor.h"
#include "gemm_blocked.h"
#include "gematrix.h"
//#include <mkl.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   256
#endif

#ifndef MIN_N
#define MIN_N   256
#endif

#ifndef MIN_K
#define MIN_K   256
#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   256
#endif

#ifndef INC_N
#define INC_N   256
#endif

#ifndef INC_K
#define INC_K   256
#endif

#ifndef ALPHA
#define ALPHA   1.5
#endif

#ifndef BETA
#define BETA    1.5
#endif

/*
namespace mkl {

void
gemm(MKL_INT m, MKL_INT n, MKL_INT k,
     double alpha,
     const double *A, MKL_INT incRowA, MKL_INT incColA,
     const double *B, MKL_INT incRowB, MKL_INT incColB,
     double beta,
     double *C, MKL_INT incRowC, MKL_INT incColC)
{
    assert(incRowC==1);

    char transA = (incRowA==1) ? 'N' : 'T';
    char transB = (incRowB==1) ? 'N' : 'T';

    MKL_INT ldA = (incRowA==1) ? incColA : incRowA;
    MKL_INT ldB = (incRowB==1) ? incColB : incRowB;
    MKL_INT ldC = incColC;

    dgemm(&transA, &transB,
          &m, &n, &k,
          &alpha,
          A, &ldA,
          B, &ldB,
          &beta,
          C, &ldC);
}

} // namespace mkl
*/

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("#%4s %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();
/*
        mkl::gemm(C1.m, C1.n, A.n, ALPHA,
                  A.data, A.incRow, A.incCol,
                  B.data, B.incRow, B.incCol,
                  BETA,
                  C1.data, C1.incRow, C1.incCol);
*/
        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");
    }
}