More on vector and matrix classes

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.

Getting the vector and matrix classes

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

theon$mkdir hpc_session14 theon$ cd  hpc_session14
theon$cp /home/numerik/pub/hpc/ws19/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/ws19/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 Vector type that allocates memory when constructed and deallocates memory when destructed. DenseVectorView Vector type that references an existing vector. No memory gets allocated/deallocated when such a matrix gets constructed/destructed DenseVectorConstView Like DenseVectorView 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 • Make sure this works for you: Compile and run the program as described above • Try to compile it without the modification to function init. 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

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

• Try to compile it without the modification to function init.