More on vector and matrix classes

Content

In this session we will work with ready-to-use vector and matrix classes. In the beginning we will focus on the vector class and explain by example how it works.

About the usage of concepts in this session

Note that some of the example will require a C++ compiler (like g++ on theon) that supports the so called concepts feature. You have to compile with the -std=c++17 and -fconcepts flags.

Here you can read more about concepts in C++.

Getting the vector and matrix classes

Create a new directory and copy /home/numerik/pub/hpc/ws18/session14/matvec.hpp to it. For example:

theon$ mkdir hpc_session14
theon$ cd  hpc_session14
theon$ cp /home/numerik/pub/hpc/ws18/session14/matvec.hpp .
theon$ 

You now have:

#ifndef MATVEC_HPP
#define MATVEC_HPP

#include <cassert>
#include <cstddef>


//------------------------------------------------------------------------------
//
// Definition of DenseVector<T>
//

template <typename T, bool isview=false>
struct DenseVector
{
    const std::size_t       length_;
    const std::ptrdiff_t    inc_;
    T * const               data_;

    DenseVector(std::size_t length)
        : length_(length), inc_(1),
          data_(new T[length_])
    {
        static_assert(!isview);
    }

    DenseVector(std::size_t length, T *data, std::ptrdiff_t inc)
        : length_(length), inc_(inc), data_(data)
    {
        static_assert(isview);
    }


    ~DenseVector()
    {
        if constexpr (!isview) {
            delete [] data_;
        }
    }

    // typedef for element type

    using ElementType = T;

    // disabled constructors and operators:

    DenseVector()                                   = delete;
    DenseVector(const DenseVector &rhs)             = delete;
    DenseVector& operator=(const DenseVector &rhs)  = delete;
    DenseVector& operator=(DenseVector &&rhs)       = delete;

    // methods needed for low level BLAS interfaces:

    std::size_t
    length() const
    {
        return length_;
    }

    std::ptrdiff_t
    inc() const
    {
        return inc_;
    }

    ElementType *
    data()
    {
        return data_;
    }

    const T *
    data() const
    {
        return data_;
    }

    // operators for element access

    ElementType &
    operator()(std::size_t i)
    {
        assert(i < length());
        return data()[i*inc()];
    }

    const ElementType &
    operator()(std::size_t i) const
    {
        assert(i < length());
        return data()[i*inc()];
    }

    // methods for creating views

    DenseVector<T, true>
    block(std::size_t i)
    {
        assert(i<=length());
        return DenseVector<T,true>(length()-i, data() + i*inc(), inc());
    }

    DenseVector<const T, true>
    block(std::size_t i) const
    {
        assert(i<=length());
        return DenseVector<const T,true>(length()-i, data() + i*inc(), inc());
    }

    DenseVector<T, true>
    dim(std::size_t n)
    {
        assert(n<=length());
        return DenseVector<T,true>(n, data(), inc());
    }

    DenseVector<const T, true>
    dim(std::size_t n) const
    {
        assert(n<=length());
        return DenseVector<const T,true>(n, data(), inc());
    }
};

template <typename T>
using DenseVectorView = DenseVector<T, true>;

template <typename T>
using DenseVectorConstView = DenseVector<const T, true>;

//------------------------------------------------------------------------------
//
// Definition of GeMatrix<T>
//

enum class Order { ColMajor, RowMajor };

template <typename T, bool isview=false>
struct GeMatrix
{
    const std::size_t       numRows_, numCols_;
    const std::ptrdiff_t    incRow_, incCol_;
    T * const               data_;

    GeMatrix(std::size_t numRows, std::size_t numCols,
             Order order = Order::ColMajor)
        : numRows_(numRows), numCols_(numCols),
          incRow_(order==Order::ColMajor ? 1 : numCols_),
          incCol_(order==Order::ColMajor ? numRows_ : 1),
          data_(new T[numRows_*numCols_])
    {
        static_assert(!isview);
    }

    GeMatrix(std::size_t numRows, std::size_t numCols,
             T *data,
             std::ptrdiff_t incRow, std::ptrdiff_t incCol)
        : numRows_(numRows), numCols_(numCols),
          incRow_(incRow), incCol_(incCol),
          data_(data)
    {
        static_assert(isview);
    }

    ~GeMatrix()
    {
        if constexpr (!isview) {
            delete [] data_;
        }
    }


    // disabled constructors and operators:

    GeMatrix()                                  = delete;
    GeMatrix(const GeMatrix &rhs)               = delete;
    GeMatrix& operator=(const GeMatrix &rhs)    = delete;
    GeMatrix& operator=(GeMatrix &&rhs)         = delete;

    // typedef for element type

    using ElementType = T;

    // methods needed for low level BLAS interfaces:

    std::size_t
    numRows() const
    {
        return numRows_;
    }

    std::size_t
    numCols() const
    {
        return numCols_;
    }

    std::ptrdiff_t
    incRow() const
    {
        return incRow_;
    }

    std::ptrdiff_t
    incCol() const
    {
        return incCol_;
    }

    ElementType *
    data()
    {
        return data_;
    }

    const ElementType *
    data() const
    {
        return data_;
    }

    // operators for element access

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

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

    // methods for creating views

    GeMatrix<T, true>
    block(std::size_t i, std::size_t j)
    {
        assert(i<=numRows());
        assert(j<=numCols());
        return GeMatrix<T,true>(numRows()-i, numCols()-j,
                                data() + i*incRow() + j*incCol(),
                                incRow(), incCol());
    }

    GeMatrix<const T, true>
    block(std::size_t i, std::size_t j) const
    {
        assert(i<=numRows());
        assert(j<=numCols());
        return GeMatrix<const T,true>(numRows()-i, numCols()-j,
                                      data() + i*incRow() + j*incCol(),
                                      incRow(), incCol());
    }

    GeMatrix<T, true>
    dim(std::size_t m, std::size_t n)
    {
        assert(m<=numRows());
        assert(n<=numCols());
        return GeMatrix<T,true>(m, n, data(), incRow(), incCol());
    }

    GeMatrix<const T, true>
    dim(std::size_t m, std::size_t n) const
    {
        assert(m<=numRows());
        assert(n<=numCols());
        return GeMatrix<const T,true>(m, n, data(), incRow(), incCol());
    }

    DenseVector<T, true>
    row(std::size_t i)
    {
        assert(i<numRows());
        return DenseVector<T, true>(numCols(),
                                    data() + i*incRow(),
                                    incCol());
    }

    DenseVector<const T, true>
    row(std::size_t i) const
    {
        assert(i<numRows());
        return DenseVector<const T, true>(numCols(),
                                          data() + i*incRow(),
                                          incCol());
    }

    DenseVector<T, true>
    col(std::size_t j)
    {
        assert(j<numCols());
        return DenseVector<T, true>(numRows(),
                                    data() + j*incCol(),
                                    incRow());
    }

    DenseVector<const T, true>
    col(std::size_t j) const
    {
        assert(j<numCols());
        return DenseVector<const T, true>(numRows(),
                                          data() + j*incCol(),
                                          incRow());
    }
};

template <typename T>
using GeMatrixView = GeMatrix<T, true>;

template <typename T>
using GeMatrixConstView = GeMatrix<const T, true>;


#endif // MATVEC_HPP

Also get the low level BLAS functions from session 13

theon$ cp /home/numerik/pub/hpc/ws18/session14/ulmblas_level*.hpp .
theon$ 

There are actually three kinds of vector classes

In matvec.hpp we define three kinds of vector classes (with a template parameter for the element type):

DenseVector<T>

Vector type that allocates memory when constructed and deallocates memory when destructed.

DenseVectorView<T>

Vector type that references an existing vector. No memory gets allocated/deallocated when such a matrix gets constructed/destructed

DenseVectorConstView<T>

Like DenseVectorView<T> but only read only access to elements is provided.

Vector interface

Assume in the following that x denotes some vector that has one of the types described above. Then you can use the following methods and operators:

x.length()

returns the length

x.data()

returns a pointer to the first element

x.inc()

returns the stride between elements

x(i)

returns a reference to the element with index i

x.block(i)

creates and returns a vector view. The view references all elements beginning at index i.

x.dim(n)

creates and returns a vector view. The view references only the first n elements.

Example (without using vector views) and exercise

#include <matvec.hpp>
#include <printf.hpp>

template <typename Vector>
void
init(Vector &x)
{
    for (std::size_t i=0; i<x.length(); ++i) {
        x(i) = i + 1;
    }
}

template <typename  Vector>
void
print(const Vector &x)
{
    for (std::size_t i=0; i<x.length(); ++i) {
        fmt::printf("%7.2f ", x(i));
    }
    fmt::printf("\n");
}

int
main()
{
    DenseVector<double>     x(10);

    init(x);

    fmt::printf("x=\n");
    print(x);
}

You can compile and run the program by

theon$ g++ -Wall -I. -std=c++17 -fconcepts -o example01 example01.cpp
theon$ ./example01
x=
   1.00    2.00    3.00    4.00    5.00    6.00    7.00    8.00    9.00   10.00 
theon$ 

Exercise part

Make sure this works for you: Compile and run the program as described above

Example (using vector views) and exercise

We modify the code to create a vector view, modify the view and show that this modified the original vector:

#include <matvec.hpp>
#include <printf.hpp>

template <typename Vector>
void
init(Vector &x)
{
    for (std::size_t i=0; i<x.length(); ++i) {
        x(i) = i + 1;
    }
}

template <typename  Vector>
void
print(const Vector &x)
{
    for (std::size_t i=0; i<x.length(); ++i) {
        fmt::printf("%7.2f ", x(i));
    }
    fmt::printf("\n");
}

int
main()
{
    DenseVector<double>     x(10);

    init(x);

    fmt::printf("x=\n");
    print(x);

    DenseVectorView<double>  y = x.block(3).dim(4);

    fmt::printf("y=\n");
    print(y);

    fmt::printf("calling: init(y);\n");
    init(y);

    fmt::printf("y=\n");
    print(y);

    fmt::printf("x=\n");
    print(x);

}

You can compile and run the program by

theon$ g++ -Wall -I. -std=c++17 -fconcepts -o example02 example02.cpp
theon$ ./example02
x=
   1.00    2.00    3.00    4.00    5.00    6.00    7.00    8.00    9.00   10.00 
y=
   4.00    5.00    6.00    7.00 
calling: init(y);
y=
   1.00    2.00    3.00    4.00 
x=
   1.00    2.00    3.00    1.00    2.00    3.00    4.00    8.00    9.00   10.00 
theon$ 

Exercise part

Make sure this works for you: Compile and run the program as described above

Example (using so called rvalue references) and exercise

We change function init such that is also accepts a rvalue reference:

template <typename Vector>
void
init(Vector &&x)
{
    for (std::size_t i=0; i<x.length(); ++i) {
        x(i) = i + 1;
    }
}

now we can pass in main()a temporary view object to function init:

init(x.block(2).dim(4));

Here the complete code:

#include <matvec.hpp>
#include <printf.hpp>

template <typename Vector>
void
init(Vector &&x)
{
    for (std::size_t i=0; i<x.length(); ++i) {
        x(i) = i + 1;
    }
}

template <typename  Vector>
void
print(const Vector &x)
{
    for (std::size_t i=0; i<x.length(); ++i) {
        fmt::printf("%7.2f ", x(i));
    }
    fmt::printf("\n");
}

int
main()
{
    DenseVector<double>     x(10);

    init(x);

    fmt::printf("x=\n");
    print(x);

    init(x.block(2).dim(4));

    fmt::printf("x=\n");
    print(x);
}

You can compile and run the program by

theon$ g++ -Wall -I. -std=c++17 -fconcepts -o example03 example03.cpp
theon$ ./example03
x=
   1.00    2.00    3.00    4.00    5.00    6.00    7.00    8.00    9.00   10.00 
x=
   1.00    2.00    1.00    2.00    3.00    4.00    7.00    8.00    9.00   10.00 
theon$ 

Exercise part

Matrix interface

Analogously to vectors we have matrix type GeMatrix<T>, GeMatrixView<T> and GeMatrixConstView<T>.

Assume in the following that A denotes some matrix that has one of the types described above. Then you can use the following methods and operators:

A.numRows()

returns the number of rows

A.numCols()

returns the number of columns

A.data()

returns a pointer to the first element

A.incRow()

returns the stride between elements in the same column

A.incCol()

returns the stride between elements in the same row.

A(i, j)

returns a reference to the element with index i, j

A.block(i, j)

creates and returns a matrix view. The view references all elements beginning at index i, j.

A.dim(m, n)

creates and returns a matrix view. The view references only the first m rows and ncolumns.

A.col(j)

creates and returns a vector view from the column with index j.

A.row(i)

creates and returns a vector view from the row with index i.

Example (without using matrix views) and exercise

#include <matvec.hpp>
#include <printf.hpp>

template <typename  Matrix>
void
init(Matrix &A)
{
    for (std::size_t i=0; i<A.numRows(); ++i) {
        for (std::size_t j=0; j<A.numCols(); ++j) {
            A(i,j) = i*A.numCols() + j + 1;
        }
    }
}

template <typename  Matrix>
void
print(const Matrix &A)
{
    for (std::size_t i=0; i<A.numRows(); ++i) {
        for (std::size_t j=0; j<A.numCols(); ++j) {
            fmt::printf("%7.2f ", A(i,j));
        }
        fmt::printf("\n");
    }
}

int
main()
{
    GeMatrix<double>     A(4,10);

    init(A);

    fmt::printf("A=\n");
    print(A);
}

You can compile and run the program by

theon$ g++ -Wall -I. -std=c++17 -fconcepts -o example04 example04.cpp
theon$ ./example04
A=
   1.00    2.00    3.00    4.00    5.00    6.00    7.00    8.00    9.00   10.00 
  11.00   12.00   13.00   14.00   15.00   16.00   17.00   18.00   19.00   20.00 
  21.00   22.00   23.00   24.00   25.00   26.00   27.00   28.00   29.00   30.00 
  31.00   32.00   33.00   34.00   35.00   36.00   37.00   38.00   39.00   40.00 
theon$ 

Exercise part

Make sure this works for you: Compile and run the program as described above

Example (using vector views) and exercise

We modify the code to create a vector view, modify the view and show that this modified the original vector:

#include <matvec.hpp>
#include <printf.hpp>

template <typename  Matrix>
void
init(Matrix &A)
{
    for (std::size_t i=0; i<A.numRows(); ++i) {
        for (std::size_t j=0; j<A.numCols(); ++j) {
            A(i,j) = i*A.numCols() + j + 1;
        }
    }
}

template <typename  Matrix>
void
print(const Matrix &A)
{
    for (std::size_t i=0; i<A.numRows(); ++i) {
        for (std::size_t j=0; j<A.numCols(); ++j) {
            fmt::printf("%7.2f ", A(i,j));
        }
        fmt::printf("\n");
    }
}

int
main()
{
    GeMatrix<double>     A(4,10);

    init(A);

    fmt::printf("A=\n");
    print(A);

    GeMatrixView<double>  B = A.block(1, 4).dim(2, 4);

    fmt::printf("B=\n");
    print(B);

    fmt::printf("calling: init(B);\n");
    init(B);

    fmt::printf("B=\n");
    print(B);

    fmt::printf("A=\n");
    print(A);

}

You can compile and run the program by

theon$ g++ -Wall -I. -std=c++17 -fconcepts -o example05 example05.cpp
theon$ ./example05
A=
   1.00    2.00    3.00    4.00    5.00    6.00    7.00    8.00    9.00   10.00 
  11.00   12.00   13.00   14.00   15.00   16.00   17.00   18.00   19.00   20.00 
  21.00   22.00   23.00   24.00   25.00   26.00   27.00   28.00   29.00   30.00 
  31.00   32.00   33.00   34.00   35.00   36.00   37.00   38.00   39.00   40.00 
B=
  15.00   16.00   17.00   18.00 
  25.00   26.00   27.00   28.00 
calling: init(B);
B=
   1.00    2.00    3.00    4.00 
   5.00    6.00    7.00    8.00 
A=
   1.00    2.00    3.00    4.00    5.00    6.00    7.00    8.00    9.00   10.00 
  11.00   12.00   13.00   14.00    1.00    2.00    3.00    4.00   19.00   20.00 
  21.00   22.00   23.00   24.00    5.00    6.00    7.00    8.00   29.00   30.00 
  31.00   32.00   33.00   34.00   35.00   36.00   37.00   38.00   39.00   40.00 
theon$ 

Exercise part

Make sure this works for you: Compile and run the program as described above

Example (using so called rvalue references) and exercise

We again change function init such that is also accepts a rvalue reference:

#include <matvec.hpp>
#include <printf.hpp>

template <typename  Matrix>
void
init(Matrix &&A)
{
    for (std::size_t i=0; i<A.numRows(); ++i) {
        for (std::size_t j=0; j<A.numCols(); ++j) {
            A(i,j) = i*A.numCols() + j + 1;
        }
    }
}

template <typename  Matrix>
void
print(const Matrix &A)
{
    for (std::size_t i=0; i<A.numRows(); ++i) {
        for (std::size_t j=0; j<A.numCols(); ++j) {
            fmt::printf("%7.2f ", A(i,j));
        }
        fmt::printf("\n");
    }
}

int
main()
{
    GeMatrix<double>     A(4,10);

    init(A);

    fmt::printf("A=\n");
    print(A);

    init(A.block(1, 4).dim(2, 4));

    fmt::printf("A=\n");
    print(A);

}

You can compile and run the program by

theon$ g++ -Wall -I. -std=c++17 -fconcepts -o example06 example06.cpp
theon$ ./example06
A=
   1.00    2.00    3.00    4.00    5.00    6.00    7.00    8.00    9.00   10.00 
  11.00   12.00   13.00   14.00   15.00   16.00   17.00   18.00   19.00   20.00 
  21.00   22.00   23.00   24.00   25.00   26.00   27.00   28.00   29.00   30.00 
  31.00   32.00   33.00   34.00   35.00   36.00   37.00   38.00   39.00   40.00 
A=
   1.00    2.00    3.00    4.00    5.00    6.00    7.00    8.00    9.00   10.00 
  11.00   12.00   13.00   14.00    1.00    2.00    3.00    4.00   19.00   20.00 
  21.00   22.00   23.00   24.00    5.00    6.00    7.00    8.00   29.00   30.00 
  31.00   32.00   33.00   34.00   35.00   36.00   37.00   38.00   39.00   40.00 
theon$ 

Exercise part