============= Dreieckslöser ============= Für eine untere $n\times n$ Dreiecksmatrix $L$ soll mit der Funktion `dtrlsv` (double triangular lower solve vector) die Operation ---- LATEX ----------------------------------------------------------- x \leftarrow L^{-1} x ---------------------------------------------------------------------- durchgeführt werden. Der Vektor $x$ soll also mit der Lösung $L^{-1} x$ überschrieben werden. Bei der Herleitung unterscheiden wir aber formal zwischen dem Vektor vor und nach der Operation und definieren ---- LATEX ----------------------------------------------------------- z := L^{-1} x ---------------------------------------------------------------------- als Lösung des Problems. Dann muss nach entsprechender Partitionierung von $x$ und $z$ gelten: ---- LATEX ----------------------------------------------------------- \left(\begin{array}{c} x_0 \\ \hline x_k \end{array}\right) = \left(\begin{array}{c|c} L_{0,0} & \\ \hline L_{k,0} & L_{k,k} \end{array}\right) \left(\begin{array}{c} z_0 \\ \hline z_k \end{array}\right) = \left(\begin{array}{c} L_{0,0} z_0 \\ \hline L_{k,0} z_0 + L_{k,k} z_k \end{array}\right) \quad \Longleftrightarrow\quad \left\{\begin{array}{lcl} z_0 &=& L_{0,0}^{-1} \; x_0 \\ L_{k,k} z_k &=& x_k - L_{k,0}\, z_0 \end{array}\right. ---------------------------------------------------------------------- Wir leiten daraus zwei Varianten her, mit denen die Operation iterativ mit den Zwischenergebnissen ---- LATEX ----------------------------------------------------------- x = x^{(0)} \leadsto x^{(1)} \leadsto \dots x^{(n)} = z ---------------------------------------------------------------------- berechnet wird: - In der ersten Variante (Row-Major) soll für $x^{(k)}$ gelten: ---- LATEX ----------------------------------------------------------- x^{(k)} = \left(\begin{array}{c} x_0^{(k)} \\ \hline x_k^{(k)} \end{array}\right) = \left(\begin{array}{c} L_{0,0}^{-1} \; x_0 \\ \hline x_k \end{array}\right) = \left(\begin{array}{c} z_0 \\ \hline x_k \end{array}\right) ---------------------------------------------------------------------- - In der zweiten Variante (Col-Major) soll für $x^{(k)}$ gelten: ---- LATEX ----------------------------------------------------------- x^{(k)} = \left(\begin{array}{c} x_0^{(k)} \\ \hline x_k^{(k)} \end{array}\right) = \left(\begin{array}{c} L_{0,0}^{-1} x_0 \\ \hline x_k - L_{k,0}\, z_0 \end{array}\right) = \left(\begin{array}{c} z_0 \\ \hline x_k - L_{k,0}\, z_0 \end{array}\right) ---------------------------------------------------------------------- Dass die Bezeichnungen Row-Major und Col-Major berechtigt sind, wird allerdings erst im Nachhinein deutlich. Row-Major (Unblocked) ===================== Wir betrachten die gemeinsamen Partitionierungen ---- LATEX ----------------------------------------------------------- x^{(k)} = \left(\begin{array}{c} x_0^{(k)} \\ \hline x_k^{(k)} \\ \hdashline x_{k+1}^{(k)} \end{array}\right) \qquad\text{und}\qquad x^{(k+1)} = \left(\begin{array}{c} x_0^{(k+1)} \\ \hdashline x_k^{(k+1)} \\ \hline x_{k+1}^{(k+1)} \end{array}\right) ---------------------------------------------------------------------- sowie ---- LATEX ----------------------------------------------------------- L = \left(\begin{array}{c|c|c} L_{0,0} & & \\ \hline \ell_{k,0}^T & \ell_{k,k} & \\ \hline L_{k+1,0} & \ell_{k+1,k} & L_{k+1,k+1} \\ \end{array}\right) ---------------------------------------------------------------------- Für die Durchführung einer Iteration rechnen wir nun nach: - Laut Voraussetzung gilt für $x^{(k)}$: ---- LATEX ----------------------------------------------------------- L_{0,0} x_0^{(k)} = x_0 \qquad\Longleftrightarrow\qquad x_0^{(k)} = z_0 ---------------------------------------------------------------------- - Für $x^{(k+1)}$ soll nach der Iteration ---- LATEX ----------------------------------------------------------- \begin{array}{lclcl} L_{0,0} x_0^{(k+1)} & & &=& x_0 \\ \ell_{k,0}^T x_0^{(k+1)} &+& \ell_{k,k} x_k^{(k+1)} &=& x_k \end{array} ---------------------------------------------------------------------- gelten. Offensichtlich ist also $x_0^{(k+1)} = x_0^{(k)} = z_0$ und ---- LATEX ----------------------------------------------------------- x_k^{(k+1)} = \frac{x_k - \ell_{k,0}^T x_0^{(k+1)}}{\ell_{k,k}}. ---------------------------------------------------------------------- Um $x$ mit $z$ zu überschreiben, erhalten wir als Verfahren ---- BOX ------------------------------------------------------------- - For $k=0, \dots, n-1$: - $x_k \leftarrow x_k - \ell_{k,0}^T x_0$ (DOT) - $x_k \leftarrow \frac{x_k}{\ell_{k,k}}$ ---------------------------------------------------------------------- Wie es die Prophezeiung besagte, wird hier nur auf Zeilen von $L$ zugegriffen. Col-Major (Unblocked) ===================== - Laut Voraussetzung gilt für $x^{(k)}$: ---- LATEX ----------------------------------------------------------- L_{0,0} x_0^{(k)} = x_0 \qquad\Longleftrightarrow\qquad x_0^{(k)} = z_0 ---------------------------------------------------------------------- und ---- LATEX ----------------------------------------------------------- \left(\begin{array}{c} x_k^{(k)} \\ \hdashline x_{k+1}^{(k)} \end{array}\right) = \left(\begin{array}{c} x_k \\ \hdashline x_{k+1} \end{array}\right) - \left(\begin{array}{c} \ell_{k,0}^T \\ \hdashline L_{k+1,0} \end{array}\right) z_0 \qquad\Longleftrightarrow\qquad \left\{\begin{array}{lclcl} x_k^{(k)} &=& x_k &-& \ell_{k,0}^T z_0 \\ x_{k+1}^{(k)} &=& x_{k+1} &-& L_{k+1,0} z_0 \end{array}\right. ---------------------------------------------------------------------- - Für $x^{(k+1)}$ soll nach der Iteration ---- LATEX ----------------------------------------------------------- \left(\begin{array}{c:c} L_{0,0} & \\ \hdashline \ell_{k,0}^T & \ell_{k,k} \\ \end{array}\right) \left(\begin{array}{c} x_0^{(k+1)} \\ \hdashline x_k^{(k+1)} \end{array}\right) = \left(\begin{array}{c} x_0 \\ \hdashline x_k \end{array}\right) \qquad\Longleftrightarrow\qquad \left\{\begin{array}{lclcl} L_{0,0} x_0^{(k+1)} & & &=& x_0 \\ \ell_{k,0}^T x_0^{(k+1)} &+& \ell_{k,k} x_k^{(k+1)} &=& x_k \end{array}\right. \qquad\Longleftrightarrow\qquad \left\{\begin{array}{lcl} x_0^{(k+1)} &=& z_0 \\ \ell_{k,k} x_k^{(k+1)} &=& x_k - \ell_{k,0}^T z_0 \end{array}\right. ---------------------------------------------------------------------- und ---- LATEX ----------------------------------------------------------- \begin{array}{lcl} x_{k+1}^{(k+1)} &=& x_{k+1} - \left(\begin{array}{c:c} L_{k+1,0} & \ell_{k+1,k} \end{array}\right) \left(\begin{array}{c} x_0^{(k+1)}\\ \hdashline x_k^{(k+1)} \end{array}\right) \\ &=& x_{k+1} - L_{k+1,0} x_0^{(k+1)} - \ell_{k+1,k} x_k^{(k+1)}\\ &=& x_{k+1} - L_{k+1,0} z_0 - \ell_{k+1,k} x_k^{(k+1)}\\ &=& x_{k+1}^{(k)} - \ell_{k+1,k} x_k^{(k+1)}\\ \end{array} ---------------------------------------------------------------------- gelten. Wegen $x_0^{(k+1)} = x_0^{(k)} = z_0$ muss also nur ---- LATEX ----------------------------------------------------------- \begin{array}{lcl} x_k^{(k+1)} &=& \frac{1}{\ell_{k,k}}\left(x_k - \ell_{k,0}^T z_0\right) = \frac{1}{\ell_{k,k}} x_k^{(k)} \\ x_{k+1}^{(k+1)} &=& x_{k+1}^{(k)} - \ell_{k+1,k} x_k^{(k+1)} \end{array} ---------------------------------------------------------------------- berechnet werden. Als Algorithmus erhalten wir somit: ---- BOX ------------------------------------------------------------- - For $k=0, \dots, n-1$: - $x_k \leftarrow \frac{1}{\ell_{k,k}} x_k$ - $x_{k+1} \leftarrow x_{k+1} - \ell_{k+1,k} x_k$ (AXPY) ---------------------------------------------------------------------- Aufgabe ======= Implementiert die Prozeduren `dtrlsv_row` und `dtrlsv_col` in der untenstehenden Vorlage. - Testet eure Implementierung jeweils für zeilen- und spaltenweise gespeicherte Matrizen. Übersetzt dazu wiederum das Programm jeweils mit `-DROWMAJOR=1` und `-DROWMAJOR=0`. - Verwendet verschiedene Optimierungen: `-O1`, `-O3` und `-Ofast`. :import: session06/example01/bench_trlsv.c :navigate: up -> doc:index next -> doc:session06/page02