Content

Vektor-Views und Matrizen

Ist A eine Matrix vom Typ GeMatrix, GeMatrixView oder GeMatrixConstView, so soll mit A.row(i) bzw. A.col(j) jeweils eine Vektor-View erzeugt werden, die die \(i\)-te Zeile bzw. \(j\)-te Spalte referenziert.

Vorlage

Als Vorlage könnt Ihr Euch an folgendem Code-Segment orientieren. Es zeigt allerdings nur die Klasse GeMatrix. Analog sollen die anderen Varianten GeMatrixView und GeMatrixConstView erweitert werden.

/*
   ...
*/

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

    typedef GeMatrix<T,Index>              NoView;
    typedef GeMatrixConstView<T,Index>     ConstView;
    typedef GeMatrixView<T,Index>          View;

    typedef DenseVectorConstView<T,Index>  ConstVectorView;
    typedef DenseVectorView<T,Index>       VectorView;


    GeMatrix(Index numRows, Index numCols,
             StorageOrder order=StorageOrder::ColMajor)
        : numRows(numRows), numCols(numCols),
          incRow(order==StorageOrder::ColMajor ? 1: numCols),
          incCol(order==StorageOrder::RowMajor ? 1: numRows),
          data(new T[numRows*numCols])
    {
    }

    ~GeMatrix()
    {
        delete[] data;
    }

    const ElementType &
    operator()(Index i, Index j) const
    {
        assert(i<numRows && j<numCols);
        return data[i*incRow + j*incCol];
    }

    ElementType &
    operator()(Index i, Index j)
    {
        assert(i<numRows && j<numCols);
        return data[i*incRow + j*incCol];
    }

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

    View
    operator()(Index i, Index j, Index m, Index n);

    ConstVectorView
    row(Index i) const;

    VectorView
    row(Index i);

    ConstVectorView
    col(Index j) const;

    VectorView
    col(Index j);

    const Index     numRows, numCols, incRow, incCol;
    ElementType*    data;
};

/*
   ...
*/

Traits

Die IsGematrix müssen ergänzt werden, wenn auch Rvalue-Referenzen verwendet werden sollen. Diese könnt ihr direkt übernehmen:

#ifndef HPC_MATVEC_ISGEMATRIX_H
#define HPC_MATVEC_ISGEMATRIX_H 1

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

namespace hpc { namespace matvec {

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

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

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

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

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


template <typename Matrix>
struct IsRealGeMatrix
{
    typedef typename Matrix::ElementType    T;

    static constexpr bool value = IsGeMatrix<Matrix>::value
                               && !aux::IsComplex<T>::value;
};

template <typename Matrix>
struct IsComplexGeMatrix
{
    typedef typename Matrix::ElementType    T;

    static constexpr bool value = IsGeMatrix<Matrix>::value
                               && aux::IsComplex<T>::value;
};

} } // namespace matvec, hpc

#endif // HPC_MATVEC_ISGEMATRIX_H

Aufgabe

Folgendes Test-Programm soll benutzt werden, um die Implementierung zu testen. Dabei sollt Ihr wieder feststellen, welchen konkreten Typ die mit auto deklarierten Variablen tatsächlich besitzen. Wie zuvor soll das Programm auch so abgeändert werden, dass alle Fälle getestet werden.

#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 MA>
typename std::enable_if<hpc::matvec::IsGeMatrix<MA>::value,
                        void>::type
randomInit(MA &A)
{
    randomInit(A.numRows, A.numCols, A.data, A.incRow, A.incCol);
}


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

template <typename MA>
typename std::enable_if<hpc::matvec::IsGeMatrix<MA>::value,
                        void>::type
foo(const MA &A)
{
    std::printf("Entering foo\n");

    auto m = A.numRows;
    auto n = A.numCols;

    auto B = A(0, 0, m, n);

    auto x = B.row(0);
    auto y = B.col(0);

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

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

    typedef double       T;
    typedef std::size_t  Index;

    GeMatrix<T, Index> A(8,10);

    auto x = A.row(3);
    auto y = A.col(4);

    randomInit(A);

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

    auto B = A(1, 2, 6, 5);

    auto v = B.row(3);
    auto w = B.col(4);

    print(B, "B");
    print(v, "v");
    print(w, "w");

    foo(B);
}