========================== Matrix-Vektor Produkt TRMV [TOC] ========================== Bei der BLAS Level 2 Operation TRMV (*tr*iangular *m*atrix *v*ector product) wird eine im Full-Storage-Format gespeicherte Matrix als obere oder untere Dreieicksmatrix interpretiert. Zusätzlich kann angegeben werden, ob bei der Operation angenommen werden soll, dass die Elemente auf der Diagonale alle Eins sind. Das ist eine reine Vereinbarung: - Soll die Matrix als untere Dreiecksmatrix interpretiert werden, dann müssen in der oberen Dreiecksmatrix keine Nullen stehen. Die dort liegenden Werte werden einfach ignoriert. Bei gewissen LAPACK Funktionen wird die obere Dreiecksmatrix auch mit Zwischenergebnissen überschrieben. - Soll angenommen werden, dass auf der Diagonale Einsen stehen, dann gilt gleiches. Auch hier wird nicht vorausgesetzt, dass im Speicher wirklich Einsen stehen. - Analoges gilt für obere Dreiecksmatrizen. Es geht hier also nicht darum eine Dreiecksmatrix möglichst kompakt zu speichern (wie beispielsweise beim sogenanten Packed-Storage-Format). Der Vorteil ist ein schneller Zugriff auf Elemente und Matrix-Blöcke. Also insbesondere auch auf Spalten und Zeilen. Wir werden hier nur die TRMV-Operation für den Fall einer unteren Dreiecksmatrix realisieren (der andere Fall geht dann analog). In der Prozedur `trlmv` (triangular lower matrix vector product) wird die Operation ---- LATEX ------------------------------------------------------------------- x \longleftarrow L(A) x \quad\text{oder} \quad x \longleftarrow L_1(A) x \quad\text{mit} \quad L(A), L_1(A) \in \mathbb{M}^{n \times n} ----------------------------------------------------------------------------- durchgeführt. Durch einen zusätzlichen Parameter kann angegeben werden, dass Diagonalelemente als Eins angenommen werden sollen. Zu beachten ist, dass bei der Operation der Vektor $x$ mit dem Ergebnis überschrieben werden muss. Aufgabe: Col-Major-Variante für TRLMV ===================================== Von $L(A)$ betrachten wir die Partitionierung bezüglich der Zeile und Spalte mit Index $k$: ---- LATEX ------------------------------------------------------------------- \left(\begin{array}{c} \vec{x_0} \\[0.2cm] \hline x_k \\[0.2cm] \hline \vec{x_{k+1}} \end{array}\right) \longleftarrow \left(\begin{array}{c|c|c} L_{0,0} & & \\ \hline {\vec{l_{k,0}}}^T & l_{k,\,k} & \\[0.2cm] \hline L_{k+1,0} & \vec{l_{k+1,\,k}} & L_{k+1,\,k+1} \end{array}\right) \left(\begin{array}{c} \vec{x_0} \\[0.2cm] \hline x_k \\[0.2cm] \hline \vec{x_{k+1}} \end{array}\right) = \left(\begin{array}{c} L_{0,0} \\[0.2cm] \hline {\vec{l_{k,0}}}^T \\[0.2cm] \hline L_{k+1,0} \end{array}\right) \vec{x_0} + \left(\begin{array}{c} \\[0.2cm] \hline l_{k,\,k} \\[0.2cm] \hline \vec{l_{k+1,\,k}} \end{array}\right) x_k + \left(\begin{array}{c} \\[0.2cm] \hline \\[0.2cm] \hline L_{k+1,\,k+1} \end{array}\right) \vec{x_{k+1}} ----------------------------------------------------------------------------- Dies kann man in die drei Teil-Operationen ---- LATEX ------------------------------------------------------------------ \begin{array}{|c|c|c|} \hline \\ \text{pre}(k) & \text{inv}(k) & \text{post}(k) \\[0.8cm] \hline \\ \quad \vec{x_{k+1}} \longleftarrow L_{k+1,\,k+1} \vec{x_{k+1}} \quad & \quad \left\{ \begin{array}{lccl} (1) & \vec{x_{k+1}} & \leftarrow & \vec{x_{k+1}} + \vec{l_{k+1,\,\,k}} x_k\;, \\[0.5cm] (2) & x_k & \leftarrow & l_{k,\,k} x_k\;. \end{array} \right. \quad & \quad \left\{ \begin{array}{lccl} (1) & \vec{x_{k+1}} & \leftarrow & \vec{x_{k+1}} + L_{k+1,0} \vec{x_0}\;, \\[0.5cm] (2) & x_k & \leftarrow & x_k + {\vec{l_{k,0}}}^T \vec{x_0}\;, \\[0.5cm] (3) & \vec{x_0} & \leftarrow & L_{0,0} \vec{x_0}\;. \end{array} \right. \quad \\[1cm] \\ \hline \end{array} ----------------------------------------------------------------------------- zerlegen. Da $x$ überschrieben wird müssen diese in der Reihenfolge $\text{pre}(k)$, $\text{inv}(k)$, $\text{post}(k)$ durchgeführt werden. Wir erhalten damit folgenden iterativen Algorithmus: ---- LATEX ------------------------------------------------------------------- \begin{array}{l} \text{for}\; k=n-1, \dots, 0 \\ \qquad \begin{array}{lcl} \vec{x_{k+1}} & \longleftarrow & \vec{x_{k+1}} + \vec{l_{k+1,\,k}} \cdot x_k \\ x_k & \longleftarrow & l_{k,\,k} \cdot x_k \\ \end{array} \\ \\ \end{array} ----------------------------------------------------------------------------- Dabei muss aber zusätzlich der Fall behandelt werden, dass eventuell $l_{k,\,k}=1$ gelten soll. Implementiert diesen Algorithmus mit der Prozedur ---- CODE(type=c) -------------------------------------------------------------- void dtrlmv_col(int n, int unitDiag, const double *A, int incRowA, int incColA, double *x, int incX); -------------------------------------------------------------------------------- Im Fall $\text{unitDiag} \neq 0$ soll angenommen werden, dass alle Diagonalelemente den Wert Eins besitzen. Aufgabe: Row-Major-Variante für TRLMV ===================================== Wir betrachten wieder die gleiche Partionierung von $L(A)$. Beim Ausmultiplizieren fassen wir das Ergebnis aber zeilenweise zusammen: ---- LATEX ------------------------------------------------------------------- \left(\begin{array}{c} \vec{x_0} \\[0.2cm] \hline x_k \\[0.2cm] \hline \vec{x_{k+1}} \end{array}\right) \longleftarrow \left(\begin{array}{c|c|c} L_{0,0} & & \\[0.2cm] \hline {\vec{l_{k,0}}}^T & l_{k,\,k} & \\[0.2cm] \hline L_{k+1,0} & \vec{l_{k+1,\,k}} & L_{k+1,\,k+1} \end{array}\right) \left(\begin{array}{c} \vec{x_0} \\[0.2cm] \hline x_k \\[0.2cm] \hline \vec{x_{k+1}} \end{array}\right) = \left(\begin{array}{lclcl} L_{0,0} \,\vec{x_0} && &&\\[0.2cm] \hline {\vec{l_{k,0}}}^T \,\vec{x_0} &+& l_{k,\,k} \,x_k & & \\[0.2cm] \hline L_{k+1,0} \,\vec{x_0} &+& l_{k+1,\,k} \,x_k &+& L_{k+1,\,k+1} \,\vec{x_{k+1}} \end{array}\right) ----------------------------------------------------------------------------- Das ist also äquivalent zu den drei Teil-Operationen ---- LATEX ------------------------------------------------------------------ \begin{array}{|c|c|c|} \hline \\ \text{pre}(k) & \text{inv}(k) & \text{post}(k) \\[0.8cm] \hline \\ \qquad \vec{x_{k+1}} \longleftarrow L_{k+1,0} \,\vec{x_0} + l_{k+1,\,k} \,x_k + L_{k+1,\,k+1} \,\vec{x_{k+1}} \qquad & \qquad x_k \longleftarrow {\vec{l_{k,0}}}^T \,\vec{x_0} + l_{k,\,k} \,x_k \qquad & \qquad \vec{x_0} \longleftarrow L_{0,0} \,\vec{x_0} \qquad \\[1cm] \\ \hline \end{array} ----------------------------------------------------------------------------- Wir erhalten damit folgenden iterativen Algorithmus: ---- LATEX ------------------------------------------------------------------- \begin{array}{l} \text{for}\; k=n-1, \dots, 0 \\ \qquad \begin{array}{lcl} x_k & \longleftarrow & {\vec{l_{k,0}}}^T \,\vec{x_0} + l_{k,\,k} \,x_k \end{array}\\ \\ \end{array} ----------------------------------------------------------------------------- Auch hier muss zusätzlicher der eventuelle Fall $l_{k,\,k}=1$ gesondert behandelt werden. Implementiert diesen Algorithmus mit der Prozedur ---- CODE(type=c) -------------------------------------------------------------- void dtrlmv_row(int n, int unitDiag, const double *A, int incRowA, int incColA, double *x, int incX); -------------------------------------------------------------------------------- Grundgerüst für Test ==================== :import: session05/example04/trmv_test.c Aufgabe: Geblockte TRMV Varinaten ================================= Die Operation soll nun bezüglich einem Blockfaktor $\text{bs}$ quasi-äquidistant partitioniert werden. Bezüglich der Blockzeile und -spalte erhält man damit ---- LATEX ------------------------------------------------------------------- \left(\begin{array}{c} \vec{x_0} \\[0.2cm] \hline \vec{x_k} \\[0.2cm] \hline \vec{x_{k+\text{bs}}} \end{array}\right) \longleftarrow \left(\begin{array}{l|l|l} L_{0,0} & & \\[0.2cm] \hline L_{k,0} & L_{k,\,k} & \\[0.2cm] \hline L_{k+\text{bs},0} & L_{k+\text{bs},\,k} & L_{k+\text{bs},\,k+\text{bs}} \end{array}\right) \left(\begin{array}{c} \vec{x_0} \\[0.2cm] \hline \vec{x_k} \\[0.2cm] \hline \vec{x_{k+\text{bs}}} \end{array}\right) = \left(\begin{array}{l} L_{0,0} \\[0.2cm] \hline L_{k,0} \\[0.2cm] \hline L_{k+\text{bs},0} \\ \end{array}\right) ----------------------------------------------------------------------------- Col-Major Algorithmus --------------------- Für einen Block-Spalten orientierten Algorithmus betrachtet man ---- LATEX ------------------------------------------------------------------- \left(\begin{array}{l|l|l} L_{0,0} & & \\[0.2cm] \hline L_{k,0} & L_{k,\,k} & \\[0.2cm] \hline L_{k+\text{bs},0} & L_{k+\text{bs},\,k} & L_{k+\text{bs},\,k+\text{bs}} \end{array}\right) \left(\begin{array}{c} \vec{x_0} \\[0.2cm] \hline \vec{x_k} \\[0.2cm] \hline \vec{x_{k+\text{bs}}} \end{array}\right) = \left(\begin{array}{l} L_{0,0} \\[0.2cm] \hline L_{k,0} \\[0.2cm] \hline L_{k+\text{bs},0} \\ \end{array}\right) \vec{x_0} + \left(\begin{array}{l} \\[0.2cm] \hline L_{k,\,k} \\[0.2cm] \hline L_{k+\text{bs},\,k} \\ \end{array}\right) \vec{x_k} + \left(\begin{array}{l} \\[0.2cm] \hline \\[0.2cm] \hline L_{k+\text{bs},\,k+\text{bs}} \\ \end{array}\right) \vec{x_{k+\text{bs}}} ----------------------------------------------------------------------------- Row-Major Algorithmus --------------------- Für einen Block-Zeilen orientierten Algorithmus betrachtet man ---- LATEX ------------------------------------------------------------------- \left(\begin{array}{l|l|l} L_{0,0} & & \\[0.2cm] \hline L_{k,0} & L_{k,\,k} & \\[0.2cm] \hline L_{k+\text{bs},0} & L_{k+\text{bs},\,k} & L_{k+\text{bs},\,k+\text{bs}} \end{array}\right) \left(\begin{array}{c} \vec{x_0} \\[0.2cm] \hline \vec{x_k} \\[0.2cm] \hline \vec{x_{k+\text{bs}}} \end{array}\right) = \left(\begin{array}{lclcl} L_{0,0} \;\vec{x_0} && && \\[0.2cm] \hline L_{k,0} \;\vec{x_0} &+& L_{k,\,k} \;\vec{x_k} \\[0.2cm] \hline L_{k+\text{bs},0} \;\vec{x_0} &+& L_{k+\text{bs},\,k} \;\vec{x_k} &+& L_{k+\text{bs},\,k+\text{bs}} \;\vec{x_{k+\text{bs}}} \\ \end{array}\right) ----------------------------------------------------------------------------- Leitet für Col-Major und Row-Major Matrizen jeweils einen geblockten Algorithmus für die TRLMV-Operation her. Diese sollen die speziellen GEMV-Operationen für Matrix-Paneele und die ungeblockten TRLMV-Operationen verwenden. Implementiert die Algorithmen in einer gemeinsamen Prozedur: ---- CODE(type=c) -------------------------------------------------------------- void dtrlmv_ulm(int n, int unitDiag, const double *A, int incRowA, int incColA, double *x, int incX); -------------------------------------------------------------------------------- Durch eine Fallunterscheidung soll stets der günstigere Algorithmus verwendet werden. Grundgerüst für Benchmark ------------------------- :import: session05/example04/trmv_bench.c :navigate: up -> doc:index back -> doc:session05/page05 next -> doc:session05/page07