How to write a CFD library: The vector class

The vector class provides a fundamental building block in our linear system of equations solver library. We have already seen in the previous article how it will be used by our conjugate gradient algorithm, in this week’s article, we look at how to implement the vector class itself. Next week, we’ll close the loop by looking at the matrix implementation as well. For now, though, we concentrate on how to define the vector class interface (its header file), which arithmetic operations will be required, and then, we look at how to implement all of these. By the end, you should have a firm understanding of how to implement a mathematical vector class in C++.

Download Resources

All developed code and resources in this article are available for download. If you are encountering issues running any of the scripts, please refer to the instructions for running scripts downloaded from this website.

• Download: linearAlgebraLib-Part-5-The-vector-class.zip

In this series

• Part 1: Understanding static, dynamic, and header-only C++ libraries
• Part 2: How to write a CFD library: Discretising the model equation
• Part 3: How to write a CFD library: Basic library structure
• Part 4: How to write a CFD library: The Conjugate Gradient class
• Part 5: How to write a CFD library: The vector class
• Part 6: How to write a CFD library: The sparse matrix class
• Part 7: How to write a CFD library: Compiling and testing
• Part 8: How to integrate external CFD libraries in your code
• Part 9: How to handle C++ libraries with a package manager

In this article

How we got here

Thus far in our quest to write a linear algebra solver to solve \mathbf{Ax}=\mathbf{b}, we have looked at how to discretise the 1D heat diffusion equation and we arrived at the discretised equation using the finite volume approximation. Having done that, and using an implicit time integration, we arrived at the need to solve \mathbf{Ax}=\mathbf{b} and we have subsequently looked at writing our own library for that.

Thus far, we have looked at the basic library structure, as well as how to implement the Conjugate Gradient method, along with a brief overview of the theory behind it. We saw during the Conjugate Gradient implementation, that we have a few vector and matrix operations to solve, and in this week’s article, we look at the vector class, and specifically, how to define specific overloaded operators to carry out arithmetic operations.

In terms of our library structure, this is where we are now:

root
├── build
├── linearAlgebraLib
│   ├── linearAlgebraLib.hpp
│   └── src
│       ├── vector.(hpp|cpp)                   // current article
│       ├── sparseMatrixCSR.(hpp|cpp)
│       ├── linearAlgebraSolverBase.(hpp|cpp)  // previous article
│       └── conjugateGradient.(hpp|cpp)        // previous article
├── main.cpp
└── buildAndRun.(sh|bat)

Let’s recap the arithmetic operations required for our vector class, and then follow that up with the header file (interface) and source file (implementation) for our vector class.

Required arithmetic operators for the Vector class

To start us off, let’s repeat the solve() function for the Conjugate Gradient algorithm below. Our goal is to identify any vector operations we need to implement.

Vector ConjugateGradient::solve(unsigned maxIterations, double convergenceThreshold) {
  Vector r0(_rightHandSide.size());
  Vector r1(_rightHandSide.size());
  Vector p0(_rightHandSide.size());
  Vector p1(_rightHandSide.size());
  Vector x(_rightHandSide.size());

  auto &A = _coefficientMatrix;
  auto &b = _rightHandSide;
  
  double alpha = 0.0;
  double beta = 0.0;

  unsigned iteration = 0;

  p1 = b - A * x;
  r1 = b - A * x;

  do {
    r0 = r1;
    p0 = p1;

    alpha = (r0.transpose() * r0) / (p0.transpose() * (A * p0));
    x = x + alpha * p0;
    r1 = r0 - alpha * A * p0;
    beta = (r1.transpose() * r1) / (r0.transpose() * r0);
    p1 = r1 + beta * p0;

    ++iteration;
  } while (iteration < maxIterations && r1.getL2Norm() > convergenceThreshold);

  return x;
}

We can ignore lines 1-15 here, and we are only interested in counting operations that occur. The list below summarises all of these

• Matrix-vector multiplication (returns a vector):
  • Line 16: A * x
  • Line 23: A * p0
  • Line 25: A * p0
• Vector-vector multiplication (returns a scalar):
  • Line 23: r0.transpose() * r0
  • Line 23: p0.transpose() * (A * p0) ((A * p0) returns a vector)
  • Line 26: r1.transpose() * r1
  • Line 26: r0.transpose() * r0
• Scalar-vector multiplication (returns a vector):
  • Line 24: alpha * p0
  • Line 25: alpha * A * p0 ((A * p0) returns a vector)
  • Line 27: beta * p0
• Vector-vector addition (returns a vector):
  • Line 24: x + alpha * p0 ((alpha * p0) returns a vector)
  • Line 27: r1 + beta * p0 ((beta * p0) returns a vector)
• Vector-vector subtraction (returns a vector):
  • Line 16/17: b - A * x ((A * x) returns a vector)
  • Line 25: r0 - alpha * A * p0 ((alpha * A * p0) returns a vector)
• Vector-vector equality
  • Line 20: r0 = r1
  • Line 21: p0 = p1

That’s quite a list, but it essentially boils down to 5 operations that we have to implement: Matrix-vector multiplication, vector-vector multiplication, scalar-vector multiplication, vector-vector addition as well as subtraction. Why do I leave off the equality operation here? Because C++ will provide that implementation for us automatically.

Also note the order of the operation, having scalar-vector and vector-scalar multiplications, for example, will provide the same result in a mathematical sense but would need to be implemented with two separate functions (in our case overloaded operators). So it is important to either keep the scalar always to the left / right of the vector and then provide the corresponding implementation, or two provide both. Otherwise, you get a compilation error saying the function you thought you just implemented does not exist.

With these operators defined, we are in a good position to perform operations with our vector, so let’s look at the class interface next.

Vector class interface

Apart from the arithmetic operators listed above, we also want to be able to calculate the L2-norm for our vector, which will essentially tell us the magnitude of a vector. We’ll use that to check convergence, as seen in the code above on line 30. We want to be able to transpose the vector, which will not change its underlying data but just denote it as either a row or column vector. We’ll also need a function that tells us the size of the vector, otherwise, we would be unable to loop over it. Finally, it would be good to have a function to print the vector to the console, so we’ll provide a function for that as well.

The vector class, thus, is given below with all the required functions that we have just looked at.

#pragma once

#include <iostream>
#include <vector>
#include <cassert>
#include <cmath>

namespace linearAlgebraLib {

class Vector {
public:
  using VectorType = std::vector<double>;

public:
  Vector(const unsigned &size);

  unsigned size() const;

  Vector transpose();
  double getL2Norm() const;

  const double &operator[](unsigned index) const;
  double &operator[](unsigned index);

  Vector operator+(const Vector &other);
  Vector operator-(const Vector &other);

  Vector &operator*(const double &scaleFactor);
  friend Vector operator*(const double &scaleFactor, Vector vector);  
  double operator*(const Vector &other);
  
  friend std::ostream &operator<<(std::ostream &out, const Vector &vector);

private:
  VectorType _vector;
  bool _isRowVector = false;
};

} // namespace linearAlgebraLib

On line 12, we are defining the underlying vector type as a std::vector that contains doubles. This definition is then used on line 35 to hold the values of our vector. On line 36, we’ll also see that the default construction of a vector assumes that it is a column vector. We have a constructor on line 15 that will resize the underlying vector to the correct size. The size() function itself will return the size of the vector, the transpose() function will simply change the _isRowVector variable on line 36, and the getL2Norm() function calculates the L2-norm.

Afterwards, we define all our overloaded operators on lines 22-32. So we see vector-vector addition and subtraction defined on lines 25 and 26, respectively, vector-scalar multiplication on line 28, scalar-vector multiplication on line 29, and vector-vector multiplication (dot-product) on line 30. We define a friend function to the operator<< on line 32, which will allow us to print our vector to the console. We looked at this function already in somewhat more detail in our operating overloading article, which you may want to check for further details.

Here is the short version of friend functions: We see on line 29 that we also have to befriend the scalar-vector multiplication. Remember that this function should allow us to write code like vector = scalar * vector;. However, we could also write this statement as vector = scalar.operator*(vector);. The problem is that the scalar (here of type double) does not have an operator*() function defined. Therefore, by using the friend keyword, we are saying that the Vector and, well, double class (it is not actually a class) are friends and can define functions for each other.

The same holds true for the matrix-vector multiplication; we haven’t defined it in this class but writing a statement like vector = matrix * vector; could also be written as vector = matrix.operator*(vector);. So it makes sense to define this multiplication in the matrix class (and this is where we’ll implement it), but we could have also defined it in the vector class as a friend function and then provided the implementation here. At least technically speaking this would be possible, if you do it in this particular example, you get a nasty compilation error, which will not be obvious.

Not to go too deep into the issue (but I figured I mention it here in case someone is trying to do exactly that), but we would get a circular dependency. The matrix class depends on the vector class (for the matrix-vector multiplication), but we’ll also introduce the matrix class into the vector class, therefore, the vector class depends on the matrix class as well. We have to include the header of the matrix in the vector class, and the header of the vector in the matrix class. To resolve this issue, we have to start using pointers where possible and use forward declarations of classes instead of includes.

A final note on the operator[] that is defined on lines 22-23; we specify two versions here to allow to get, but also, to set values. The version on line 22 returns a const reference (can’t be modified), but the version on line 23 returns just a reference, which can be modified, and hence be used to set values. If you want to understand this in more depth, there is a good explanation about this on learncpp.

Having looked at all these definitions, let’s look at the corresponding implementations.

This website exists to create a community of like-minded CFD enthusiasts and I’d love to start a discussion with you. If you would like to be part of it, sign up using the link below and you will receive my OpenFOAM quick reference guide, as well as my guide on Tools every CFD developer needs for free.

Join now

Vector class implementation

After looking at the interface (header file), let’s have a look at the function implementations. We’ll start by including the header file and wrapping all implementations in a namespace, as shown below:

#include "linearAlgebraLib/src/vector.hpp"

namespace linearAlgebraLib {

// function implementations go here ...

} // namespace linearAlgebraLib

All subsequent implementations would go between lines 3-7.

Required class utility functions

Vector::Vector(const unsigned &size) {
  _vector.resize(size, 0.0);
}

The first function we need to implement is the constructor, which simply takes the size of the vector as an argument and then resizes the underlying _vector variable, which is of type std::vector<double> as defined before.

unsigned Vector::size() const {
  return _vector.size();
}

The next function is the size() function which simply returns the size of the vector.

To access variables within the _vector variable, we overload the operator[] twice, which allows us to set and get values, as mentioned above. This may seem a bit strange if you see it for the first time and it will take some time to digest this, I can only recommend the link given above to get a good understanding of why this is possible.

double &Vector::operator[](unsigned index) {
  return _vector[index];
}

const double &Vector::operator[](unsigned index) const {
  return _vector[index];
}

Both functions actually do exactly the same thing, but the compiler will figure out which version to call and then either allows to only read the content or also to set the content.

We also want to provide an overloaded operator for the output stream operator (operator<<). We don’t use it in the code anywhere, it is just a useful way to quickly see the content of the vector during debugging, at which point this function was used quite heavily. We discussed how this operator overloading works when we talked about operator overloading by itself.

std::ostream &operator<<(std::ostream &out, const Vector &vector) {
  out << "( ";
  for (unsigned index = 0; index < vector.size(); ++index)
    out << vector[index] << " ";
  out << ")";
  return out;
}

Mathematical utility functions

Vector Vector::transpose() {
  Vector transposedVector(_vector.size());
  transposedVector._vector = _vector;
  transposedVector._isRowVector = _isRowVector == true ? false : true;
  return transposedVector;
}

The transpose function returns a new vector where the _isRowVector variable is either set to true or false, depending on what the value is for the current vector. Now, you may be saying that this is wasteful, let’s say we have a domain size of 10 million cells, then our vector (for example, of the temperature here) will have 10 million elements as well and we are attempting to create a new vector with 10 million cells. So, instead, you may be tempted to give the following implementation:

Vector &Vector::transpose() {
  _isRowVector = _isRowVector == true ? false : true;
  return *this;
}

In this version, we set the _isRowVector variable of the vector that is calling the transpose() function directly, and then return a reference to itself with the *this return value. The reason this won’t work (in our case), is that we are using the following line in the conjugate gradient algorithm:

alpha = (r0.transpose() * r0) / (p0.transpose() * (A * p0));

If we were to return a reference to r0, which is making the call to transpose(), then we are not solving r^Tr, but actually r^Tr^T, since we have changed r0 directly. This is not a mathematical operation for vectors and if you were to run this, you’ll get an assertion error that you need to have a column and row vector for this type of multiplication.

The remedy is that if you wanted to improve performance here, you first have to realise that r^Tr can be replaced by the dot product r\cdot r, and then use that instead (for which we’ll see the corresponding implementation in just a second).

The next function we’ll need to add is the one to calculate the L2-norm, which for vectors is defined as

\lvert\lvert\mathbf{x}\rvert\rvert_2=\sqrt{\sum_{i=0}^{size(\mathbf{x})}x_i^2}

In code, this looks like the code below:

double Vector::getL2Norm() const {
  double L2Norm = 0.0;
  for (const auto &entry : _vector)
    L2Norm += std::pow(entry, 2);
  L2Norm = std::sqrt(L2Norm);
  return L2Norm;
}

Arithmetic overloaded operators

In the next few sections, we’ll look at the remaining overloaded operators for the arithmetic operations with our vector.

Vector-Vector addition

Vector-vector addition simply adds the individual components of the vector together, resulting in a new vector given as

\begin{pmatrix}a_1 \\ a_2 \\ a_3 \\\vdots \\ a_n\end{pmatrix} + \begin{pmatrix}b_1 \\ b_2 \\ b_3 \\\vdots \\ b_n\end{pmatrix} = \begin{pmatrix}a_1 + b_1 \\ a_2 + b_2 \\ a_3 + b_3 \\\vdots \\ a_n + b_n\end{pmatrix}

In code, we can translate this into

Vector Vector::operator+(const Vector &other) {
  assert(_vector.size() == other._vector.size() && "vectors must have the same dimension");
  Vector resultVector(_vector.size());
  resultVector._vector.resize(_vector.size());
  for (unsigned i = 0; i < _vector.size(); ++i)
    resultVector._vector[i] = _vector[i] + other._vector[i];
  return resultVector;
}

It is worth pointing out here, that we use an assert() statement at the beginning. This is coming from the fail-fast programming philosophy. The idea here is that we check if certain conditions, that we know must be met, are correct. If they are not, then we will stop execution, before we get the wrong output and then try to figure out what went wrong. In this way, we know exactly at what point the execution went wrong and we can start our debugging quest from here. The advantage of using assert() statements is that we can turn them off globally, which we typically want for production runs after we made sure the code is correct.

Vector-Vector subtraction

Vector-vector subtraction works pretty much the same as the addition case, we just exchange the plus operator for the minus operator

\begin{pmatrix}a_1 \\ a_2 \\ a_3 \\\vdots \\ a_n\end{pmatrix} - \begin{pmatrix}b_1 \\ b_2 \\ b_3 \\\vdots \\ b_n\end{pmatrix} = \begin{pmatrix}a_1 - b_1 \\ a_2 - b_2 \\ a_3 - b_3 \\\vdots \\ a_n - b_n\end{pmatrix}

This results again in very similar code

Vector Vector::operator-(const Vector &other) {
  assert(_vector.size() == other._vector.size() && "vectors must have the same dimension");
  Vector resultVector(_vector.size());
  resultVector._vector.resize(_vector.size());
  for (unsigned i = 0; i < _vector.size(); ++i)
    resultVector._vector[i] = _vector[i] - other._vector[i];
  return resultVector;
}

Vector-Vector multiplication

This case is a bit more interesting, as we have a few different cases to discuss. First, let’s start with the dot product, which we have already briefly touched upon above. It returns a scalar and is defined as

\mathbf{a}\cdot \mathbf{b} = a_1 b_1 +a_2 b_2 + a_3 b_3 + \dots + a_n b_n

In code, we would translate the above equation in the following form

double Vector::operator*(const Vector &other) {
  assert(_isRowVector && !other._isRowVector && "first vector must be a row vector, second must be a column vector");
  double sum = 0.0;
  for (unsigned i = 0; i < _vector.size(); ++i)
    sum += _vector[i] * other._vector[i];
  return sum;
}

As discussed before, we check again that the vectors have the same size in line with the fail-fast philosophy and after that condition is true, we continue with the dot product calculation.

Scalar-Vector multiplication

This is the aforementioned more interesting case. For scalar-vector multiplication, we have to distinguish between scalar-vector and vector-scalar multiplication, both of them result in different code, which we discussed above as well when we looked at the class interface (header file). Let’s start with the simple vector-scalar multiplication, which is given as

\mathbf{a}\cdot s=\begin{pmatrix}a_1 * s_1 \\ a_2 * s_2 \\ a_3 * s_3 \\\vdots \\ a_n * s_n\end{pmatrix}

This is put into code as

Vector &Vector::operator*(const double &scaleFactor) {
  for (unsigned index = 0; index < _vector.size(); ++index)
    _vector[index] *= scaleFactor;
  return *this;
}

The reverse case, i.e. scalar-vector multiplication is slightly more complicated as we saw above, as we now need to make this multiplication implementation a friend function. The actual implementation looks pretty much the same, we just need to change the arguments to the function. The code becomes

Vector operator*(const double &scaleFactor, Vector vector) {
  for (unsigned index = 0; index < vector.size(); ++index)
    vector[index] *= scaleFactor;
  return vector;
}

In the second code, we get a vector passed to the function (scalar-vector multiplication, remember, in this case, we could write scalar * vector as scalar.operator*(vector), hence we get the vector as an argument). In the first code, on the other hand, we don’t, and we use the class vector data (stored in the _vector variable). This is because vector * scalar can also be written as vector.operator*(scalar), so we are calling the operation from the vector object, which can access its own state (data).

That’s it, if you followed everything up until now, you have a good understanding of how we can implement our vector class. Next, we’ll look at the matrix class and see how we can provide an implementation for it.

Summary

There was quite a bit of text here but I hope you get the idea for how to write a vector class. This is half of the puzzle done towards creating our linear algebra solver library. Once we have defined the matrix class, which we do in the next article, you’ll have all the pieces required to use it for a real CFD application (or stick with the 1D heat equation we have defined for the moment).