Allgemeines Matrix-Matrix Produkt

Für eine \(m \times k\) Matrix \(A\), eine \(k\times n\) Matrix \(B\), eine \(m \times n\) Matrix \(C\) und Skalare \(\alpha\) und \(\beta\) ist die GEMM-Operation (general matrix-matrix product) definiert als

\[C \leftarrow \beta C + \alpha A \, B.\]

Diese Matrix-Matrix Operation kann man auf drei verschiedene Arten auf Matrix-Vektor Operationen zurückführen:

Betrachtet man \(B\) und damit auch zwangsläufig \(C\) spaltenweise gilt
\[\beta C + \alpha A \, B=\beta\left(\begin{array}{c|c|c} c_0 & \dots & c_{n-1} \end{array}\right)+\alpha\left(\begin{array}{c|c|c} A\,b_0 & \dots & A\,b_{n-1} \end{array}\right)\]
In \(C\) muss die Spalte mit Index \(j\) somit überschrieben werden mit
\[c_j \leftarrow \beta c_j + \alpha A b_j\]
Als Algorithmus erhält man dann:
- For \(j=0, \dots, n-1\)
  - \(c_j \leftarrow \beta c_j + \alpha A b_j\) (GEMV)
Betrachtet man \(A\) und damit auch zwangsläufig \(C\) zeilenweise gilt
\[\beta C + \alpha A \, B=\beta\left(\begin{array}{c} c_0^T \\ \hline \vdots \\ \hline c_{m-1}^T \end{array}\right)+\alpha\left(\begin{array}{c} a_0^T B \\ \hline \vdots \\ \hline a_{m-1}^T B \\ \end{array}\right)\]
In \(C\) muss die Zeile mit Index \(i\) somit überschrieben werden mit
\[c_i^T \leftarrow \beta c_i^T + \alpha a_i^T B\]
Als Algorithmus erhält man dann:
- For \(i=0, \dots, m-1\)
  - \(c_i^T \leftarrow \beta c_i^T + \alpha a_i^T B\) (GEMV)
Betrachtet man \(A\) spaltenweise und damit zwangsläufig \(B\) zeilenweise gilt
\[\beta C + \alpha A \, B=\beta\,C+\alpha\,\left(\begin{array}{c|c|c} a_0 & \dots & a_{k-1} \end{array}\right)\left(\begin{array}{c} b_0^T \\ \hline \vdots \\ \hline b_{k-1}^T \end{array}\right)=\beta\,C+\alpha\,a_0 b_0^T+\dots+\alpha\,a_{k-1} b_{k-1}^T\]
Als Algorithmus erhält man dann:
- \(C \leftarrow \beta\,C\) (GESCAL)
- For \(\ell=0, \dots, k-1\)
  - \(C \leftarrow C + \alpha a_\ell b_\ell^T\) (GER)

Übersicht über alle Varianten

Jede der Matrix-Vektor Operationen kann wiederum zeilen- oder spaltenweise ausgeführt werden. Somit gibt es insgesamt sechs verschiedene Varianten.

Varianten mit Matrix-Vektor Produkt

For \(j=0, \dots, n-1\)
- \(c_j \leftarrow \beta c_j + \alpha A b_j\) (GEMV)

GEMV Row-Major: GEMM_jil

Betrachtet man \(A\) und damit auch \(c_j\) zeilenweise erhält man für \(c_j \leftarrow \beta c_j + \alpha A b_j\) zunächst

\[\left(\begin{array}{c} c_{0,j} \\ \hline \vdots \\ \hline c_{m-1,j} \end{array}\right)\leftarrow\beta\left(\begin{array}{c} c_{0,j} \\ \hline \vdots \\ \hline c_{m-1,j} \end{array}\right)+ \alpha\left(\begin{array}{c} a_{0}^T b_j \\ \hline \vdots \\ \hline a_{m-1}^T b_j \\ \end{array}\right)\]

Und somit einen Algorithmus

For \(j=0, \dots, n-1\)
- For \(i=0, \dots, m-1\)
  - \(c_{i,j} \leftarrow \beta c_{i,j} + \alpha a_{i}^T b_{\ell, j}\)

Betrachtet man jetzt noch \(a_{i}^T\) komponentenweise erhält man

For \(j=0, \dots, n-1\)
- For \(i=0, \dots, m-1\)
  - \(c_{i,j} \leftarrow \beta c_{i,j}\)
  - For \(\ell=0, \dots, k-1\)
    - \(c_{i,j} \leftarrow c_{i,j} + \alpha a_{i,\ell} b_{\ell, j}\)

GEMV Col-Major: GEMM_jli

Geht man analog vor, so erhält man den Algorithmus

For \(j=0, \dots, n-1\)
- For \(i=0, \dots, m-1\)
  - \(c_{i,j} \leftarrow \beta c_{i,j}\)
- For \(\ell=0, \dots, k-1\)
  - For \(i=0, \dots, m-1\)
    - \(c_{i,j} \leftarrow c_{i,j} + \alpha a_{i,\ell} b_{\ell, j}\)

wir optimieren diesen zu

For \(j=0, \dots, n-1\)
- For \(\ell=0, \dots, k-1\)
  - For \(i=0, \dots, m-1\)
    - Falls \(\ell = 0 \)
      - \(c_{i,j} \leftarrow \beta c_{i,j}\)
    - \(c_{i,j} \leftarrow c_{i,j} + \alpha a_{i,\ell} b_{\ell, j}\)

wobei siche Unterschiede in der Laufzeit wahrschelinch gar nicht messen lassen. Allerdings wird sich die Idee dahinter bei der geblockten Variante auszahlen.

Varianten mit Vektor-Matrix Produkt

For \(i=0, \dots, m-1\)
- \(c_i \leftarrow \beta c_i + \alpha a_i^T B\) (GEMV)

GEMV Row-Major: GEMM_ilj

For \(i=0, \dots, m-1\)
- For \(\ell=0, \dots, k-1\)
  - For \(j=0, \dots, n-1\)
    - Falls \(\ell = 0 \)
      - \(c_{i,j} \leftarrow \beta c_{i,j}\)
    - \(c_{i,j} \leftarrow c_{i,j} + \alpha a_{i,\ell} b_{\ell, j}\)

GEMV Col-Major: GEMM_ijl

For \(i=0, \dots, m-1\)
- For \(j=0, \dots, n-1\)
  - \(c_{i,j} \leftarrow \beta c_{i,j}\)
  - For \(\ell=0, \dots, k-1\)
    - \(c_{i,j} \leftarrow c_{i,j} + \alpha a_{i,\ell} b_{\ell, j}\)

Varianten mit Rang-1 Update

\(C \leftarrow \beta\,C\) (GESCAL)
For \(\ell=0, \dots, k-1\)
- \(C \leftarrow C + \alpha a_\ell b_\ell^T\) (GER)

GER Row-Major: GEMM_lij

For \(\ell=0, \dots, k-1\)
- For \(i=0, \dots, m-1\)
  - For \(j=0, \dots, n-1\)
    - Falls \(\ell = 0 \)
      - \(c_{i,j} \leftarrow \beta c_{i,j}\)
    - \(c_{i,j} \leftarrow c_{i,j} + \alpha a_{i,\ell} b_{\ell, j}\)

GER Col-Major: GEMM_lji

For \(\ell=0, \dots, k-1\)
- For \(j=0, \dots, n-1\)
  - For \(i=0, \dots, m-1\)
    - Falls \(\ell = 0 \)
      - \(c_{i,j} \leftarrow \beta c_{i,j}\)
    - \(c_{i,j} \leftarrow c_{i,j} + \alpha a_{i,\ell} b_{\ell, j}\)

Aufgaben

Implementiert die GESCAL Operation in `ulmblas/level1.c
Implementiert diese GEMM Varianten in `ulmblas/level3.c

Unter http://www.mathematik.uni-ulm.de/~lehn/hpc_project.tgz findet ihr eine Vorlage:

Mit wget laden.
Mit tar xfvz hpc_project.tgz entpacken.
In ulmblas/level3.c die Funktionen implementieren.

Übersetzen der ulmBLAS Bibliothek

In ulmblas mit make die Bibliothek erzeugen bzw. aktualisieren.

Durchführen der Benchmarks

In bench mit make die verschiedenen Benchmarks erzeugen bzw. aktualisieren:

bench_gemm_mv
bench_gemm_vm
bench_gemm_vv