====================== Geblockte LU-Zerlegung ====================== Ziel: Implementiert in der Funktion `hpc::ulmlapack::getrf` die geblockte LU-Zerlegung mit partieller Pivotisierung. Idee für den Algorithmus ohne Pivotisierung =========================================== Wir wählen eine bestimmte Blockgröße `bs` (zum Beispiel `bs=64`). In ersten Schritt betrachten wir für eine Matrix $A$ die Partitionierung entlang der ersten `bs` Spalten: ---- LATEX --------------------------------------------------------------------- A = \left(\begin{array}{c|c} A_1 & A_2 \end{array}\right) -------------------------------------------------------------------------------- Für den Block $A_1$ bestimmen wir mit der ungeblockten LU-Zerlegung ---- LATEX --------------------------------------------------------------------- A_1 = L_1 U_1 -------------------------------------------------------------------------------- mit ---- LATEX --------------------------------------------------------------------- L_1 = \left(\begin{array} L_{1,1} \\ \hline L_{2,1} \end{array}\right), \quad L_{1,1} \in \mathbb{R}^{bs \times bs} -------------------------------------------------------------------------------- und ---- LATEX --------------------------------------------------------------------- U_1 = \left(\begin{array} U_{1,1} \\ \hline \mathbf{0} \end{array}\right), \quad U_{1,1} \in \mathbb{R}^{bs \times bs} -------------------------------------------------------------------------------- Analog kann man nun die Partitionierung ---- LATEX --------------------------------------------------------------------- A = \left(\begin{array}{c|c} A_1 & A_2 \end{array}\right) = \left(\begin{array}{c|c} A_{1,1} & A_{1,2} \\ \hline A_{2,1} & A_{2,2} \end{array}\right) \quad\text{mit} A_{1,1} \in \mathbb{R}^{bs \times bs} -------------------------------------------------------------------------------- betrachten. Dann folgt ---- LATEX --------------------------------------------------------------------- A = L U \quad\Leftrightarrow\quad \left(\begin{array}{c|c} A_{1,1} & A_{1,2} \\ \hline A_{2,1} & A_{2,2} \end{array}\right) = \left(\begin{array}{c|c} L_{1,1} & \mathbf{0} \\ \hline L_{2,1} & L_{2,2} \end{array}\right) \left(\begin{array}{c|c} U_{1,1} & U_{1,2} \\ \hline \mathbf{0} & U_{2,2} \end{array}\right) -------------------------------------------------------------------------------- Aufgaben -------- Leitet damit zunächst einen Algorithmus für eine geblockte LU-Zerlegung _ohne_ Pivotisierung her. Dieser soll als Operationen lediglich die ungeblockte LU-Zerlegung, das Matrix-Produkt und das Vorwärtseinsetzen benutzen. Haltet außerdem fest, welche Dimensionen die einzelen Blöcke besitzen. Beachtet dabei, dass $A$ nicht quadratisch sein muss. Die Blockgröße `bs` soll dabei beliebig aber fest sein! Pivotisierung einbauen ====================== Die geblockte LU-Zerlegung soll bei der Zerlegung $A = PLU$ die Permutationsmatrix durch einen Pivotvektor $p$ repräsentieren. Wir betrachten die Partitionierung ---- LATEX --------------------------------------------------------------------- A = \left(\begin{array}{c|c} A_1 & A_2 & A_3 \end{array}\right) = \left(\begin{array}{c|c} A_{1,1} & A_{1,2} & A_{1,3} \\ \hline A_{2,1} & A_{2,2} & A_{2,3} \\ \hline A_{3,1} & A_{3,2} & A_{3,3} \\ \end{array}\right) \quad\text{mit} A_{2,2} \in \mathbb{R}^{bs \times bs} -------------------------------------------------------------------------------- Aufgaben -------- - Wie sind die Dimensionen der einzelnen Blöcke? - Die Blöcke $A_{1,1}$, $A_{2,1}$, $A_{3,1}$, $A_{1,2}$ und $A_{1,3}$ seien bereits mit der LU-Zerlegung überschrieben: - Welche Zeilen könnten dann bereits vertauscht worden sein? - Welcher Teil des Pivotvektors enthält dann schon die endgültigen Werte? - Wie könnte man bei der Restmatrix ---- LATEX ----------------------------------------------------------------- A = \left(\begin{array}{c|c} \overline{A}_1 & \overline{A}_2 \end{array}\right) = \left(\begin{array}{c|c} A_{2,2} & A_{2,3} \\ \hline A_{3,2} & A_{3,3} \\ \end{array}\right) ---------------------------------------------------------------------------- durch Anwenden der ungeblockten LU-Zerlegung auf $\overline{A}_1$ weiter vorgehen? Vorlage für die Implementierung =============================== Als Vorlage für die geblockte LU-Zerlegung `hpc::ulmlapack::getrf` im Verzeichnis `hpc/ulmlapack/` könnt ihr folgendes Gerüst verwenden: ---- CODE (type=cc) ------------------------------------------------------------ #ifndef HPC_ULMLAPACK_GETRF_H #define HPC_ULMLAPACK_GETRF_H 1 namespace hpc { namespace ulmlapack { template Index getrf(Index m, Index n, TA *A, Index incRowA, Index incColA, TP *p, Index incP) { Index mn = std::min(m,n); Index bs = 64; Index info = -1; if (bs<=1 || bs>mn) { /* Matrix zu klein -> ungeblockte Variante aufrufen */ } else { /* Geblockte LU-Zerlegung */ } return info; } } } // namespace ulmblas, hpc #endif // HPC_ULMLAPACK_GETRF_H -------------------------------------------------------------------------------- Test für die geblockte LU-Zerlegung =================================== Um die Implementierung zu testen und einen Benchmark zu erzeugen, kann folgendes Programm benutzt werden: :import: session17/hpc/tests/lu_blocked.cc Aufgaben ======== - Implementiert die Funktion `hpc::ulmlapack::getrf` - Erzeugt Benchmarks, um die Effektivität von `getf2` und `getrf`zu vergleichen. Natürlich sollen diese mit `gnuplot` schön dargestellt werden. - Verwendet auch die Assembler-Kernel für die Benchmarks: - Auf Rechnern mit AVX mit `-DAVX -DD_BLOCKSIZE_MR=4 -DD_BLOCKSIZE_NR=8` und - auf Rechnern mit SSE mit `-DSSE` übersetzen. :navigate: up -> doc:index back -> doc:session17/page04