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.

#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:

#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.

#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);
}