=============
Dreieckslöser                                                           [TOC]
=============

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.

---- CODE (type=cc) ------------------------------------------------------------
#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.

---- CODE (type=cc) ------------------------------------------------------------
#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.

---- CODE (type=cc) ------------------------------------------------------------
#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.

:import: session17/hpc/tests/sv.cc

:navigate: up    -> doc:index
           back  -> doc:session17/page03
           next  -> doc:session17/page05