========================== Matrix-Vektor Produkt GEMV [TOC] ========================== In BLAS Level 2 wird mit der GEMV-Operation (*ge*neral *m*atrix *v*ector product) ---- LATEX --------------------------------------------------------------------- y \leftarrow \beta y + \alpha A x, \quad A \in \mathbb{M}^{m \times n}, \; x \in \mathbb{M}^{n}, \; y \in \mathbb{M}^{m} \quad\text{und}\quad \alpha,\; \beta \in \mathbb{M} -------------------------------------------------------------------------------- das Matrix-Vektor Produkt für eine General-Matrix $A$ definiert. Wir stellen zwei algorithmische Varianten vor, die anschliessend implementiert werden sollen. Col-Major-Variante für GEMV =========================== Die $m \times n$ Matrix $A$ wird mit ---- LATEX --------------------------------------------------------------------- A = \left(a_1, \dots, a_n \right) -------------------------------------------------------------------------------- spaltenweise betrachtet. Für eine konsitente Partitionierung des Produktes wird entsprechend $x$ komponentenweise als ---- LATEX --------------------------------------------------------------------- x = \left( \begin{array}{c} x_1 \\ \vdots \\ x_n \end{array} \right) -------------------------------------------------------------------------------- betrachtet. Die Berechnung soll durch ---- LATEX --------------------------------------------------------------------- \begin{array}{lcl} y &\leftarrow& \beta y \\ y &\leftarrow& y + \alpha\, a_1\, x_1\\ &\vdots & \\ y &\leftarrow& y + \alpha\,a_n\, x_n\\ \end{array} -------------------------------------------------------------------------------- erfolgen. Dazu sollen die BLAS-Level 1 Funktionen für die SCAL und AXPY Operationen benutzt werden. Als Signatur ist ---- CODE(type=c) -------------------------------------------------------------- void dgemv_col(int m, int n, double alpha, const double *A, int incRowA, int incColA, const double *x, int incX, double beta, double *y, int incY); -------------------------------------------------------------------------------- zu verwenden. Row-Major-Variante für GEMV =========================== Die $m \times n$ Matrix $A$ wird mit ---- LATEX --------------------------------------------------------------------- A = \left( \begin{array}{c} a_1^T \\ \vdots \\ a_m^T \end{array} \right) -------------------------------------------------------------------------------- zeilenweise betrachtet. Für eine konsitente Partitionierung des Produktes wird entsprechend $y$ komponentenweise als ---- LATEX --------------------------------------------------------------------- y = \left( \begin{array}{c} y_1 \\ \vdots \\ y_m \end{array} \right) -------------------------------------------------------------------------------- betrachtet. Die Berechnung soll durch ---- LATEX --------------------------------------------------------------------- \begin{array}{lcl} y_1 &\leftarrow& \beta\, y_1 + \alpha\,a_1^T\, x \\ &\vdots & \\ y_m &\leftarrow& \beta\, y_m + \alpha\,a_m^T\, x \\ \end{array} -------------------------------------------------------------------------------- erfolgen. Dazu soll die BLAS-Level 1 Funktion DOT für das Skalarprodukt verwendet werden. Als Signatur ist ---- CODE(type=c) -------------------------------------------------------------- void dgemv_row(int m, int n, double alpha, const double *A, int incRowA, int incColA, const double *x, int incX, double beta, double *y, int incY); -------------------------------------------------------------------------------- zu verwenden. Test nach Augenmaß ================== Testet eure Implementation zunächst an einfachen Beispielen, die man einfach nachrechnen kann. Ein etwas rigoroserer Test erfolgt dann bei den Benchmarks. Ihr könnt nachstehendes Grundgerüst verwenden. Natürlich könnt ihr die BLAS Level 1 Funktionen und Prozeduren (ddot, dscal,...) durch eure eigenen Implementierungen ersetzen. :import: session05/example01/gemv_test.c Test und Benchmark ================== Der Benchmark vergleicht beide Implementierungen mit einer Referenz-Implementierung: - Die Fehlerschranke wird berechnet durch ---- LATEX ------------------------------------------------------------------- \text{err} = \frac{ \|y_1 - y_0 \|_\infty }{ \max\{m, n\} \, |\alpha| \, |\beta|\,\|A\|_1 \, \|x\|_\infty \|y_0\|_\infty } ------------------------------------------------------------------------------ wobei mit $y_0$ das Ergebnis der Referenzimplemnetierung und mit $y_1$ das Ergebnis einer zu testenden Implementierung bezeichnet wird. - Die Rechenleistung (Rechenoperationen pro Zeit) berechnet sich bei der GEMV-Operation als ---- LATEX ------------------------------------------------------------------- \text{FLOPS} = \frac{\text{#OP}_{\text{add}} + \text{#OP}_{\text{mult}}}{t} = \frac{m \cdot n + m\cdot (n+1)}{t} ------------------------------------------------------------------------------ und ---- LATEX ------------------------------------------------------------------- \text{MFLOPS} = \frac{\text{#OP}_{\text{add}} + \text{#OP}_{\text{mult}}}{t} 10^{-6} ------------------------------------------------------------------------------ Das entspricht genau genommen nur dem Fall $\beta = 1$. Außerdem wird hier stets die Mindestanzahl der Operationen benutzt. Die Implementierung von `dgemv_col` führt bei einem Aufruf mehr Multiplikationen aus. Nach dem von-Neumann-Modell müsste diese Implementierung also immer schlechter sein als `dgemv_row`. Aufgaben -------- Fügt in dem Grundgerüst eure Implementierung für `dgemv_col` und `dgemv_row` ein und führt Tests für verschiedene Parameter und Compiler-Optimierungen durch: - Führt jeden Benchmark einmal für Row-Major (`-DROWMAJOR=1`) und Col-Major (`-DROWMAJOR=0`) Matrizen aus. - Übersetzt jeweils mit den Optimierungsflags `-O1`, `-O2`, `-O3` und `-Ofast`. *Achtet bei der Ausgabe zuerst darauf, ob eure Implementierung korrekt ist und in den `err`-Spalten Werte im Bereich der Maschinengenauigkeit stehen. Erst dann auf die MFLOPS schauen!* Wenn alles klaptt, dann sollen die Ergebnisse in hübschen Plots aufbereitet werden. Folgendes Vorgehen kann man später leicht durch Makefiles automatisieren: - Die ausführbaren Dateien bekommen einen Namensuffix `_row`, falls mit `-DROWMAJOR=1` übersetzt wird, und sonst den Zusatz `_col`. Ein weiteres Suffix `_O1`, ... `_O3` oder `_Ofast` markiert die verwendete Optimierung. Insgesamt also übersetzen wir mit ---- CODE(type=txt) ---------------------------------------------------------- gcc -Wall -O1 -DROWMAJOR=1 -o gemv_bench_row_O1 gemv_bench.c gcc -Wall -O2 -DROWMAJOR=1 -o gemv_bench_row_O2 gemv_bench.c gcc -Wall -O3 -DROWMAJOR=1 -o gemv_bench_row_O3 gemv_bench.c gcc -Wall -Ofast -DROWMAJOR=1 -o gemv_bench_row_Ofast gemv_bench.c gcc -Wall -O1 -DROWMAJOR=0 -o gemv_bench_col_O1 gemv_bench.c gcc -Wall -O2 -DROWMAJOR=0 -o gemv_bench_col_O2 gemv_bench.c gcc -Wall -O3 -DROWMAJOR=0 -o gemv_bench_col_O3 gemv_bench.c gcc -Wall -Ofast -DROWMAJOR=0 -o gemv_bench_col_Ofast gemv_bench.c ------------------------------------------------------------------------------ - Die Ausgaben lenken wir anschließend in Report-Dateien mit jeweils entsprechendem Suffix um: ---- CODE(type=txt) ---------------------------------------------------------- gemv_bench_row_O1 > report.gemv.row_O1 gemv_bench_row_O2 > report.gemv.row_O2 gemv_bench_row_O3 > report.gemv.row_O3 gemv_bench_row_Ofast > report.gemv.row_Ofast gemv_bench_col_O1 > report.gemv.col_O1 gemv_bench_col_O2 > report.gemv.col_O2 gemv_bench_col_O3 > report.gemv.col_O3 gemv_bench_col_Ofast > report.gemv.col_Ofast ------------------------------------------------------------------------------ - Mit zwei GNUPLOT-Skripten beschreiben wir, wie daraus Plots erzeugt werden sollen: - Plots für RowMajor-Benchmarks: :import: session05/example01/plot.gemv_row - Plots für ColMajor-Benchmarks: :import: session05/example01/plot.gemv_col Erzeugt werden die Plots dann durch ---- CODE(type=txt) ---------------------------------------------------------- gnuplot plot.gemv_row gnuplot plot.gemv_col ------------------------------------------------------------------------------ Grundgerüst ----------- :import: session05/example01/gemv_bench.c :navigate: up -> doc:index next -> doc:session05/page02