Content |
Dreieckslöser
Ist die Systemmatrix \(A\) eine Dreiecksmatrix, so kann die Matrixgleichung \(AX=B\) durch Vorwärts- oder Rückwärtseinsetzen gelöst werden. Diese Löser werden ein wichtiger Baustein für die geblockte LU-Zerlegung sein. Zunächst motivieren wir diese aber als Hilfsmittel, um ein Gleichungssystem zu lösen, nachdem für dies Systemmatrix die LU-Zerlegung bestimmt wurde.
-
Die Funktion hpc::ulmblas::trlsm führt für eine untere Dreiecksmatrix \(L\) die Operation \(X \leftarrow \alpha L^{-1} B\) durch. Die Funktion greift bei der \(m \times m\) Matrix \(L\) nur auf Elemente der unteren Dreiecksmatrix zu. Es wird deshalb nicht verlangt, dass im oberen Dreiecksteil tatsächlich Nullen stehen.
-
Die Funktion hpc::ulmblas::trusm führt für eine obere Dreiecksmatrix \(U\) die Operation \(X \leftarrow \alpha U^{-1} B\) durch. Wie bei trlsm wird nur auf Elemente der relevanten oberen Dreiecksmatrix zugegriffen.
Dabei können beide Funktionen den Speziallfall behandeln, dass angenommen werden soll, dass die Diagonal nur Einsen enthält. In diesem Fall wird auch nicht auf Diagonalelemente zugegriffen.
Um ein lineaeres Gleichungssystem \(AX=B\) mit Hilfe der LU-Zerlegung zu lösen, muss die Permuation \(P\) (durch einen Pivot-Vektor dargestellt) auf die Matrix \(B\) angewandt werden. Anschließend wird die Gleichung durch Vorwärts- und Rückwärtseinsetzen gelöst. Offensichtlich kann dabei \(B\) mit der Lösung \(X\) überschrieben werden:
-
\(B \leftarrow P^T B\)
-
\(B \leftarrow L^{-1} B\)
-
\(B \leftarrow U^{-1} B\)
Für die Operation \(B \leftarrow P^T B\) müssen die bei der LU-Zerlegung durchgeführten Zeilenvertauschungen in der gleichen Reihenfolge auf \(B\) angewandt werden. Dies soll mit Hilfe des Pivotvektors durch die Funktion hpc::ulmlapack::swap ermöglicht werden.
Im Folgenden die Signaturen der Funktionen:
hpc::ulmblas::trlsm
Hat unitDiag den Wert true wird beim Vorwärtseinsetzen \(X \leftarrow \alpha L^{-1} B\) angenommen, dass alle Diagonalelemente den Wert \(1\) besitzen.
#ifndef HPC_ULMBLAS_TRLSM_H #define HPC_ULMBLAS_TRLSM_H 1 namespace hpc { namespace ulmblas { template <typename Index, typename Alpha, typename TA, typename TB> void trlsm(Index m, Index n, const Alpha &alpha, bool unitDiag, const TA *A, Index incRowA, Index incColA, TB *B, Index incRowB, Index incColB) { /* ... */ } } } // namespace ulmblas, hpc #endif // HPC_ULMBLAS_TRLSM_H 1
hpc::ulmblas::trusm
Hat unitDiag den Wert true wird beim Rückwärtseinsetzen \(X \leftarrow \alpha U^{-1} B\) angenommen, dass alle Diagonalelemente den Wert \(1\) besitzen.
#ifndef HPC_ULMBLAS_TRUSM_H #define HPC_ULMBLAS_TRUSM_H 1 namespace hpc { namespace ulmblas { template <typename Index, typename Alpha, typename TA, typename TB> void trusm(Index m, Index n, const Alpha &alpha, bool unitDiag, const TA *A, Index incRowA, Index incColA, TB *B, Index incRowB, Index incColB) { /* ... */ } } } // namespace ulmblas, hpc #endif // HPC_ULMBLAS_TRUSM_H 1
hpc::ulmlapack::swap
Der Vektor \(p\) mit Länge \(m\) enthält die Pivotindizes, die bei der LU-Zerlegung verwendet wurden. Für \(k\) mit \(k_1 \leq k\) soll die \(k\)-te Zeile von \(A\) mit der \(p_k\)-ten Zeile vertauscht werden. Dabei soll \(k\) den Bereich aufsteigend durchlaufen.
#ifndef HPC_ULMLAPACK_SWAP_H #define HPC_ULMLAPACK_SWAP_H 1 #include <cassert> #include <hpc/ulmblas/swap.h> namespace hpc { namespace ulmlapack { template <typename Index, typename TA, typename TP> void swap(Index m, Index n, TA *A, Index incRowA, Index incColA, Index k0, Index k1, TP *p, Index incP) { assert(0<=k0); assert(k0<=k1); assert(k1<=m); /* ... */ } } } // namespace ulmlapack, hpc #endif // HPC_ULMLAPACK_SWAP_H
Aufgabe
Implementiert diese Funktionen. Zum Test kann unten stehender Code benutzt werden. Überlegt euch, wie man bei trlsm und trusm durch kleine Änderungen des Test auch Fälle für \(\alpha \neq 1\) testen kann.
#include <cassert> #include <random> #include <hpc/matvec/copy.h> #include <hpc/matvec/gematrix.h> #include <hpc/matvec/mm.h> #include <hpc/matvec/print.h> #include <hpc/ulmblas/trlsm.h> #include <hpc/ulmblas/trusm.h> #include <hpc/ulmlapack/getf2.h> #include <hpc/ulmlapack/swap.h> // // Random initializer for general matrices: real and complex valued // template <typename Index, typename T> void randomInit(Index m, Index n, T *A, Index incRowA, Index incColA) { std::random_device random; std::default_random_engine mt(random()); std::uniform_real_distribution<T> uniform(-100,100); for (Index i=0; i<m; ++i) { for (Index j=0; j<n; ++j) { A[i*incRowA+j*incColA] = uniform(mt); } } } int main() { using namespace hpc::matvec; typedef double T; typedef std::size_t Index; const Index m = 600; const Index n = 800; GeMatrix<T> A(m, m); GeMatrix<T> X(m, n); GeMatrix<T> B(m, n); GeMatrix<Index> p(m, 1); randomInit(m, m, A.data, A.incRow, A.incCol); randomInit(m, n, X.data, X.incRow, X.incCol); // // A*X = B // mm(T(1), A, X, T(0), B); // // Factorize A=P*L*U // std::size_t info = hpc::ulmlapack::getf2(A.numRows, A.numCols, A.data, A.incRow, A.incCol, p.data, A.incRow); std::printf("getf2 returned: info = %ld\n", info); // // Solve P*L*U*X = B // hpc::ulmlapack::swap(B.numRows, B.numCols, B.data, B.incRow, B.incCol, Index(0), B.numRows, p.data, p.incRow); hpc::ulmblas::trlsm(m, n, T(1), true, A.data, A.incRow, A.incCol, B.data, B.incRow, B.incCol); hpc::ulmblas::trusm(m, n, T(1), false, A.data, A.incRow, A.incCol, B.data, B.incRow, B.incCol); // // Compute residual // double res = 0; for (Index j=0; j<n; ++j) { for (Index i=0; i<m; ++i) { res += std::abs(B(i,j)-X(i,j)); } } std::printf("res=%e\n", res); }