Home     Back     Up     Next Browse Files

Content

GEMM für unterschiedliche Datentypen

Die Referenz-Implementierung der GEMM Operation zu verallgemeinern ist relativ einfach. Dies geschieht analog zu den vorigen Beispielen gecopy, geaxpy und gescal. Lediglich die Anzahl der Template-Parameter erhöht sich:

namespace refblas {

template <typename Index, typename Alpha, typename TA, typename TB,
          typename Beta, typename TC>
void
gemm(Index m, Index n, Index k,
     Alpha alpha,
     const TA *A, Index incRowA, Index incColA,
     const TB *B, Index incRowB, Index incColB,
     Beta beta,
     TC *C, Index incRowC, Index incColC)
{
    hpc::ulmblas::gescal(m, n, beta, C, incRowC, incColC);
    for (Index j=0; j<n; ++j) {
        for (Index l=0; l<k; ++l) {
            for (Index i=0; i<m; ++i) {
                C[i*incRowC+j*incColC] += alpha*A[i*incRowA+l*incColA]
                                               *B[l*incRowB+j*incColB];
            }
        }
    }
}

Etwas anspruchsvoller ist es, die geblockte GEMM Operation aufzubohren ohne die Effizienz zu verlieren. Insbesondere die optimierten Assembler Micro-Kernel setzen voraus, dass die Elemente von Paneelen aus A und B den gleichen Typ besitzen.

Andererseits ist es doch recht einfach, denn die Matrix-Multiplikation wird blockweise durchgeführt. Und bevor zwei Blöcke miteinander multipliziert werden, müssen diese in Puffer gepackt werden. Den Typ der beiden Puffer legen wir als den aus Session 10, Page 6 bekannten common_type für A, B und alpha fest. Wie bei gecopy wird beim Packen also eventuell ein upcast durchgeführt, wenn sich der Typ der Matrix-Elemente vom Typ des Puffer unterscheiden. Für die gepackten Blöcke kann also ein Micro-Kernel verwendet werden, der vorraussetzt, dass Paneele homogene Datentypen enthalten.

Zu beachten ist allerdings noch, dass sich die Elemente von C noch im Typ unterscheiden können. Auf diesen Fall gehen wir weiter unten ein. Wir betrachten die Anpassung des geblockten GEMM Algorithmus in folgender Reihenfolge:

Diese einzelnen Komponenten pflegen wir dabei in unsere neue Verzeichnisstruktur ein.

Festlegung von Blockgrößen durch Macros

Wir verwenden Macros, um für Blockgößen Default-Werte festzulegen, die wir aber beim Übersetzen eventuell abändern können. Dies geschieht nach dem Muster

//
// double precision(double)
//
#ifndef D_BLOCKSIZE_MR
#define D_BLOCKSIZE_MR 4
#endif

#ifndef D_BLOCKSIZE_NR
#define D_BLOCKSIZE_NR 4
#endif

#ifndef D_BLOCKSIZE_MC
#define D_BLOCKSIZE_MC 256
#endif

#ifndef D_BLOCKSIZE_KC
#define D_BLOCKSIZE_KC 256
#endif

#ifndef D_BLOCKSIZE_NC
#define D_BLOCKSIZE_NC 4096
#endif

Dabei steht D für double. Analog benutzen wir gemäß der BLAS Konvention Macros mit Prefix S für float, C für std::complex<float> und Z für std::complex<double>. Für andere Datentypen legen wir gemeinsame Default-Werte mit Prefix DEFAULT fest.

Aufgabe

Übernehmt das fertige config.h und kopiert es in hpc/ulmblas.

Festlegung von Blockgrößen durch Traits

Die oben definierten Macro verwenden wir wieder indirekt mit Hilfe einer Traits-Klasse. Damit kann die Blockgröße in Abhängigkeit eines Template-Parameter zur Compile-Zeit festgelegt werden:

template <typename Index, typename T>
void
foo(Index m, Index k, T *A)
{
    Index MR = hpc::ulmblas::BlockSize<T>::MR;
    /*
       ...
    */
}

Aufgabe

Die Traits-Klasse hpc::ulmblas::BlockSize packen wir in hpc/ulmblas/blocksize.h. Dabei sollt ihr folgende Vorlage um Spezialisierungen für float, double, std::complex<float> und std::complex<double> ergänzen:

#ifndef HPC_ULMBLAS_BLOCKSIZE_H
#define HPC_ULMBLAS_BLOCKSIZE_H 1

#include <complex>
#include <hpc/ulmblas/config.h>

namespace hpc { namespace ulmblas {

template <typename T>
struct BlockSize
{
    static const int MC = DEFAULT_BLOCKSIZE_MC;
    static const int KC = DEFAULT_BLOCKSIZE_KC;
    static const int NC = DEFAULT_BLOCKSIZE_NC;
    static const int MR = DEFAULT_BLOCKSIZE_MR;
    static const int NR = DEFAULT_BLOCKSIZE_NR;


    static_assert(MC>0 && KC>0 && NC>0 && MR>0 && NR>0, "Invalid block size.");
    static_assert(MC % MR == 0, "MC must be a multiple of MR.");
    static_assert(NC % NR == 0, "NC must be a multiple of NR.");
};

/*

   Special cases of BlockSize for float, double, ...

 */

} } // namespace ulmblas, hpc

#endif // HPC_ULMBLAS_BLOCKSIZE_H

Funktionen zum Packen von Blöcken

Als Vorlage zum Packen von Blöcken dient:

template <typename Index, typename T>
void
pack_A(Index mc, Index kc,
       const T *A, Index incRowA, Index incColA,
       T *p)
{
    Index MR = BlockSize<T>::MR;
    Index mp = (mc+MR-1) / MR;

    for (Index j=0; j<kc; ++j) {
        for (Index l=0; l<mp; ++l) {
            for (Index i0=0; i0<MR; ++i0) {
                Index i  = l*MR + i0;
                Index nu = l*MR*kc + j*MR + i0;
                p[nu]   = (i<mc) ? A[i*incRowA+j*incColA]
                                 : T(0);
            }
        }
    }
}

template <typename Index, typename T>
void
pack_B(Index kc, Index nc,
       const TB *B, Index incRowB, Index incColB,
       T *p)
{
    Index NR = BlockSize<T>::NR;
    Index np = (nc+NR-1) / NR;

    for (Index l=0; l<np; ++l) {
        for (Index j0=0; j0<NR; ++j0) {
            for (Index i=0; i<kc; ++i) {
                Index j  = l*NR+j0;
                Index nu = l*NR*kc + i*NR + j0;
                p[nu]   = (j<nc) ? B[i*incRowB+j*incColB]
                                 : T(0);
            }
        }
    }
}

Aufgaben

Etwas mehr Überblick und Weitsicht ist hier notwendig:

Macro-Kernel

Bei folgender Vorlage wurde ein Bug beseitigt. Da der Puffer C_ auf dem Stack liegt kann diese Werte wie NaNund inf beihalten. In diesem Fall wird der Inhalt beim Aufruf des Micro-Kernel mit beta=0 nicht mit Null Initialisiert. Mit std::fill_n(C_, MR*NR, T(0)) wird der Puffer vor jedem Aufruf des Micro-Kernel mit Null überschrieben. Mehr zu std::fill findet ihr hier.

template <typename Index, typename T>
void
mgemm(Index mc, Index nc, Index kc,
      T alpha,
      const T *A, const T *B,
      T beta,
      T *C, Index incRowC, Index incColC)
{
    Index MR = BlockSize<T>::MR;
    Index NR = BlockSize<T>::NR;
    T C_[BlockSize<T>::MR*BlockSize<T>::NR];

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

    for (Index j=0; j<np; ++j) {
        Index nr = (j!=np-1 || nr_==0) ? NR : nr_;

        for (Index i=0; i<mp; ++i) {
            Index mr = (i!=mp-1 || mr_==0) ? MR : mr_;

            if (mr==MR && nr==NR) {
                ugemm(kc, alpha,
                      &A[i*kc*MR], &B[j*kc*NR],
                      beta,
                      &C[i*MR*incRowC+j*NR*incColC],
                      incRowC, incColC);
            } else {
                // Fill buffer with zeros in case of infs and nans
                std::fill_n(C_, MR*NR, T(0));
                ugemm(kc, alpha,
                      &A[i*kc*MR], &B[j*kc*NR],
                      T(0),
                      C_, Index(1), MR);
                gescal(mr, nr, beta,
                       &C[i*MR*incRowC+j*NR*incColC],
                       incRowC, incColC);
                geaxpy(mr, nr, T(1), C_, Index(1), MR,
                       &C[i*MR*incRowC+j*NR*incColC],
                       incRowC, incColC);
            }
        }
    }
}

Aufgabe

Der Macro-Kernel erhält mit A und B gepackte Blöcke. Das Packen geschieht im Frame-Algorithmus und zwar so, dass A und B den gleichen Typ besitzen. Allerdings kann sich dieser vom Typ für beta und C unterscheiden. Ändert die Implementierung so ab, dass dies möglich ist. Überlegt euch dabei bereits welche Signatur vom Micro-Kernel im Allgemeinen unterstützt werden muss.

Der Macro-Kernel soll in hpc/ulmblas/mgemm.h mit entsprechenden Include-Guards und Namesräumen implementiert werden.

Micro-Kernel

Für den Fall, dass der Typ der Puffer-Elemente mit dem Typ für beta und Puffer-Elementen identisch ist (homogene Typen), kann die bisherige Implementierung verwendet werden:

template <typename Index, typename T>
void
ugemm(Index kc, T alpha,
      const T *A, const T *B,
      T beta,
      T *C, Index incRowC, Index incColC)
{
    const Index MR = BlockSize<T>::MR;
    const Index NR = BlockSize<T>::NR;
    T P[BlockSize<T>::MR*BlockSize<T>::NR];

    for (Index l=0; l<MR*NR; ++l) {
        P[l] = 0;
    }
    for (Index l=0; l<kc; ++l) {
        for (Index j=0; j<NR; ++j) {
            for (Index i=0; i<MR; ++i) {
                P[i+j*MR] += A[i+l*MR]*B[l*NR+j];
            }
        }
    }
    for (Index j=0; j<NR; ++j) {
        for (Index i=0; i<MR; ++i) {
            C[i*incRowC+j*incColC] *= beta;
            C[i*incRowC+j*incColC] += alpha*P[i+j*MR];
        }
    }
}

Für spezielle Element-Typen und Blöckgrößen MR und NR stellen wir zudem optimierte Assembler Implementierungen zur Verfügung.

Fälle mit inhomogenen Datenen Typen (bei denen beta oder Elemente von C einen anderen Typ als alpha, A oder B besitzen) führen wir deshalb auf den Fall von homogenen Typen zurück. Diese Fälle werden von einer verallgemeinerten hpc::ulmblas::ugemm Variante behandelt, die einen Puffer P der Dimension MR x NR auf dem Stack anlegt. Dieser ist vom gleichen Typ wie die Elemente von A und B. Die Operation wird dann wie folgt durchgeführt:

Entscheidend ist dabei, dass hier eventuell ein optimierter Micro-Kernel verwendet werden kann.

Aufgabe

Die Micro-Kernel landen im Verzeichnis hpc/ulmblas/kernels/:

Die Liste an Micro-Kernel wird sich regelmäßig erweitern. Deshalb soll es möglich sein, dass mit einem einzigen #include alle verfügbaren Micro-Kernel eingebunden werden. Dazu dient folgendes Include-File:

session12/page01/hpc/ulmblas/kernels/ugemm.h 
#ifndef HPC_ULMBLAS_KERNELS_UGEMM_H
#define HPC_ULMBLAS_KERNELS_UGEMM_H 1

#include <hpc/ulmblas/config.h>
#include <hpc/ulmblas/kernels/ugemm_ref.h>
#include <hpc/ulmblas/kernels/ugemm_buf.h>

#endif // HPC_ULMBLAS_KERNELS_UGEMM_H

Frame-Algorithmus

Als Vorlage für den Frame-Algorithmus verwenden wir:

template <typename Index, typename T>
void
gemm(Index m, Index n, Index k,
     T alpha,
     const T *A, Index incRowA, Index incColA,
     const T *B, Index incRowB, Index incColB,
     T beta,
     T *C, Index incRowC, Index incColC)
{
    Index MC = BlockSize<T>::MC;
    Index NC = BlockSize<T>::NC;
    Index KC = BlockSize<T>::KC;

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

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

    T *A_ = new T[MC*KC];
    T *B_ = new T[KC*NC];

    if (alpha==T(0) || k==0) {
        gescal(m, n, beta, C, incRowC, incColC);
        return;
    }

    for (Index j=0; j<nb; ++j) {
        Index nc = (j!=nb-1 || nc_==0) ? NC : nc_;

        for (Index l=0; l<kb; ++l) {
            Index   kc    = (l!=kb-1 || kc_==0) ? KC : kc_;
            T beta_ = (l==0) ? beta : T(1);

            pack_B(kc, nc,
                   &B[l*KC*incRowB+j*NC*incColB],
                   incRowB, incColB,
                   B_);

            for (Index i=0; i<mb; ++i) {
                Index mc = (i!=mb-1 || mc_==0) ? MC : mc_;

                pack_A(mc, kc,
                       &A[i*MC*incRowA+l*KC*incColA],
                       incRowA, incColA,
                       A_);

                mgemm(mc, nc, kc,
                      alpha, A_, B_, beta_,
                      &C[i*MC*incRowC+j*NC*incColC],
                      incRowC, incColC);
            }
        }
    }
    delete [] A_;
    delete [] B_;
}

Aufgaben

Test

Da wir die Matrix-Klassen noch noch angepasst haben und auch andere Hilfsmittel noch übertragen werden müssen, sieht der low-level Benchmark wieder etwas umständlicher aus:

session12/page01/hpc/tests/gemm.cc 
#include <chrono>
#include <complex>
#include <cmath>
#include <cstdio>
#include <limits>
#include <random>
#include <hpc/ulmblas/gecopy.h>
#include <hpc/ulmblas/gemm.h>

//
//  Reference implementation for GEMM
//
namespace refblas {

template <typename Index, typename Alpha, typename TA, typename TB,
          typename Beta, typename TC>
void
gemm(Index m, Index n, Index k,
     Alpha alpha,
     const TA *A, Index incRowA, Index incColA,
     const TB *B, Index incRowB, Index incColB,
     Beta beta,
     TC *C, Index incRowC, Index incColC)
{
    hpc::ulmblas::gescal(m, n, beta, C, incRowC, incColC);
    for (Index j=0; j<n; ++j) {
        for (Index l=0; l<k; ++l) {
            for (Index i=0; i<m; ++i) {
                C[i*incRowC+j*incColC] += alpha*A[i*incRowA+l*incColA]
                                               *B[l*incRowB+j*incColB];
            }
        }
    }
}

} // namespace refblas


//
//  Random initializer for general matrices: real and complex valued
//
template <typename Index, typename T>
void
randomInit(Index m, Index n, T *A, Index incRowA, Index incColA)
{
    std::random_device                  random;
    std::default_random_engine          mt(random());
    std::uniform_real_distribution<T>   uniform(-100,100);

    for (Index i=0; i<m; ++i) {
        for (Index j=0; j<n; ++j) {
            A[i*incRowA+j*incColA] = uniform(mt);
        }
    }
}

template <typename Index, typename T>
void
randomInit(Index m, Index n, std::complex<T> *A, Index incRowA, Index incColA)
{
    std::random_device                  random;
    std::default_random_engine          mt(random());
    std::uniform_real_distribution<T>   uniform(-100,100);

    for (Index i=0; i<m; ++i) {
        for (Index j=0; j<n; ++j) {
            A[i*incRowA+j*incColA] = std::complex<T>(uniform(mt), uniform(mt));
        }
    }
}

//
// Timer
//
template <typename T>
struct WallTime
{
    void
    tic()
    {
        t0 = std::chrono::high_resolution_clock::now();
    }

    T
    toc()
    {
        using namespace std::chrono;

        elapsed = high_resolution_clock::now() - t0;
        return duration<T,seconds::period>(elapsed).count();
    }

    std::chrono::high_resolution_clock::time_point t0;
    std::chrono::high_resolution_clock::duration   elapsed;
};


//
//  Compute 1-norm of difference A - B
//
template <typename Index, typename T>
double
asumDiff(Index m, Index n,
         const T *A, Index incRowA, Index incColA,
         const T *B, Index incRowB, Index incColB)
{
    double diff = 0;

    for (Index i=0; i<m; ++i) {
        for (Index j=0; j<n; ++j) {
            diff += std::abs(A[i*incRowA+j*incColA] - B[i*incRowB+j*incColB]);
        }
    }
    return diff;
}

//
// Compute 1-norm
//
template <typename Index, typename T>
double
asum(Index m, Index n,
     const T *A, Index incRowA, Index incColA)
{
    double result = 0;

    for (Index i=0; i<m; ++i) {
        for (Index j=0; j<n; ++j) {
            result += std::abs(A[i*incRowA+j*incColA]);
        }
    }
    return result;
}

//
//  Compute residual for matrix-product
//
template <typename Index, typename Alpha, typename TA, typename TB,
          typename Beta, typename TC>
double
res_gemm(Index m, Index n, Index k,
         Alpha alpha,
         const TA *A, Index incRowA, Index incColA,
         const TB *B, Index incRowB, Index incColB,
         Beta beta,
         TC *C1, Index incRowC1, Index incColC1,
         TC *C2, Index incRowC2, Index incColC2)
{
    double aNorm = asum(m, k, A, incRowA, incColA) * std::abs(alpha);
    double bNorm = asum(k, n, B, incRowB, incColB);
    double cNorm = asum(m, n, C2, incRowC2, incColC2);
    double aDiff = asumDiff(m, n,
                            C1, incRowC1, incColC1,
                            C2, incRowC2, incColC2);
    // Using eps for double gives upper bound in case elements have lower
    // precision.
    double eps = std::numeric_limits<double>::epsilon();
    double res = aDiff/(aNorm*bNorm*cNorm*eps*std::max(std::max(m,n),k));
    return res;
}

int
main()
{
    // Element type for A, B and C
    typedef float                TA;
    typedef float                TB;
    typedef std::complex<double> TC;

    // Type for scalar factors
    typedef float   Alpha;
    typedef float   Beta;

    // Max dimensions of A, B and C 
    long max_m  = 7000;
    long max_n  = 7000;
    long max_k  = 7000;

    // Storage order of A
    long incRowA = 1;
    long incColA = max_m;

    // Storage order of B
    long incRowB = max_n;
    long incColB = 1;

    // Storage order of C
    long incRowC = max_n;
    long incColC = 1;

    // Allocate A, B, C1, C2
    TA *A  = new TA[max_m*max_k];
    TB *B  = new TB[max_k*max_n];
    TC *C1 = new TC[max_m*max_n];
    TC *C2 = new TC[max_m*max_n];

    // Init all matrices
    randomInit(max_m, max_k, A, incRowA, incColA);
    randomInit(max_k, max_n, B, incRowB, incColB);
    randomInit(max_m, max_n, C1, incRowC, incColC);
    hpc::ulmblas::gecopy(max_m, max_n,
                         C1, incRowC, incColC,
                         C2, incRowC, incColC);

    // Init scalar factors
    const Alpha alpha(1.5);
    const Beta  beta(2.5);

    // Header for benchmark
    printf("%5s %5s %5s ", "m", "n", "k");
    printf("%20s %9s", "refColMajor: t", "MFLOPS");
    printf("%20s %9s %9s", "blocked GEMM: t", "MFLOPS", "diff");
    printf("\n");

    WallTime<double> wallTime;

    for (long m=200, n=200, k=200;
         m <=max_m && n<=max_n && k<=max_k;
         m+=100, n+=100, k+=100)
    {
        printf("%5ld %5ld %5ld ", m, n, k);

        wallTime.tic();
        refblas::gemm(m, n, k, alpha,
                      A, incRowA, incColA,
                      B, incRowB, incColB,
                      beta,
                      C1, incRowC, incColC);
        double t = wallTime.toc();
        printf("%20.4lf %9.2lf", t, 2.*m/1000*n/1000*k/t);

        wallTime.tic();
        hpc::ulmblas::gemm(m, n, k, alpha,
                           A, incRowA, incColA,
                           B, incRowB, incColB,
                           beta,
                           C2, incRowC, incColC);
        t = wallTime.toc();
        double res = res_gemm(m, n, k, alpha,
                              A, incRowA, incColA,
                              B, incRowB, incColB,
                              beta,
                              C1, incRowC, incColC,
                              C2, incRowC, incColC);

        printf("%20.4lf %9.2lf %9.1e", t, 2.*m/1000*n/1000*k/t, res);
        printf("\n");
    }
}

Aufgaben

Home     Back     Up     Next