BLAS mit gemischen Datentypen

Bisher mussten bei BLAS Funktionen alle Skalare, Vektor- und Matrixelemente den gleichen Datentyp besitzen (z.B. alle den Typ double). Mittels Template Funktionen können wir dies verallgemeinern, solange dies mathematisch definiert ist. Dazu betrachten wir als Beispiel die Operation gecopy zum Kopieren einer Matrix:

\[B \leftarrow A\]

Sind die Elemente von A und B jeweils vom Typ float und double dann kann die Operation durch einen elementweisen upcast durchgeführt werden. Dabei verliert man keine Genauigkeit bei den Fließkommazahlen. In C/C++ wird deshalb bei
```
        float  x;
        double y;
        
        // ...
        
        y = x;  // implicit upcast
        
```
der Upcast automatisch durchgeführt.

Sind die Elemente von A und B jeweils vom Typ double und float dann muss durch Runden elementweise ein downcast durchgeführt werden. Bei einem Code wie

session12/page01/downcast.cc

        int
        main()
        {
            double  x = 1.2;
            float   y;
        
            y = x;     // implicit downcast: double to float
            return y;  // implicit downcast: float to int
        }
        

erhält man beim Übersetzen eine Warnung. Zumindest dann, wenn man beim g++ Compiler entsprechende Warnungen aktiviert:

$shell> g++ -Wall downcast.cc
$shell> g++ -Wall -Wconversion downcast.cc
downcast.cc: In function 'int main()':
downcast.cc:7:7: warning: conversion to 'float' from 'double' may alter its value [-Wfloat-conversion]
     y = x;     // implicit downcast: double to float
       ^
downcast.cc:8:12: warning: conversion to 'int' from 'float' may alter its value [-Wfloat-conversion]
     return y;  // implicit downcast: float to int
            ^

So viel zum Thema -Wall bedeutet, dass alle Warnungen ausgegeben werden sollen.

Eine Konvertierung wie beispielsweis von std::complex<double> nach double ist natürlich im Allgemeinen nicht möglich:

session12/page01/downcast2.cc

        #include <complex>
        
        int
        main()
        {
            std::complex<double>  z(1,2);
            double                y;
        
            y = z;     // error
            return y;  // implicit downcast: float to int
        }
        

Entsprechend erhält man auch einen Fehler:

$shell> g++ downcast2.cc
downcast2.cc: In function 'int main()':
downcast2.cc:9:7: error: cannot convert 'std::complex' to 'double' in assignment
     y = z;     // error
       ^

Weitere Vorüberlegungen

Wie bei jedem Projekt muss in der Anfangsphase ständig re-organisiert werden. Das betrifft zum Beispiel Gliederung in Unterverzeichnisse und einzelne Source-/Header Files:
- Das gesamte Projekt bekommt ein eigenes Verzeichnis hpc. Dies hat zunächst folgende Unterverzeichnisse:
  - hpc/ulmblas wird alle low-level BLAS Funktionen enthalten. Dies sind die BLAS Funktionen, bei denen ohne Klassen sondern direkt mit Zeigern gearbeitet wird. Hier sollen sich C-Programmierier wohlfühlen, selbst dann wenn diese C++ Funktion Templates verwenden.
  - hpc/tests wird verschiedene Tests enthalten.
  Natürlich sollen alle Include-Guards entsprechend angepasst werden. Diese sollen den Pfad relativ zum Projektverzeichnis enthalten und natrülich nur aus Großbuchstaben bestehen. Die Schrägstriche ersetzen wir durch Unterstriche. Für hpc/ulmblas/gecopy.h erhält man so
  #ifndef HPC_ULMBLAS_GECOPY_H #define HPC_ULMBLAS_GECOPY_H 1 // ... #endif // HPC_ULMBLAS_GECOPY_H
- Um Namenskonflikte zu vermeiden, verwenden wir nach wie vor Namensräume. Ein einziger Namensraum wird aber nicht reichen. Wir werden BLAS Funktionen mit low-level Interface (mit Dimensionen und Zeiger) und BLAS Funktionen mit high-level Interface (mit Matrix/Vektor Klassen) implementieren. Wir werden deshlab verschachtelte Namensräume verwenden. Funktionen in hpc/ulmblas sollen im Namensraum hpc::ulmblas liegen. Für hpc/ulmblas/gecopy.h erhält man so
  #ifndef HPC_ULMBLAS_GECOPY_H #define HPC_ULMBLAS_GECOPY_H 1 namespace hpc { namespace ulmblas { // Definition of function gecopy ... } } // namespace ulmblas, hpc #endif // HPC_ULMBLAS_GECOPY_H
- Das gesamte Projekt soll später eine eigenständige Bibliothek werden. Da wir Templates verwenden, wird dies eine sogenannte header-only Library. Es wird also keine statische libhpc.a oder dynamische libhpc.so erzeugt. Damit die Bibliothek von einem Benutzer verwendet werden kann, muss dieser beim übersetzen lediglich den Include-Path beim Übersetzen anpassen:
```
$shell> g++ -std=c++11 -Wall -I/home/numerik/hpc/ws15/uebungen/session12/page01/ gecopy.cc
```
  Dabei ist /home/numerik/hpc/ws15/uebungen/session12/page01/ in diesem Fall das Verzeichnis, in dem die hpc Bibliothek (bestehend aus den Unterverzeichnissen hpc/, hpc/ulmblas, usw.) liegt.
  
  Bei Include-Direktiven soll ab jetzt ausgegangen werden, dass der Include-Pfad das Verzeichnis hpc enthält. Ein Header gecopy.h in hpc/ulmblas kann deshalb stets relativ zum Projektverzeichnis mit
  #include <hpc/ulmblas/gecopy.h>
  eingebunden werden. Andernfalls müsste ein Test gecopy.cc in hpc/tests diesen mit
  #include "../ulmblas/gecopy.h"
  einbinden.

Aufgaben

Legt die Verzeichnis hpc, hpc/tests und hpc/ulmblas an.
Implementiert in hpc/ulmblas eine low-level BLAS Funktion gecopy, die Matrizen kopiert, deren Elemente verschieden sein können. Dabei sollen nur implizite Casts verwendet werden. Warnungen und Fehler sollen also wie oben produziert werden.

Eine minimale Optimierung soll an dieser Stelle ebenfalls eingebaut werden. Zumindest dann, wenn beide Matrizen zeilen- oder spaltenweise im Speicher liegen, sollen die Matrizen dementsprechend auch abgelaufen werden.

Wir werden sehen, dass unser bisheriges Design für eine Matrix-Klasse und zugehörige high-level BLAS Funktionen noch nicht flexibel genug sind. Bevor wir diese verbessern benutzen wir primitive Tests, die ohne Klassen arbeiten. Um euer gecopy zu testen, könnt ihr nachfolgende, primitive Vorlage verwenden. Natürlich sollet ihr dieses Programm auch verstehen und als Gelegenheit nutzen, um zum Beispiel die Verwendung von Lambdas zu wiederholen:

session12/page01/hpc/tests/gecopy.cc

        #include <cstdio>
        #include <hpc/ulmblas/gecopy.h>
        
        template <typename T, typename Func>
        void
        geinit(long m, long n, T *A, long incRowA, long incColA, Func func)
        {
            for (long i=0; i<m; ++i) {
                for (long j=0; j<n; ++j) {
                    A[i*incRowA+j*incColA] = func(m, n, i, j);
                }
            }
        }
        
        void
        print_value(float x)
        {
            printf(" %6.1f", x);
        }
        
        void
        print_value(double x)
        {
            printf(" %6.1lf", x);
        }
        
        template <typename T>
        void
        geprint(long m, long n, const T *A, long incRowA, long incColA)
        {
            for (long i=0; i<m; ++i) {
                for (long j=0; j<n; ++j) {
                    print_value(A[i*incRowA+j*incColA]);
                }
                printf("\n");
            }
        }
        
        int
        main()
        {
            // Element type for A and B
            typedef float   TA;
            typedef double  TB;
        
            // Dimensions of A and B 
            long m       = 7;
            long n       = 8;
        
            // Storage order of A
            long incRowA = 1;
            long incColA = m;
        
            // Storage order of B
            long incRowB = n;
            long incColB = 1;
        
            // Allocate A and B
            TA *A = new TA[m*n];
            TB *B = new TB[m*n];
        
            // Initialize A
            auto fooInit = [](long m, long n, long i, long j) -> long
                           {
                               return m*j+i+1;
                           };
            geinit(m, n, A, incRowA, incColA, fooInit);
        
            // Copy A to B
            hpc::ulmblas::gecopy(m, n, A, incRowA, incColA, B, incRowB, incColB);
        
            // Print B
            geprint(m, n, B, incRowB, incColB);
        }
        

Testet das Verhalten für (mindestens) folgende Fälle:

Typ der Elemente von `A`	Typ der Elemente von `B`
float	double
double	float
float	std::complex<float>
std::complex<float>	float
float	std::complex<double>
std::complex<float>	std::complex<double>