Geblocktes GEMM

Implementiert den geblockten GEMM Algorithmus aus der Vorlesung. Als Vorlage könnt ihr folgenden Code verwenden oder den eigenen Code entsprechend ergänzen. Zu implementieren sind:

Die verschiedenen GEMM Algorithmen

Alle Algorithmen berechnen die Operation

\[C \leftarrow \beta\,C + \alpha\, A B\quad\text{mit}\;A \in \mathbb{M}^{m \times k},\;B \in \mathbb{M}^{k \times n}\;\text{und}\;C \in \mathbb{M}^{m \times n}.\]

Unterschiede ergeben sich aus verschiedenen Gründen:

Die Algorithmen werden also immer spezieller.

Micro-Kernel

Macro-Kernel

Vorlage

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <sys/times.h>
#include <unistd.h>

#ifndef COLMAJOR
#define COLMAJOR 1
#endif

#ifndef MAXDIM_M
#define MAXDIM_M    4000
#endif

#ifndef MAXDIM_N
#define MAXDIM_N    4000
#endif

#ifndef MAXDIM_K
#define MAXDIM_K    4000
#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


//==============================================================================
// bench
//==============================================================================

namespace bench {

//------------------------------------------------------------------------------
// Auxiliary data for benchmarking
//------------------------------------------------------------------------------

double A[MAXDIM_M*MAXDIM_K];
double B[MAXDIM_K*MAXDIM_N];

double C1[MAXDIM_M*MAXDIM_N];
double C2[MAXDIM_M*MAXDIM_N];
double C3[MAXDIM_M*MAXDIM_N];

//------------------------------------------------------------------------------
// Auxiliary functions for benchmarking
//------------------------------------------------------------------------------

double
walltime()
{
   struct tms    ts;
   static double ClockTick=0.0;

   if (ClockTick==0.0) {
        ClockTick = 1.0 / ((double) sysconf(_SC_CLK_TCK));
   }
   return ((double) times(&ts)) * ClockTick;
}

void
initMatrix(long m, long n, double *A, long incRowA, long incColA)
{
    for (long j=0; j<n; ++j) {
        for (long i=0; i<m; ++i) {
            A[i*incRowA+j*incColA] = ((double)rand() - RAND_MAX/2)*200/RAND_MAX;
        }
    }
}

double
asumDiffMatrix(long m, long n,
               const double *A, long incRowA, long incColA,
               double *B, long incRowB, long incColB)
{

    double asum = 0;

    for (long j=0; j<n; ++j) {
        for (long i=0; i<m; ++i) {
            asum += fabs(B[i*incRowB+j*incColB] - A[i*incRowA+j*incColA]);
        }
    }
    return asum;
}

} // namespace bench


//==============================================================================
// ulmBLAS
//==============================================================================

namespace ulmBLAS {

//------------------------------------------------------------------------------
// A <- B
//------------------------------------------------------------------------------

void
gecopy(long m, long n,
       const double *A, long incRowA, long incColA,
       double *B, long incRowB, long incColB)
{
    for (long j=0; j<n; ++j) {
        for (long i=0; i<m; ++i) {
            B[i*incRowB+j*incColB] = A[i*incRowA+j*incColA];
        }
    }
}

//------------------------------------------------------------------------------
// Y <- Y + alpha*Y#
//------------------------------------------------------------------------------

void
geaxpy(long m, long n,
       double alpha,
       const double *X, long incRowX, long incColX,
       double       *Y, long incRowY, long incColY)
{
    for (long i=0; i<m; ++i) {
        for (long j=0; j<n; ++j) {
            Y[i*incRowY+j*incColY] += alpha*X[i*incRowX+j*incColX];
        }
    }
}

//------------------------------------------------------------------------------
// A <- alpha * A
//------------------------------------------------------------------------------

void
gescal(long m, long n,
       double alpha,
       double *X, long incRowX, long incColX)
{
    if (alpha!=1.0) {
        for (long i=0; i<m; ++i) {
            for (long j=0; j<n; ++j) {
                X[i*incRowX+j*incColX] *= alpha;
            }
        }
    }
}

} // namespace ulmBLAS

//==============================================================================
// refColMajor
//==============================================================================

namespace refColMajor {

void
gemm(long m, long n, long k,
     double alpha,
     const double *A, long incRowA, long incColA,
     const double *B, long incRowB, long incColB,
     double beta,
     double *C, long incRowC, long incColC)
{
    ulmBLAS::gescal(m, n, beta, C, incRowC, incColC);
    for (long j=0; j<n; ++j) {
        for (long l=0; l<k; ++l) {
            for (long i=0; i<m; ++i) {
                C[i*incRowC+j*incColC] += alpha*A[i*incRowA+l*incColA]
                                               *B[l*incRowB+j*incColB];
            }
        }
    }
}

} // namespace refColMajor

//==============================================================================
// simpleBuffer
//==============================================================================

#ifndef SIMPLE_PUFFER_M_C
#define SIMPLE_PUFFER_M_C 256
#endif

#ifndef SIMPLE_PUFFER_K_C
#define SIMPLE_PUFFER_K_C 256
#endif

#ifndef SIMPLE_PUFFER_N_C
#define SIMPLE_PUFFER_N_C 1024
#endif

namespace simpleBuffer {

double A_[SIMPLE_PUFFER_M_C*SIMPLE_PUFFER_K_C];
double B_[SIMPLE_PUFFER_K_C*SIMPLE_PUFFER_N_C];
double C_[SIMPLE_PUFFER_M_C*SIMPLE_PUFFER_N_C];

void
gemm(long m, long n, long k,
     double alpha,
     const double *A, long incRowA, long incColA,
     const double *B, long incRowB, long incColB,
     double beta,
     double *C, long incRowC, long incColC)
{
    long i, j, l;

    long M_C = SIMPLE_PUFFER_M_C;
    long N_C = SIMPLE_PUFFER_N_C;
    long K_C = SIMPLE_PUFFER_K_C;

    long mb = m / M_C;
    long nb = n / N_C;
    long kb = k / K_C;

    long mr = m % M_C;
    long nr = n % N_C;
    long kr = k % K_C;

    using ulmBLAS::geaxpy;
    using ulmBLAS::gecopy;
    using ulmBLAS::gescal;
    using refColMajor::gemm;

    gescal(m, n, beta, C, incRowC, incColC);

    for (j=0; j<=nb; ++j) {
        long N = (j<nb || nr==0) ? N_C : nr;

        for (l=0; l<=kb; ++l) {
            long K = (l<kb || kr==0) ? K_C : kr;

            gecopy(K, N,
                   &B[l*K_C*incRowB+j*N_C*incColB], incRowB, incColB,
                   B_, 1, K_C);

            for (i=0; i<=mb; ++i) {
                long M = (i<mb || mr==0) ? M_C : mr;

                gecopy(M, K,
                       &A[i*M_C*incRowA+l*K_C*incColA], incRowA, incColA,
                       A_, 1, M_C);

                gemm(M, N, K,
                     1.0,
                     A_, 1, M_C,
                     B_, 1, K_C,
                     0.0,
                     C_, 1, M_C);

                geaxpy(M, N, alpha,
                       C_, 1, M_C,
                       &C[i*M_C*incRowC+j*N_C*incColC], incRowC, incColC);
            }
        }
    }
}

} // namespace simpleBuffer

//==============================================================================
// blocked
//==============================================================================


#ifndef BLOCKED_MC
#define BLOCKED_MC 256
#endif

#ifndef BLOCKED_KC
#define BLOCKED_KC 256
#endif

#ifndef BLOCKED_NC
#define BLOCKED_NC 1024
#endif

#ifndef BLOCKED_MR
#define BLOCKED_MR 2
#endif

#ifndef BLOCKED_NR
#define BLOCKED_NR 8
#endif

namespace blocked {

double A_[BLOCKED_MC*BLOCKED_KC];
double B_[BLOCKED_KC*BLOCKED_NC];

void
pack_A(long mc, long kc,
       const double *A, long incRowA, long incColA,
       double *p)
{
    long MR = BLOCKED_MR;
    long mp = (mc+MR-1) / MR;

    // ...
}

void
pack_B(long kc, long nc,
       const double *B, long incRowB, long incColB,
       double *p)
{
    long NR = BLOCKED_NR;
    long np = (nc+NR-1) / NR;

    // ...
}

void
ugemm(long kc, double alpha,
      const double *A, const double *B,
      double beta,
      double *C, long incRowC, long incColC)
{
    double P[BLOCKED_MR*BLOCKED_NR];
    long MR = BLOCKED_MR;
    long NR = BLOCKED_NR;

    // ...
}

void
mgemm(long mc, long nc, long kc,
      double alpha,
      const double *A, const double *B,
      double beta,
      double *C, long incRowC, long incColC)
{
    double C_[BLOCKED_MR*BLOCKED_NR];
    long MR = BLOCKED_MR;
    long NR = BLOCKED_NR;

    long mp  = (mc+MR-1) / MR;
    long np  = (nc+NR-1) / NR;
    long mr_ = mc % MR;
    long nr_ = nc % NR;

    // ...
}

void
gemm(long m, long n, long k,
     double alpha,
     const double *A, long incRowA, long incColA,
     const double *B, long incRowB, long incColB,
     double beta,
     double *C, long incRowC, long incColC)
{
    long MC = BLOCKED_MC;
    long NC = BLOCKED_NC;
    long KC = BLOCKED_KC;

    long mb = (m+MC-1) / MC;
    long nb = (n+NC-1) / NC;
    long kb = (k+KC-1) / KC;

    long mc_ = m % MC;
    long nc_ = n % NC;
    long kc_ = k % KC;

    // ...
}

} // namespace blocked

int
main()
{
    using namespace bench;
    using namespace ulmBLAS;
    using namespace std;

    initMatrix(MAXDIM_M, MAXDIM_K, A, 1, MAXDIM_M);
    initMatrix(MAXDIM_K, MAXDIM_N, B, 1, MAXDIM_K);
    initMatrix(MAXDIM_M, MAXDIM_N, C1, 1, MAXDIM_M);
    gecopy(MAXDIM_M, MAXDIM_N, C1, 1, MAXDIM_M, C2, 1, MAXDIM_M);
    gecopy(MAXDIM_M, MAXDIM_N, C1, 1, MAXDIM_M, C3, 1, MAXDIM_M);

    // Header-Zeile fr die Ausgabe
    printf("%5s %5s %5s ", "m", "n", "k");
    printf("%5s %5s ", "IRA", "ICA");
    printf("%5s %5s ", "IRB", "ICB");
    printf("%5s %5s ", "IRC", "ICC");
    printf("%20s %9s", "refColMajor: t", "MFLOPS");
    printf("%20s %9s %9s", "refSimpleBuffer: t", "MFLOPS", "diff");
    printf("\n");

    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 t0, t1, diff;

        long incRowA = (COLMAJOR) ? 1 : k;
        long incColA = (COLMAJOR) ? m : 1;

        long incRowB = (COLMAJOR) ? 1 : n;
        long incColB = (COLMAJOR) ? k : 1;

        long incRowC = (COLMAJOR) ? 1 : n;
        long incColC = (COLMAJOR) ? m : 1;

        printf("%5ld %5ld %5ld ", m, n, k);
        printf("%5ld %5ld ", incRowA, incColA);
        printf("%5ld %5ld ", incRowB, incColB);
        printf("%5ld %5ld ", incRowC, incColC);

        t0 = walltime();
        refColMajor::gemm(m, n, k, ALPHA,
                          A, incRowA, incColA,
                          B, incRowB, incColB,
                          BETA,
                          C1, incRowC, incColC);
        t1 = walltime() - t0;
        printf("%20.4lf %9.2lf", t1, 2.*m/1000*n/1000*k/t1);

        t0 = walltime();
        simpleBuffer::gemm(m, n, k, ALPHA,
                           A, incRowA, incColA,
                           B, incRowB, incColB,
                           BETA,
                           C2, incRowC, incColC);
        t1 = walltime() - t0;
        diff = asumDiffMatrix(m, n,
                              C1, incRowC, incColC,
                              C2, incRowC, incColC);
        printf("%20.4lf %9.2lf %9.1e", t1, 2.*m/1000*n/1000*k/t1, diff);

        t0 = walltime();
        blocked::gemm(m, n, k, ALPHA,
                      A, incRowA, incColA,
                      B, incRowB, incColB,
                      BETA,
                      C3, incRowC, incColC);
        t1 = walltime() - t0;
        diff = asumDiffMatrix(m, n,
                              C1, incRowC, incColC,
                              C3, incRowC, incColC);
        printf("%20.4lf %9.2lf %9.1e", t1, 2.*m/1000*n/1000*k/t1, diff);
        printf("\n");
    }
    return 0;
}