========================
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:

  - Mit `git clone https://github.com/michael-lehn/ulmBLAS.git` könnt ihr euch
    das `ulmBLAS` Repository holen.
  - Mit `git tag -l` bekommt ihr eine Liste aller Etiketten, engl. Tags.  Jeder
    Lösungsvorschlag wurde von mir mit einer Etikette/Tag versehen.
  - Mit `git checkout tags/v0.0.10` kann man jetzt den Zustand des Tags
    "v0.0.10" herstellen.  Das ist der Zustand des Repositories als der Tag
    erzeugt wurde.

ulmBLAS Level 1: Step by Step
=============================
Wir beginnen im Zustand von Tag `v0.0.1`:
 - In `level1/` sind stubs zu allen notwendigen BLAS Routinen enthalten.
 - Es gibt ein Top-Level Makefile, das Makefiles in den Unterverzeichnissen
   aufruft:
     - `make` ruft in allen Unterverzeichnissen `make` auf.
     - `make clean` ruft in allen Unterverzeichnissen `make clean` auf.
     - `make check` ruft im Unterverzeichnis `test` nur `make check_ulm` auf.
       Wenn man in `test` die Netlib Referenz-Implementierung testen möchte
       muss man das von Hand tun.

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:

 - Portiert `ddot.f` (Lösungsvorschlag in `v0.0.2`)
 - Portiert `dasum.f` (Lösungsvorschlag in `v0.0.3`)
 - Portiert `daxpy.f` (Lösungsvorschlag in `v0.0.4`)
 - Portiert `dcopy.f` (Lösungsvorschlag in `v0.0.5`)
 - Portiert `dnrm2.f` (Lösungsvorschlag in `v0.0.6`)
 - Portiert `dscal.f` (Lösungsvorschlag in `v0.0.7`)
 - Portiert `idamax.f` (Lösungsvorschlag in `v0.0.8`)
 - Und die übrigen Given Rotationen (Lösungsvorschlag in `v0.0.9`)

Hinweise
--------
Ihr werdet verschieden Funktionen aus `math.h` benötigen:
 - `fabs`, `sqrt`, `pow`. Schaut euch dazu ggf. die Manual Seiten an.  Ihr
   werdet dann merken, dass zum Beispiel `pow(double, int)` nicht unterstützt
   wird.
 - Es lohnt sich, dass bei jeder Funktion zu Beginn die skalaren Argumente
   einmal explizit dereferenziert werden.  Damit wird der Code lesbarer.


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:

---- CODE(type=c) --------------------------------------------------------------
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:
 - Der C-Schnittstelle wird ein `ULM_` Präfix verpasst und
 - der Fortran-Schnittstelle wird ein `_` Suffix verpasst.

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:

---- CODE(type=c) --------------------------------------------------------------
#define F77BLAS(x) x##_
#define ULMBLAS(x) ULM_##x
--------------------------------------------------------------------------------

Damit sieht der Aufbau eine ulmBLAS Routine so aus:

---- CODE(type=c) --------------------------------------------------------------
#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:

---- CODE(type=c) --------------------------------------------------------------
#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:*
  - Passt die Funktionen in `level1/` obigem Muster an.
  - Wie muss im Makefile `CFLAGS` angepasst werden damit alles übersetzt wird?

:links: GitHub -> https://github.com/michael-lehn

:navigate: __up__    -> doc:index
           __back__  -> doc:day06/page01
           __next__  -> doc:day06/page03