======================================
LU-Zerlegung mit Matrix/Vektor-Klassen
======================================

Bisher wurde die LU-Zerlegung mit den Low-Level BLAS Funktionen im Namensraum
`hpc::ulmblas` implementiert.  Dies soll nun abstrakter möglich sein, so dass
bei der Implementierung technische Details wie Zeiger und Inkremente versteckt
werden.  Wie dies möglich ist, soll am Beispiel der ungeblockten LU-Zerlegung
gezeigt werden.  Dazu werden folgende Schritte notwendig sein:

- Implementierung einer Vektor Klasse `DenseVector` mit zugehörigen View-Klassen
  `DenseVectorView` und `DenseVectorConstView`.
- Anpassen der Matrix-Klassen, so dass Vektor-Views erzeugt werden, die eine
  Zeile oder Spalte referenzieren.
- High-Level BLAS Wrapper für die BLAS Funktionen `scal`, `swap`, `iamax` und
  `ger`.


Die DenseVector-Klassen
=======================

Benutzt nachfolgende Vorlage für die Vektor-Klassen.  Mit dem Funktions-Operator
`x(i, num)` soll eine View erzeugt werden, die die Elemente mit Indices
$i, \dots, i+num-1$ referenziert.

==== CODE (type=cc) ============================================================
#ifndef HPC_MATVEC_DENSEVECTOR_H
#define HPC_MATVEC_DENSEVECTOR_H 1

#include <cassert>
#include <cstdlib>

namespace hpc { namespace matvec {

template <typename T, typename I>
    struct DenseVectorConstView;

template <typename T, typename I>
    struct DenseVectorView;

template <typename T, typename I=std::size_t>
struct DenseVector
{
    typedef T                               ElementType;
    typedef I                               Index;

    typedef DenseVector<T,Index>            NoView;
    typedef DenseVectorConstView<T,Index>   ConstView;
    typedef DenseVectorView<T,Index>        View;

    DenseVector(Index length);

    ~DenseVector();

    const ElementType &
    operator()(Index i) const;

    ElementType &
    operator()(Index i);

    ConstView
    operator()(Index i, Index num, Index stride=1) const;

    View
    operator()(Index i, Index num, Index stride=1);

    const Index     length, inc;
    ElementType*    data;

};

template <typename T, typename I=std::size_t>
struct DenseVectorView
{
    typedef T                           ElementType;
    typedef I                           Index;

    typedef DenseVector<T,Index>            NoView;
    typedef DenseVectorConstView<T,Index>   ConstView;
    typedef DenseVectorView<T,Index>        View;

    DenseVectorView(Index length, ElementType *data, Index inc);

    const ElementType &
    operator()(Index i) const;

    ElementType &
    operator()(Index i);

    ConstView
    operator()(Index i, Index num, Index stride=1) const;

    View
    operator()(Index i, Index num, Index stride=1);

    const Index     length, inc;
    ElementType*    data;
};

template <typename T, typename I=std::size_t>
struct DenseVectorConstView
{
    typedef T                           ElementType;
    typedef I                           Index;

    typedef DenseVector<T,Index>            NoView;
    typedef DenseVectorConstView<T,Index>   ConstView;
    typedef DenseVectorView<T,Index>        View;

    DenseVectorConstView(Index length, const ElementType *data, Index inc);

    const ElementType &
    operator()(Index i) const;

    ConstView
    operator()(Index i, Index num, Index stride=1) const;

    const Index         length, inc;
    const ElementType*  data;
};


} } // namespace matvec, hpc

#endif // HPC_MATVEC_DENSEVECTOR_H
================================================================================

Traits
======
Analog zu den Matrix-Klassen benötigen wir eine Traits-Klasse, die zur
Compile-Zeit feststellt, ob es sich bei einem Template-Parameter um eine
Vektor-Klasse handelt.  Diese könnt ihr direkt übernehmen:


==== CODE (type=cc) ============================================================
#ifndef HPC_MATVEC_ISDENSEVECTOR_H
#define HPC_MATVEC_ISDENSEVECTOR_H 1

#include <cassert>
#include <type_traits>
#include <hpc/aux/iscomplex.h>
#include <hpc/matvec/densevector.h>

namespace hpc { namespace matvec {

template <typename Any>
struct IsDenseVector_
{
    static constexpr bool value = false;
};

template <typename T, typename I>
struct IsDenseVector_<DenseVector<T,I> >
{
    static constexpr bool value = true;
};

template <typename T, typename I>
struct IsDenseVector_<DenseVectorView<T,I> >
{
    static constexpr bool value = true;
};

template <typename T, typename I>
struct IsDenseVector_<DenseVectorConstView<T,I> >
{
    static constexpr bool value = true;
};

template <typename Any_>
struct IsDenseVector
{
    typedef typename std::remove_reference<Any_>::type  Any;
    static constexpr bool value = IsDenseVector_<Any>::value;
};

template <typename Vector>
struct IsRealDenseVector
{
    typedef typename Vector::ElementType    T;

    static constexpr bool value = IsDenseVector<Vector>::value
                               && !aux::IsComplex<T>::value;
};

template <typename Vector>
struct IsComplexDenseVector
{
    typedef typename Vector::ElementType    T;

    static constexpr bool value = IsDenseVector<Vector>::value
                               && aux::IsComplex<T>::value;
};

} } // namespace matvec, hpc

#endif // HPC_MATVEC_ISDENSEVECTOR_H
================================================================================


Aufgabe
=======
Folgendes Test-Programm soll benutzt werden, um die Implementierung zu testen.
Dabei sollt Ihr feststellen, welchen konkreten Typ die mit `auto` deklarierten
Variablen tatsächlich besitzen.  Werden alle Operatoren zum Erzeugen einer
Vektor-View aufgerufen?  Falls nicht, ergänzt das Programm entsprechend.


==== CODE (type=cc) ============================================================
#include <cassert>
#include <random>
#include <type_traits>
#include <hpc/matvec/densevector.h>
#include <hpc/matvec/print.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);
        }
    }
}

template <typename VX>
typename std::enable_if<hpc::matvec::IsDenseVector<VX>::value,
                        void>::type
randomInit(VX &x)
{
    typedef typename VX::Index  Index;

    randomInit(x.length, Index(1), x.data, x.inc, Index(1));
}

//------------------------------------------------------------------------------

template <typename VX>
typename std::enable_if<hpc::matvec::IsDenseVector<VX>::value,
                        void>::type
foo(const VX &x)
{
    std::printf("Entering foo\n");
    auto y = x(3, 7);
    auto z = y(0, 4, 2);

    print(x, "x");
    print(y, "y");
    print(z, "z");
    std::printf("Leaving foo\n");
}

int
main()
{
    using namespace hpc::matvec;

    typedef double       T;
    typedef std::size_t  Index;

    DenseVector<T, Index> x(10);

    auto y = x(3, 7);
    auto z = y(0, 4, 2);

    randomInit(x);

    print(x, "x");
    print(y, "y");
    print(z, "z");

    foo(x);
}
================================================================================

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