ulmBLAS: Tests Bestehen!

Das Grundgerüst unserer BLAS Level 1 Implementierung sollte jetzt stehen. Diese besteht aber hauptsächlich aus Funktionsrümpfen (“Stubs”) die noch nicht tun was sie sollen. Das sollt ihr jetzt ändern. Schritt für Schritt soll jetzt die Referenz-Implementierung portiert werden. Um die Optimierung kümmern wir uns später. Lösungsvorschläge für die einzelnen Schritte kann man sich über GitHub besorgen.

Lösungsvorschläge von GitHub beziehen

Wir werden uns später genauer mit dem git Tool beschäftigen. Für den Moment sollte eine Benutzung nach Rezept genügen:

ulmBLAS Level 1: Step by Step

Wir beginnen im Zustand von Tag v0.0.1:

Aufgaben

Portiert nun eine BLAS Routine nach der anderen. Ihr werdet merken, dass bestimmte Strukturen wiederkehren. Zum Beispiel werden immer die Spezialfälle bei denen alle Strides gleich Eins sind speziell behandelt. Das sind die Fälle in denen eine weitere Optimierung möglich ist. Euer Port sollte deshalb diese Fälle ebenfalls gesondert behandeln.

Nach jedem Port sollte ein zusätzlicher Test bestanden werden:

Hinweise

Ihr werdet verschieden Funktionen aus math.h benötigen:

ulmBLAS Schnittstellen

Im Moment hat ulmBLAS nur eine Fortran 77 Schnittstelle. Das heisst alle Argumente werden als Zeiger übergeben. Intern werden dann skalare Argumente explizit dereferenziert. Anschliessen wird dann nur mit den dereferenzierten Argumenten weitergearbeitet. Es wird sich lohnen hier eine Schnittstelle dazwischen einzuziehen. Dabei entsteht eine zusätzliche C-Schnittstelle:

double
ULM_ddot(const int     n,
         const double  *x,
         const int     incX,
         const double  *y,
         const int     incY)
{
    // Implementierung der eigentlichen Funktionalität
}

double
ddot_(const int     *_n,
      const double  *x,
      const int     *_incX,
      const double  *y,
      const int     *_incY)
{
    // Dereferenziere _n, _incX, _incY und rufe dann ULMBLAS(ddot) auf
}

Die Namensgebung ist hier handgestrickt:

Mit Hilfe eines Macros kann man das flexibler gestalten. Wir definieren dies in einem Include File im obersten Verzeichnis. CFLAGS müssen dann den eventuell in den Makefiles angepasst werden). Aufgabe: findet heraus was das ## tut:

#define F77BLAS(x) x##_
#define ULMBLAS(x) ULM_##x

Damit sieht der Aufbau eine ulmBLAS Routine so aus:

#inlcude <ulmblas.h>

double
ULMBLAS(ddot)(const int     n,
              const double  *x,
              const int     incX,
              const double  *y,
              const int     incY)
{
    // ...
}

double
F77BLAS(ddot)(const int     *_n,
              const double  *x,
              const int     *_incX,
              const double  *y,
              const int     *_incY)
{
    // Dereferenziere _n, _incX, _incY und rufe dann ULMBLAS(ddot) auf
}

Und die fertige Funktion zum Beispiel so:

#include <ulmblas.h>

double
ULMBLAS(ddot)(const int     n,
              const double  *x,
              const int     incX,
              const double  *y,
              const int     incY)
{
//
//  Local scalars
//
    double result = 0.0;

    int    i, m;

//
//  Quick return if possible
//
    if (n==0) {
        return 0;
    }
    if (incX==1 && incY==1) {
//
//      Code for both increments equal to 1
//
        m = n % 5;
        if (m!=0) {
            for (i=0; i<m; ++i) {
                result += x[i] * y[i];
            }
            if (n<5) {
                return result;
            }
        }
        for (i=m; i<n; i+=5) {
            result += x[i  ] * y[i  ] ;
            result += x[i+1] * y[i+1];
            result += x[i+2] * y[i+2];
            result += x[i+3] * y[i+3];
            result += x[i+4] * y[i+4];
        }
    } else {
//
//      Code for unequal increments or equal increments not equal to 1
//
        if (incX<0) {
            x -= incX*(n-1);
        }
        if (incY<0) {
            y -= incY*(n-1);
        }
        for (i=0; i<n; ++i, x+=incX, y+=incY) {
            result += (*x) * (*y);
        }
    }
    return result;
}

double
F77BLAS(ddot)(const int     *_n,
              const double  *x,
              const int     *_incX,
              const double  *y,
              const int     *_incY)
{
//
//  Dereference scalar parameters
//
    int n    = *_n;
    int incX = *_incX;
    int incY = *_incY;

    return ULMBLAS(ddot)(n, x, incX, y, incY);
}

Aufagben: