===================================
The Art of Designing Matrix Classes
===================================

==== BOX ====================================
Natürlich ist die Teilnahme an den Übungen
freiwillig.  Alternativ zum Votieren im E44
kann auch das _submit_-Kommando benutzt
werden.  Hinweise dazu auf der __Seite 5__.
=============================================


Wir werfen zunächst einen kritischen Blick auf die low-level BLAS Funktionen im
Namensraum `hpc::ulmblas`.  Hier werden Matrizen ähnlich wie in C behandelt:

- Eine Matrix wird durch die Dimensionen, einen Zeiger auf das erste Element und
  den Inkrementen beschrieben:

  ==== CODE (type=cc) ==========================================================
  typename <typename Index, typename Alpha, typename TX>
  void
  gescal(Index m, Index n, const Alpha &alpha,
         TX *X, Index incRowX, Index incColX)
  {
      /* ...*/
  }
  ==============================================================================

  Funktionen haben deshlab eine recht große Anzahl an Parametern.

- Der Zugriff auf Elemente erfolgt durch Pointer Arithmetik:

  ==== CODE (type=cc) ==========================================================
  /* ... */

  X[i*incRowX+j*incColX] *= alpha;

  /* ... */
  ==============================================================================

  Dies verlangt ein Verständnis dafür, wie Matrizen im Speicher liegen.  Der
  Programmierer muss sich zudem sicher sein, dass nur auf zulässige Elemente
  mit $0\leq i < m$ und $0\leq j < n$ zugegriffen wird.  Copy-Paste-Fehler, die
  zum Beispiel zu Fehlern wie `incRowX` statt `incColX` führen, können nur durch
  Sorgfalt vermieden oder aufgespürt werden.

- Die Speicherverwaltung liegt im Verantwortungsbereich des Programmierer

  ==== CODE (type=cc) ==========================================================
  void
  foo()
  {
      double *A = new double[5*3];

      /* ... */

      hpc::ulmblas::gescal(5, 3, 1,5, A, 1, 5);

      /* ... */

      delete [] A;   // nicht vergessen!
  }
  ==============================================================================

Diese Zusammenfassung soll allerdings keine generelle Kritik darstellen.  Ein
Vorteil ist, dass nur eine einfache Teilmenge der C++ Programmiersprache
verwendet wird.  Im Prinzip sind Programmierkenntnisse in C ausreichend.

Für die Implementierung von numerischen Verfahren kann es aber vorteilhaft
sein, wenn diese technischen Aspekte abstrahiert werden.  Dazu wurden im
Namensraum `hpc::matvec` Matrix-Klassen eingeführt.  Die wesentlichen
Eigenschaften dieser Klassen fassen wir zunächst zusammen.  Anschliessend
überlegen wir uns, ob es möglich wäre, dies einfacher zu realisieren.


Ein Blick auf das Vorlesungsprojekt
===================================
Im Vorlesungsprojekt der letzten Session wurden für eine rechteckige Matrix
(general matrix with full storage) drei verschiedene Klassen implementiert:

- Eine Instanz/Objekt der Klasse `GeMatrix` verwaltet einen zugehörigen
  Speicherbereich.  Im Konstruktor wird dieser allokiert und im Destruktor
  wieder freigegeben.

  ==== CODE (type=cc) ==========================================================
  GeMatrix<double>   A(5, 6);   // Allokiert Speicher fuer eine 5x6 Matrix
  ==============================================================================

- Eine Instanz der Klasse `GeMatrixView` verwendet lediglich einen bereits
  allokierten Speicherbereich.  Dieser kann zum Beispiel von einer `GeMatrix`
  allokiert worden sein:

  ==== CODE (type=cc) ==========================================================
  // Matrix A allokiert Speicher fuer eine 5x6 Matrix
  GeMatrix<double>       A(5, 6);

  // Matrix B referenziert einen Teil der Matrix A:  2 Zeilen und 3 Spalten
  // beginnend bei A.data (der Adresse von A(0,0)).
  GeMatrixView<double>   B(2, 3, A.data, A.incRow, A.incCol);

  // Analog haette man diese View wie folgt erzeugen koennen:
  GeMatrixView<double>   C(2, 3, &A(0,0), A.incRow, A.incCol);
  ==============================================================================

  Ebenso kann aber auch ein beliebiger Speicherbereich verwendet werden.  Es
  muss lediglich sichergestellt sein, dass dieser große genug ist:

  ==== CODE (type=cc) ==========================================================
  double A[5*6];                   // Array auf Stack oder Datensegment

  /* ... */

  double B = new double[5*6];      // Array auf Heap


  // Views die den Speicher der obige Arrays verwenden
  GeMatrixView<double>   C(2, 3, A, 2, 3);
  GeMatrixView<double>   D(2, 3, B, 2, 3);
  ==============================================================================

  Der häufigste Anwendungsfall für Views wird aber das Referenzieren von
  Matrix-Blöcken darstellen.  Um dies zu vereinfachen wurde in allen
  Matrix-Klassen der Funktionsoperator überladen.  Zum Beispiel in `GeMatrix`:

  ==== CODE (type=cc) ==========================================================
  template <typename T, typename I=std::size_t>
  struct GeMatrix
  {
      /* ... */

      typedef GeMatrixView<T,Index>       View;

      /* ... */

      View
      operator()(Index i, Index j, Index m, Index n)
      {
          assert(i+m<=numRows);
          assert(j+n<=numCols);
          return View(m, n, &(operator()(i,j)), incRow, incCol);
      }

      /* ... */
  };
  ==============================================================================

  Die Implementierung der `GeMatrixView` Klasse unterscheidet sich von der
  `GeMatrix` Klasse im wesentlichen nur bei den Konstruktoren und Destruktoren:

  - In den Konstruktoren wird kein Speicher allokiert.  Statt dessen wird der
    Daten-Zeiger auf bereits allokierten Speicher gesetzt.

  - Der Destruktor tut *nichts*: Es wurde kein Speicher allokiert und deshalb
    muss auch keiner freigegeben werden.

- Mit `GeMatrixConstView` wurde eine dritte Klasse für Matrizen implementiert.
  Diese gleicht im wesentlichen wiederum der Klasse `GeMatrix` View.
  Allerdings wurden alle Methoden entfernt, die es ermöglichen wurde die Matrix
  zu verändern.  Zusätzlich wurde der Daten-Zeiger als `const` definiert.


Eventuelle Nachteile dieser Realisierung
========================================
Für einen numerischen Algorithmus soll es keine Rolle spielen, ob eine Matrix
ihren eigenen Speicher verwaltet oder nicht.  Um an dieser Stelle nicht auf
numerische Fragestellungen eingehen zu müssen, betrachten wir aber nur eine
Funktion `print` zur Ausgabe einer Matrix.  Wir benötigen bei unserem Ansatz
drei Versionen, die sich nur in der Signatur unterscheiden:

- Ausgabe einer `GeMatrix`:

  ==== CODE (type=cc) ==========================================================
  template <typename Index, typename T>
  void
  print(const GeMatrix<T, Index> &A, const char *name = "")
  {
      if (*name) {
          printf("%s = \n", name);
      }
      for (Index i=0; i<A.numRows; ++i) {
          for (Index j=0; j<A.numCols; ++j) {
              print_value(A(i,j));
          }
          printf("\n");
      }
      printf("\n");
  }
  ==============================================================================

- Ausgabe einer `GeMatrixView`:

  ==== CODE (type=cc) ==========================================================
  template <typename Index, typename T>
  void
  print(const GeMatrixView<T, Index> &A, const char *name = "")
  {
      if (*name) {
          printf("%s = \n", name);
      }
      for (Index i=0; i<A.numRows; ++i) {
          for (Index j=0; j<A.numCols; ++j) {
              print_value(A(i,j));
          }
          printf("\n");
      }
      printf("\n");
  }
  ==============================================================================

- Ausgabe einer `GeMatrixView`:

  ==== CODE (type=cc) ==========================================================
  template <typename Index, typename T>
  void
  print(const GeMatrixConstView<T, Index> &A, const char *name = "")
  {
      if (*name) {
          printf("%s = \n", name);
      }
      for (Index i=0; i<A.numRows; ++i) {
          for (Index j=0; j<A.numCols; ++j) {
              print_value(A(i,j));
          }
          printf("\n");
      }
      printf("\n");
  }
  ==============================================================================

Für Funktionen mit weiteren Matrix-Argumenten explodiert die Anzahl der
Varianten.  Störend und fehleranfällig ist zudem die Implementierung im Stil von
Copy-Paste.  Und noch schlimmer, ändert man anschliessend eine Variante, so muss
man analog alle anderen Varianten ändern.

Verwendung von Template-Funktionen
----------------------------------
Diese stupide Arbeit können wir durch Templates dem Compiler überlassen:

==== CODE (type=cc) ==========================================================
template <typename Matrix>
void
print(const Matrix &A, const char *name = "")
{
    typedef typename Matrix::Index      Index;

    if (*name) {
        printf("%s = \n", name);
    }
    for (Index i=0; i<A.numRows; ++i) {
        for (Index j=0; j<A.numCols; ++j) {
            print_value(A(i,j));
        }
        printf("\n");
    }
    printf("\n");
}
==============================================================================

Dies funktioniert wunderbar, solange man keine `print` Funktion für andere
Datentypen im gleichen Namensraum definieren möchte.  Stellen wir uns vor, dass
wir irgendwann auch für Vektoren entsprechende Klassen implementieren.  Dann
wäre es naheliegend eine zugehörige `print` Funktion zu implementieren:

==== CODE (type=cc) ==========================================================
template <typename Vector>
void
print(const Vector &x, const char *name = "")
{
    typedef typename Vector::Index      Index;

    if (*name) {
        printf("%s = \n", name);
    }
    for (Index i=0; i<x.length; ++i) {
        print_value(x(i));
    }
    printf("\n");
}
==============================================================================


Bei einem Aufruf wie in

==== CODE (type=cc) ==========================================================
void
foo()
{
    GeMatrix<double>  A(4,5);
    /*  ... */
    print(A);
}
==============================================================================

wäre für den Compiler nicht entscheidbar, welcher Funktionsaufruf erfolgen
soll.  Ein Lösung wäre es den Typ der Objekte in den Funktionsnamen
aufzunehmen.  Dies würde zu `printGeMatrix` und `printVector` führen.  Bei
Funktionen mit mehreren Matrix/Vector Parametern würde dies im Allgemeinen zu
sehr langen Funktionsnamen führen.  Auch diese stupide Erzeugung von
unterschiedlichen Funktionsnamen können wir dem Compiler überlassen.  Um
festzulegen welche Funktion aufgerufen werden soll verwenden wir das SFINAE
Konzept mit entsprechenden Traits:

==== CODE (type=cc) ==========================================================
template <typename MA>
typename std::enable_if<IsGeMatrix<MA>::value,
         void>::type
print(const MA &A, const char *name = "")
{
    typedef typename MA::Index    Index;

    if (*name) {
        printf("%s = \n", name);
    }
    for (Index i=0; i<A.numRows; ++i) {
        for (Index j=0; j<A.numCols; ++j) {
            print_value(A(i,j));
        }
        printf("\n");
    }
    printf("\n");
}
==============================================================================

Für jede Matrix/Vector Klasse muss dann aber ein entsprechender Trait wie hier
`IsGeMatrix` implementiert werden:

==== CODE (type=cc) ==========================================================
template <typename Any>
struct IsGeMatrix
{
    static constexpr bool value = false;
};

/* ... */
==============================================================================


Aufgaben
========
Implementiert zwei Template-Klassen `Foo` und `Dummy` und zugehörige Traits,
so dass dieses Beispiel


==== CODE (type=cc) ==========================================================
#include <cstdio>
#include <type_traits>


/* ... */


template <typename T>
typename std::enable_if<IsFoo<T>::value,
         void>::type
print(const T &)
{
        std::printf("foo\n");
}

template <typename T>
typename std::enable_if<IsDummy<T>::value,
         void>::type
print(const T &)
{
        std::printf("dummy\n");
}

int
main()
{
    Foo<double>    foo;
    Dummy<float>   dummy;

    print(foo);
    print(dummy);
}
==============================================================================

Diese Ausgabe erzeugt:

==== CODE (type=tyt) =========================================================
foo
dummy
==============================================================================

:links: Seite 5 -> doc:session14/page05

:navigate: up    -> doc:index
           next  -> doc:session14/page02