# Basic C++ for numerical computations: vectors

This short note gives a few basic tips on how to efficiently write C++ code with a numerical performance in mind, focusing on vectors.

## Introduction

Quick reminder before we begin: Your code should be as easy to understand as possible, easy to maintain and flexible. While it should be designed from the beginning with performance in mind, especially in terms of how data is stored in memory, flexibility is almost as important, especially in the early stages.

## Vectors

C++ has the Standard Templates Library (STL) that provides some basic containers for storing elements. The vector STL container stores data in the same way the C array does, but with additional functionality.

After defining a variable of type vector, you can use it pretty much as a drop-in replacement of the C array. When optimization is enabled, the performance of vector and the C array is basically the same.

vector has a few advantages:

• it remembers its size
• it allocates and deallocates memory automatically, no need to touch pointers
• fits nicely with other C++ libraries

## Example: $$l_2$$-norm of a vector

Let us first look at the code for computing the $$l_2$$-norm of a vector using the standard C arrays.

#include [itex]
int main() {
int n = 100;
double u[n];
// initialize
for (int i = 0; i < n; ++i) {
u[i] = i;
}
// l2-norm
double accum = 0.;
for (int i = 0; i < n; ++i) {
accum += u[i] * u[i];
}
double norm = sqrt(accum);
}

The same code can be easily rewritten using C++ and vector. Notice that we use u.size(), which returns the size of the vector. No need to remember the size anymore.

#include <vector>
#include <cmath>
using namespace std;
int main() {
int n = 100;
vector<double> u(n);  // create a vector with 100 elements, 0. by
// default
// initialize (note that we use u.size() instead of n)
for (int i = 0; i < u.size(); ++i) {
u[i] = i;
}
// l2-norm
double accum = 0.;
for (int i = 0; i < n; ++i) {
accum += u[i] * u[i];
}
double norm = sqrt(accum);
}

Of course, it makes sense to define a function that computes the $$l_2$$-norm. First the C code:

double l2_norm(double const* u, int n) {
double accum = 0.;
for (int i = 0; i < n; ++i) {
accum += u[i] * u[i];
}
return sqrt(accum);
}

Because the C array does not know its size, we have to pass the size as the second argument of the function.

This becomes easier with vector:

double l2_norm(vector<double> const& u) {
double accum = 0.;
for (int i = 0; i < u.size(); ++i) {
accum += u[i] * u[i];
}
return sqrt(accum);
}

The symbol & makes u a reference to a vector of doubles. This makes the code more efficient. Without &, u would be passed by value, which means that it would be copied into a local variable every time l2_norm is called. However, this is an unnecessary and slow operation. With &, if you call l2_norm in your code, for example:

vector<double> v;
...
double norm = l2_norm(v);

u inside l2_norm is just a different name for the same vector<double> that is stored in v. In particular, if you modify u inside l2_norm, v is modified as well. Hence the const keyword, which tells the caller that v will not be modified, even though it is passed by reference.

If your compiler supports the new C++11 standard (for g++ 4.6 and newer use the -std=c++0x command line option), the code can be simplified somewhat more using the new range-based loop:

double l2_norm(vector<double> const& u) {
double accum = 0.;
for (double x : u) {
accum += x * x;
}
return sqrt(accum);
}

## Example: Discretized Dirichlet energy

Let $$\Omega = (0,1)$$. We want to find the minimizer of the Dirichlet energy

$I(u):= \int_\Omega |u'(x)|^2 \;dx$

over all functions in the set

$\mathcal{K}:=\{v \in {\rm Lip}(\bar \Omega): v(0) = h_0, v(1) = h_1\}.$

Here $${\rm Lip}$$ denotes all Lipschitz continuous functions on $$\bar\Omega$$. Lipschitz continuous functions are differentiable almost everywhere with bounded derivative, and therefore the integral is well-defined and finite.

To discretize this problem, we subdivide the domain into $$N \geq 1$$ equal subintervals $$(x_i, x_{i+1})$$, where $$x_i = i/N$$, $$i = 0, \ldots, N$$. Suppose that $$u$$ is a piece-wise linear function on $$\Omega$$, linear on each $$(x_i, x_{i+1})$$. Function $$u$$ is uniquely prescribed by its values $$v_i:= u(x_i)$$ at the nodes $$x_i$$. We shall assume that the boundary values of $$u$$ at $$x_0 = 0$$ and $$x_N = 1$$ are given values $$h_0, h_1$$. We introduce the discretized Dirichlet energy functional $$J(v_1, \ldots, v_{N-1})$$ as

$J(v_1, \ldots, v_{N-1}):=\int_\Omega |u'(x)|^2 \;dx = \sum_{i=0}^{N-1}\frac{(v_{i+1}-v_i)^2}{\Delta x},$

where $$\Delta x = x_{i+1} - x_i = 1/N$$.

Minimization of $$J$$ over all vectors in $$\mathbb{R}^{N-1}$$ is equivalent to minimization of the original energy $$I$$ over all piece-wise linear functions on $$\Omega$$, linear on each subinterval $$(x_i, x_{i+1})$$, with boundary values $$u(0) = v_0 = h_0$$ and $$u(1) = v_N = h_1$$. It is a discretized version of the original problem in the sense that we restricted the search for the minimizer to a smaller set. Therefore the minimizer of the discretized problem will in general be only an approximation of the real minimizer.

An array or vector are a natural choice for storing the values of $$v_0, \ldots, v_N$$. The number of nodes (vector’s size) is $$N + 1$$.

It is easy to compute the discretized energy $$J$$, given the values $$v_0, \ldots, v_N$$.

double dirichlet_energy(vector<double> const& u) {
int N = u.size() - 1;
double accum = 0;
for (int i = 0; i < N; ++i) {
accum += pow(u[i+1] - u[i], 2.);
}
return N * accum;  // dx = 1 / N
}

The gradient of $$J$$ has $$N-1$$ components:

$\frac{\partial J}{\partial v_i} (v_1, \ldots, v_{N-1}) = -2 \frac{v_{i-1} - 2v_i + v_{i+1}}{\Delta x}.$

However, it is convenient to store it as an array of $$N+1$$ elements and pad it on both sides by zeros.

vector<double> dirichlet_energy_gradient(vector<double> const& u) {
int N = u.size() - 1;
double c = -2. * N;  // c = -2. / dx
vector<double> grad(u.size());  // vector of N+1 zeros
for (int i = 1; i < N; ++i) {
grad[i] = c * (u[i-1] - 2. * u[i] + u[i+1]);
}
}

In a high-performance code one often wants to avoid unnecessary memory allocations and deallocations. Therefore it is better to return the gradient in a preallocated vector, which can be reused multiple times.

void dirichlet_energy_gradient(vector<double>& grad, vector<double> const& u) {
// check that the vectors have the same size
"does not match");
}
int N = u.size() - 1;
double c = -2. * N;  // c = -2. / dx
for (int i = 1; i < N; ++i) {
grad[i] = c * (u[i-1] - 2. * u[i] + u[i+1]);
}
}

We can also provide the function that returns the gradient directly, for convenience.

vector<double> dirichlet_energy_gradient(vector<double> const& u) {
}

Now we write a code to minimize the discretized Dirichlet energy using the steepest descent method (also called gradient descent), for simplicity. The minimization proceeds iteratively, at each iteration $$k$$ we find the minimum of the functional $$J$$ along the line of steepest descent, that is, along the line $$v_{k-1} - t g_{k-1}$$, where $$v_{k-1}$$ is the solution after $$k-1$$ iterations, and $$g_{k-1} = \nabla J(v_{k-1})$$ is the gradient of $$J$$ at $$v_{k-1}$$. The minimum $$t_k$$ along this line must satisfy

$\nabla J(v_{k-1} - t_k g_{k-1}) \cdot g_{k-1} = 0.$

Since $$\nabla J$$ is linear in this case, we can rewrite it as

$0 = \nabla J(v_{k-1}) \cdot g_{k-1} - t_k \nabla J(g_{k-1}) \cdot g_{k-1} = g_{k-1} \cdot g_{k-1} - t_k \nabla J(g_{k-1}) \cdot g_{k-1}.$

This gives a simple formula for $$t_k$$:

$t_k = \frac{\|g_{k-1}\|^2}{\nabla J(g_{k-1}) \cdot g_{k-1}}.$

The following full algorithm is implemented in find_dirichlet_minimizer below:

Set $$v_0$$ to be the initial guess. For $$k = 1, \ldots, k_{\rm max}$$ do:

1. Compute the gradient, $$g_{k-1} = \nabla J(v_{k-1})$$.
2. If $$\|g_{k-1}\| \leq tol$$, exit. $$v_{k-1}$$ is the resulting approximation and $$k-1$$ is the number of iterations performed.
3. Compute the minimizer along the line, $$t_k = \frac{\|g_{k-1}\|^2}{\nabla J(g_{k-1}) \cdot g_{k-1}}$$.
4. Set $$v_k = v_{k-1} - t_k g_{k-1}$$.

The code follows:

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128  #include #include #include #include using namespace std; void check_size(vector const& u, vector const& v) { // check that the vectors have the same size if (u.size() != v.size()) { throw runtime_error("vector sizes do not match"); } } double dirichlet_energy(vector const& u) { int N = u.size() - 1; double accum = 0; for (int i = 0; i < N; ++i) { accum += pow(u[i+1] - u[i], 2.); } return N * accum; // dx = 1 / N } void dirichlet_energy_gradient(vector& grad, vector const& u) { check_size(u, grad); int N = u.size() - 1; double c = -2. * N; // c = -2. / dx for (int i = 1; i < N; ++i) { grad[i] = c * (u[i-1] - 2. * u[i] + u[i+1]); } grad[0] = grad[N] = 0.; // do not forget the padding } double l2_norm(vector const& u) { double accum = 0.; for (int i = 0; i < u.size(); ++i) { accum += u[i] * u[i]; } return sqrt(accum); } double dot(vector const& u, vector const& v) { check_size(u, v); double accum = 0.; for (int i = 0; i < u.size(); ++i) { accum += u[i] * v[i]; } return accum; } void add_to(vector& u, double a, vector const& v) { check_size(u, v); for (int i = 0; i < u.size(); ++i) { u[i] += a * v[i]; } } int find_dirichlet_minimizer(vector& u, double tolerance = 1e-5) { const int max_iter = 100000; vector grad(u.size()); vector ggrad(u.size()); int iter; for (iter = 0; iter < max_iter; ++iter) { dirichlet_energy_gradient(grad, u); double grad_norm = l2_norm(grad); if (grad_norm <= tolerance) { break; } // find the mimimum along the line u - t grad dirichlet_energy_gradient(ggrad, grad); double t = grad_norm * grad_norm / dot(ggrad, grad); // update add_to(u, -t, grad); } return iter; } int main() { // test the minimization int n = 128; vector u(n + 1); for (int i = 0; i <= n; ++i) { double x = static_cast(i) / n; u[i] = exp(x); } // quick test that the energy and gradient are compatible vector grad(u.size()); vector v(u); // copy of u dirichlet_energy_gradient(grad, u); double delta = 1e-3; add_to(v, delta, grad); double Ju = dirichlet_energy(u); double Jv = dirichlet_energy(v); // should be small cout << "Relative Taylor difference = " << ((Ju + delta * dot(grad, grad) - Jv) / Ju) << endl; // find the minimizer int num_iter = find_dirichlet_minimizer(u); // report results cout << "Number of iterations = " << num_iter << endl; double J = dirichlet_energy(u); double J_exact = pow(u[n] - u[0], 2.); cout << "Energy at minimizer = " << J << endl; cout << "Difference from exact energy = " << (J - J_exact) << endl; dirichlet_energy_gradient(grad, u); cout << "Norm of the residual = " << l2_norm(grad) << endl; return 0; }