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