Home     Back     Up     Next Browse Files

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.

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:

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.

session17/hpc/tests/sv.cc 
#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);
}
Home     Back     Up     Next