=======================================
GEMV: Berechnung mit Vektor-Operationen                             [TOC]
=======================================

In der letzten Vorlesung wurde die BLAS Operation _gemv_ für das allgemeine
Matrix-Vektor Produkt

---- LATEX --------------------------
y \leftarrow \beta\, y + \alpha\, A x
-------------------------------------

implementiert.  Unabhängig davon, ob die Matriz zeilen- oder spaltenweise im
Speicher liegt, kann unsere Implementierung in etwa die gleiche Performance
erreichen.  Wir haben dies erreicht, indem zwei algorithmische Varianten für
die Operation realisiert wurden:

---- BOX --------------------------------------------------
= $\text{If}\; \beta \neq 1$
    = $\text{For}\; i=1,\dots,m$
        = $y_i \leftarrow \beta\, y_i$
= $\text{If}\; \text{incRow}_\text{A} > \text{incCol}_\text{A}$
    = $\text{For}\; i=1,\dots,m$
        = $\text{For}\; j=1,\dots,n$
            = $y_i \leftarrow y_i + \alpha\, a_{i,j} x_j$
= $\text{Else}$
    = $\text{For}\; j=1,\dots,n$
        = $\text{For}\; i=1,\dots,m$
            = $y_i \leftarrow y_i + \alpha\, a_{i,j} x_j$
-----------------------------------------------------------

Vorteil einer komponentenweisen Formulierung
============================================
Der Vorteil bei der komponentenweisen Formulierung ist, dass wir direkt sehen
in welcher Reihenfolge auf einzelne Matrix oder Vektorelemente zugegriffen wird.
Dadurch war uns wiederum klar in welcher Reihenfolge auf welche Adressen in
Speicher zugegriffen wird.  Dies haben wir schließlich durch die
Fallunterscheidung für die Cache-Optimierung ausgenutzt.


Formulierung des Algorithmus mit Vektor-Operationen
===================================================
Wir behandeln die GEMV-Operation in zwei Schritten:

- Zuerst skalieren wir $y$ mit $\beta$:

  ---- LATEX ---------
  y \leftarrow \beta y
  --------------------

- Dann führen wir eine Aufdatierung mit dem Matrix-Vektor Produkt durch

  ---- LATEX ------------
  y \leftarrow \alpha A x
  -----------------------

Für die Aufdatierung 


-  *Zeilen-orientierte Variante*

   Die Matrix $A$ wird zeilenweise betrachtet:

   ---- LATEX ---------------------------------
   A = \begin{pmatrix}
       \mathbf{a}_1^T \\
       \vdots \\
       \mathbf{a}_m^T
       \end{pmatrix}
   --------------------------------------------

   Damit wird das Matrix-Vektor Produkte wohldefiniert ist, muss $y$
   komponentenweise betrachtet werden und durch Skalar-Produkte berechnet
   werden:

   ---- LATEX ---------------------------------
   \begin{pmatrix}
   y_1 \\
   \vdots \\
   y_m
   \end{pmatrix}
   \leftarrow
   \beta
   \begin{pmatrix}
   y_1 \\
   \vdots \\
   y_m
   \end{pmatrix}
   + \alpha
   \begin{pmatrix}
   \mathbf{a}_1^T x \\
   \vdots \\
   \mathbf{a}_m^T x
   \end{pmatrix}
   --------------------------------------------

-  *Spalten-orientierte Variante*

   Die Matrix $A$ wird spaltenweise betrachtet:

   ---- LATEX ---------------------------------
   A = \begin{pmatrix}
       \mathbf{a}_1 & \dots & \mathbf{a}_n
       \end{pmatrix}
   --------------------------------------------

   Damit wird das Matrix-Vektor Produkte wohldefiniert ist, muss $x$
   komponentenweise betrachtet werden.  Der Vektor $y$ wird dann als
   eine Linearkombination berechnet:

   ---- LATEX ---------------------------------
   y \leftarrow \beta\, y
   + \alpha\, \mathbf{a}_1 x_1
   + \dots
   + \alpha \,\mathbf{a}_n x_n
   --------------------------------------------

Beide Varianten können wir wieder in einem Frame-Algorithmus ausgewählt
werden:

---- BOX --------------------------------------------------------
= $y \leftarrow \beta\, y$
= $\text{If}\; \text{incRow}_\text{A} > \text{incCol}_\text{A}$
    = $\text{For}\; i=1,\dots,m$
        = $y_i \leftarrow y_i + \alpha\, \mathbf{a}_i^T x$
= $\text{Else}$
    = $\text{For}\; j=1,\dots,n$
        = $y \leftarrow y + \alpha\, \mathbf{a}_{j} x_j$
-----------------------------------------------------------------


Vorteil einer Formulierung mit Vektor-Operationen
=================================================
Der Algorithmus lässt sich kompakter aufschreiben.  Dabei haben wir
drei Kern-Operationen isoliert, die unabhängig voneinander umgesetzt
und optimiert werden können:

- Vektor-Skalierung

  ---- LATEX --------------------------------------
  \text{scal}:\quad  x \; \leftarrow \alpha\; x
  -------------------------------------------------

- Skalar-Produkt

  ---- LATEX --------------------------------------
  \text{dot}:\quad  \alpha \; \leftarrow x^T y
  -------------------------------------------------

- Vektor-Aufdatierung (_alpha x plus y_)

  ---- LATEX --------------------------------------
  \text{axpy}:\quad  y \; \leftarrow y + \alpha\, x
  -------------------------------------------------

Außerdem können wir dadurch einen Cache-optimierten Algorithmus für das
Matrix-Vektor Produkt durch eine einfache mathematische Betrachtung des
Problems herleiten:

 - Ist die Matrix zeilenweise gespeichert, dann betrachten wir die Matrix
   zeilenweise.
 - Ist die Matrix spaltenweise gespeichert, dann betrachten wir die Matrix
   spaltenweise.

Die dabei gewonnen Denk- und Vorgehensweise werden wir auf andere Probleme
im Bereich der numerischen Linearen Algebra übertragen.  Zunächst müssen wir
uns aber überzeugen, dass der zusätzliche Overhead durch Funktionsaufrufe für
die Effizienz keine Rolle spielt.


Aufgaben
========

- Fomuliert Algorithmen für _scal_, _dot_ und _axpy_.
- Implementiert die Algorithmen für Vektoren mit Elemeneten vom Typ _double_.
  Nennt diese `dscal_ulm`, `ddot_ulm` und `daxpy_ulm`.
- Implementiert für die GEMV Operation den obigen Frame-Algorithmus mit
  Vektor-Operationen.

Vorlage
=======
Im Pool E44 ist eine frei verfügbare Version der Intel MKL installiert.  Mit
dieser soll die Performance eurer Implementierung verglichen werden.  Die
Lösung von Session 4 wurde um diesen Benchmark ergänzt.  Außerdem stellen
wir euch ein Makefile zur Verfügung, so dass gegen die Intel MKL gelinkt wird:

Lösung von Session 4 mit Intel MKL Benchmark
--------------------------------------------
:import: session5/page01/gemv.c


Makefile (Angepasst an die Installation im Pool E44)
----------------------------------------------------
:import: session5/page01/Makefile


:navigate: up    -> doc:index
           next  -> doc:session5/page02