7. Arrays, Lambdas, and Parallel Loops on the CPU
Consider allocating an array in Fortran and computing the value at each grid point in a loop. In MISO, this processing is split into the following three parts.
Prepare the array’s memory.
Array1DthroughArray4Dhandle this.Take out a lightweight “window onto the array” to pass to the computation. This is a View (
Array3DViewand so on).Write the processing to be done at each grid point as a lambda expression, and pass it to
for_eachtogether with the range of indices.
The reason for splitting into three is so that the same computation can run on both the CPU and the GPU. The place to compute is selected with the tag types backend::Host (CPU) and backend::CUDA (GPU) that you learned in Chapter 6. This chapter uses only backend::Host, and goes as far as computing on arrays on the CPU and checking the results. How to run the same computation on the GPU is explained in Chapter 9.
The examples in this chapter use only MISO’s headers. If you add miso/include to the search path, you can build directly with g++ as follows (replace /path/to/MISO with the place where you put MISO).
g++ -std=c++17 -I/path/to/MISO/miso/include layout.cpp -o layout
./layout
The headers read in this chapter (array1d.hpp, array3d.hpp, execution.hpp, and utility.hpp) do not require MPI or yaml-cpp. For the CMake setup to use MISO, the end-of-chapter Practice uses the setup from Chapter 1.
Contents
1. Arrays Array1D through Array4D
1.1 Allocation and Indices
In Chapter 4, you learned that Array3D is a class that owns an array. MISO has four array classes, one for each number of dimensions.
Class |
Header |
Element access |
|---|---|---|
|
|
|
|
|
|
|
|
|
|
|
|
T is the element type, and Backend is where the array is placed. If you omit Backend, it becomes backend::Host (the declaration at the top of array3d.hpp, template <typename T, typename Backend = backend::Host> class Array3D;). Only the one-dimensional Array1D has no parentheses, and specifies elements with square brackets, a[i].
miso::Array1D<double> x(100); // 100 doubles: x[0] to x[99]
miso::Array3D<double> ro(16, 32, 8); // ro(0,0,0) to ro(15,31,7)
Indices start from 0 in every direction, and the last one is “size − 1.” This corresponds to Fortran’s real(8) :: ro(0:15, 0:31, 0:7).
Allocated elements are not initialized. The constructor of the Host version only allocates with data_ = new T[array_size]; and does not put in 0 (new from Chapter 3). As with Fortran’s allocate, always assign a value before reading.
1.2 Layout in Memory
A three-dimensional array is also laid out as a single line in memory. Where a(i, j, k) falls in that line is written in operator() in array3d.hpp.
// excerpt from Array3D<T, backend::Host> in miso/include/miso/array3d.hpp
T &operator()(int i0, int i1, int i2) const noexcept {
...
return data_[(i0 * shape_[1] + i1) * shape_[2] + i2];
}
shape_ holds the sizes {nx0, nx1, nx2} in each direction. For an array of size (nx, ny, nz), a(i, j, k) is at position (i * ny + j) * nz + k from the start. When the last index k increases by 1, you move to the next element; when j increases by 1, you advance nz; and when i increases by 1, you advance ny * nz. In other words, the last index k is contiguous in memory.
Fortran is the opposite: the first index i is contiguous (column-major). A layout in which the last index is contiguous, as in C++ multidimensional arrays and MISO’s arrays, is called row-major.
We check with the following program. a[n] is a way of directly reading the n-th element from the start. &a(0, 0, 1) - &a(0, 0, 0) expresses the difference between the addresses of two elements as a number of elements (the pointer arithmetic from Chapter 3).
#include <iostream>
#include <miso/array3d.hpp>
int main() {
miso::Array3D<double> a(2, 3, 4); // omitting the 2nd type gives backend::Host
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 3; ++j) {
for (int k = 0; k < 4; ++k) {
a(i, j, k) = 100 * i + 10 * j + k;
}
}
}
// read the first 6 elements in memory by linear index a[n]
for (int n = 0; n < 6; ++n) {
std::cout << a[n] << " ";
}
std::cout << "\n";
// how many elements in memory one index increment advances
std::cout << "k: " << &a(0, 0, 1) - &a(0, 0, 0) << "\n";
std::cout << "j: " << &a(0, 1, 0) - &a(0, 0, 0) << "\n";
std::cout << "i: " << &a(1, 0, 0) - &a(0, 0, 0) << "\n";
}
0 1 2 3 10 11
k: 1
j: 4
i: 12
From the start, the order is a(0,0,0), a(0,0,1), a(0,0,2), a(0,0,3), a(0,1,0), a(0,1,1).
The CPU reads nearby values in memory together. For that reason, loops that process contiguous elements in order are faster. In Fortran, the innermost loop was over i, but with MISO arrays you make the innermost loop over k. MISO’s Integrator::cfl (miso/include/miso/mhd_integrator_host.hpp) also nests the loops in the order i, j, k from the outside.
For Array2D, a(i, j) is at position i * ny + j, and for Array4D, the last index l of a(i, j, k, l) is contiguous.
Exercise 7-1: For miso::Array3D<double> a(4, 5, 6);, find the position of a(1, 2, 3) from the start of memory. Also give the elements immediately adjacent to a(1, 2, 3) (the one before and the one after).
The size is (nx, ny, nz) = (4, 5, 6), so the position is (i * ny + j) * nz + k = (1 * 5 + 2) * 6 + 3 = 45 (counting from 0). The last index k is the contiguous one, so the immediate neighbors are a(1, 2, 2) (position 44) and a(1, 2, 4) (position 46). a(0, 2, 3) and a(2, 2, 3) are ny * nz = 30 elements away.
1.3 Taking Out a View
As seen in Chapter 4, Array3D prohibits copying (Array3D(const Array3D &) = delete;). This is because if it were copied, two objects would try to release the same memory. On the other hand, the parallel loops you will learn in the following sections require passing the array by value. So we take out a View, a lightweight non-owning window, from the array and pass that.
miso::Array3D<double> ro(16, 32, 8);
auto rov = ro.view(); // Array3DView<double>: read/write
auto roc = ro.const_view(); // Array3DView<const double>: read-only
Array3DView holds only the start pointer data_ and the size shape_[3]. Copying it does not duplicate the elements, and the copied View points to the same elements. The index calculation is the same as for Array3D, (i0 * shape_[1] + i1) * shape_[2] + i2. A View does not release memory, so you must keep the original Array3D alive while you are using the View.
What you take out of an Array1D is an Array1DView. There are also Array2DView and Array4DView. mhd::FieldsView, which bundles MISO’s physical variables (nine, such as density, velocity, and magnetic field), is a struct with nine Array3DViews (miso/include/miso/mhd_fields.hpp).
operator() and operator[] of a View have const, but the return value is T &. As with double* const in Chapter 3, even though the View itself (the pointer and size) cannot be changed, the elements it points to can be rewritten. We use this property in the lambda expressions of the next section.
2. Lambda Expressions
To for_each, you pass “the processing to be done at each grid point.” To pass processing as if it were a value, C++ uses lambda expressions. All of MISO’s parallel computation is written with lambda expressions, so this section explains them step by step.
2.1 Unnamed Functions
In Fortran, a short function used only inside a subroutine can be written as an internal procedure after contains.
subroutine scale_all(n, a, factor)
integer, intent(in) :: n
real(8), intent(inout) :: a(n)
real(8), intent(in) :: factor
integer :: i
do i = 1, n
a(i) = scaled(a(i))
enddo
contains
real(8) function scaled(x)
real(8), intent(in) :: x
scaled = factor * x ! uses the parent procedure's factor directly
end function scaled
end subroutine scale_all
In C++, you can write a function on the spot at an expression position inside a function. This is a lambda expression.
auto square = [](double x) { return x * x; }; // takes x and returns x*x
std::cout << square(3.0) << "\n"; // 9
Part |
Meaning |
|---|---|
|
The capture specification (next section). Here it is empty because no outside variables are used |
|
The argument list. The same as an ordinary function |
|
The body of the processing. The same as an ordinary function |
The return type is not written. The compiler decides double from the return expression.
The internal procedure scaled had a name. The function of a lambda expression itself has no name. In the example above, the “callable object” made by the lambda expression is put in the variable square and called with square(3.0). The type of this object is an unnamed type given by the compiler, so auto is used to declare the variable. You can also write it on the spot as a function argument without putting it in a variable. for_each uses this style (Section 4).
2.2 Capture: [=] and [&]
The Fortran internal procedure scaled could use the parent procedure’s variable factor as is (host association). To use outside variables in a lambda expression, you specify how to take them in inside []. Taking in outside variables is called capture.
Notation |
How it captures |
|---|---|
|
Captures the outside variables used in the body as copies of the values at the time the lambda expression was created |
|
Captures the outside variables used in the body as references (Chapter 3) |
|
Captures nothing. Using an outside local variable is a compile error |
The difference shows when you change an outside variable after creating the lambda expression.
#include <iostream>
int main() {
double factor = 2.0;
auto by_value = [=](double x) { return factor * x; }; // copy the value of factor (2.0)
auto by_ref = [&](double x) { return factor * x; }; // reference factor itself
factor = 10.0; // change the outer variable after creation
std::cout << by_value(3.0) << " " << by_ref(3.0) << "\n";
}
6 30
by_value keeps holding the 2.0 from when it was created, so the result is 6, and by_ref reads the changed factor, so the result is 30. The host association of a Fortran internal procedure behaves like [&].
A copy captured with [=] cannot be rewritten inside the body.
int count = 0;
auto inc = [=]() { count += 1; }; // compile error: read-only variable 'count'
When you want to update the outside variable itself, use [&]. However, MISO’s parallel loops use [=] (Section 3). The reason you can rewrite array elements with [=] is that what is captured is a View.
miso::Array1D<double> a(3);
auto av = a.view();
auto set_one = [=](int i) { av[i] = 1.0; }; // holds a copy of av (pointer and size)
set_one(2);
std::cout << a[2] << "\n"; // 1
The copied av points to the same elements as the original av, so av[i] = 1.0 rewrites the elements of a. What cannot be changed is av itself (the pointer and size). This is the relationship between const and T & in Section 1.3.
Exercise 7-2: Give the two numbers that the following program prints.
#include <iostream>
int main() {
int n = 1;
auto f = [=]() { return n * 10; };
auto g = [&]() { return n * 10; };
n = 2;
std::cout << f() << " " << g() << "\n";
}
It prints 10 20. f holds a copy of n (1) from when it was created, so 10, and g holds a reference to n and reads n (2) at the time of the call, so 20. A lambda expression with no arguments is also called by adding ().
2.3 A Small Struct That Holds Values
If you think of a lambda expression as a small struct that the compiler creates automatically, you can trace its behavior accurately. [=](double x) { return factor * x; } works almost the same as the following struct.
#include <iostream>
// struct that behaves almost like [=](double x) { return factor * x; }
struct ScaleByValue {
double factor; // copy of the value at creation
double operator()(double x) const { return factor * x; }
};
// struct that behaves almost like [&](double x) { return factor * x; }
struct ScaleByRef {
double &factor; // reference to the outer variable
double operator()(double x) const { return factor * x; }
};
int main() {
double factor = 2.0;
ScaleByValue by_value{factor}; // initialize member factor with 2.0
ScaleByRef by_ref{factor}; // bind member factor to the outer factor
factor = 10.0;
std::cout << by_value(3.0) << " " << by_ref(3.0) << "\n"; // 6 30
}
As you learned in Chapter 4, a struct with operator() can be called like a function, as in by_value(3.0). The correspondence with the lambda expression is as follows.
Lambda expression |
Struct |
|---|---|
Captured variables |
Member variables (values for |
Arguments and body |
Arguments and body of |
The place where the lambda expression is written |
Create an object, initializing the members with the values at that point |
Copies from |
|
This model tells us the following.
When you copy a lambda expression, the members it holds are also copied. For a lambda expression that holds a View with
[=], the pointer and size are copied.A lambda expression made with
[&]holds the addresses of outside variables. Calling it after the outside variable is gone gives the dangling reference of Chapter 3.[=]captures only the variables actually used in the body. It does not copy all the variables in the function.
When computing on the GPU, this lambda expression object is copied to the GPU. The GPU cannot read CPU memory, so a lambda expression that holds addresses on the CPU, such as one made with [&], cannot be used on the GPU. The details are explained in Chapter 9.
2.4 Member Variables and this
Consider writing a lambda expression inside a member function of a class and using member variables in its body. We said [=] makes “copies of values,” but member variables are not copied. Inside a member function, the member variable gm means “the gm of the object itself,” that is, shorthand for this->gm. this is a pointer to the object that called that member function. What [=] captures is a copy of this pointer this.
We check with the following example. Gas has the ratio of specific heats gm as a member, and creates and returns a lambda expression that computes pressure. If you make the return type auto, the compiler decides the type from the return expression.
#include <iostream>
struct Gas {
double gm; // ratio of specific heats (member variable)
// create and return a lambda that computes pressure
auto pressure_function() const {
return [=](double ro, double ei) { return (gm - 1.0) * ro * ei; };
}
// copy the member variable to a local variable before using it
auto pressure_function_copy() const {
const double gm_ = gm;
return [=](double ro, double ei) { return (gm_ - 1.0) * ro * ei; };
}
};
int main() {
Gas gas{1.4};
auto f = gas.pressure_function();
auto g = gas.pressure_function_copy();
std::cout << f(1.0, 2.5) << " " << g(1.0, 2.5) << "\n";
gas.gm = 2.0; // change the member variable after creating the lambda
std::cout << f(1.0, 2.5) << " " << g(1.0, 2.5) << "\n";
}
1 1
2.5 1
f holds not the value of gm but this (the address of gas), so changing gas.gm changes the result. g holds a copy of the value 1.4 of the local variable gm_, so the result does not change. In the struct model, f is a struct with the member const Gas *self;, and g is a struct with the member double gm_;.
When computing only on the CPU, this difference is merely a matter of “which one sees a value changed later.” On the GPU, the problem is that this points to CPU memory. For that reason, MISO copies member variables into local variables just before the lambda expression and then uses them. The reason, and what happens if you do not follow this, are explained in Chapter 9.
In the real MISO code, gas_pressure of the ideal-gas equation of state eos::IdealEOS has this form.
// excerpt from IdealEOS in miso/include/miso/eos.hpp
template <typename Real> struct IdealEOS {
/// @brief Ratio of specific heats (gamma)
Real gm;
...
template <typename Backend>
void gas_pressure(Backend btag, mhd::FieldsView<const Real> qq,
Array3DView<Real> pr) const {
const Real gm_ = gm; // device lambdas must not capture `this`
Range1D range{0, qq.size()};
for_each(
btag, range,
MISO_LAMBDA(int i) { pr[i] = (gm_ - Real(1)) * qq.ro[i] * qq.ei[i]; });
}
The comment “device lambdas must not capture this” means “lambda expressions that run on the GPU must not capture this.” for_each, Range1D, and MISO_LAMBDA are explained in order in the following sections.
Exercise 7-3: Suppose you delete const Real gm_ = gm; from gas_pressure and use gm directly inside the lambda expression. Explain what value the lambda expression then holds.
gm inside the lambda expression is shorthand for this->gm, so [=] copies and holds not the value of gm but the pointer this to the IdealEOS object. Each time the body of the lambda expression is called, it reads gm from the object that this points to. pr and qq are function arguments (local variables), so their values (Views) are copied as before. On the CPU the computed result does not change, but on the GPU it becomes a problem because this points to CPU memory (Chapter 9).
2.5 Passing a Lambda to a Function
The type of a lambda expression object has no name, so the function that receives it is written as a function template from Chapter 6.
template <class F> void apply_to_all(int n, F f) {
for (int i = 0; i < n; ++i) {
f(i); // call the given lambda with i
}
}
// caller side
apply_to_all(5, [=](int i) { av[i] = 0.5 * i; });
F is decided by the compiler as the type of the lambda expression passed in. The argument F f is passed by value, so the lambda expression object (including the Views it holds) is copied and passed. The compiler creates apply_to_all for each type of lambda expression passed. In builds with optimization enabled, the contents of f(i) are expanded inside the loop, so almost no function call overhead remains. MISO’s for_each, which we see in the next two sections, has exactly this form.
3. MISO_LAMBDA
In MISO, the lambda expressions passed to for_each and others are written MISO_LAMBDA instead of [=]. MISO_LAMBDA is a macro defined in miso/include/miso/execution.hpp (the #define of Chapter 2).
// excerpt from miso/include/miso/execution.hpp
#ifdef __CUDACC__
#define MISO_LAMBDA [=] __host__ __device__
#else
#define MISO_LAMBDA [=]
#endif // __CUDACC__
__CUDACC__ is a name defined only when the CUDA compiler nvcc compiles. In CPU builds compiled with g++, MISO_LAMBDA(int i) { ... } is replaced with [=](int i) { ... }. In other words, in CPU builds, MISO_LAMBDA is exactly the value capture [=].
In CUDA builds, __host__ __device__ is added. This is a CUDA specification meaning “this function can run on both the CPU (host) and the GPU (device)” (Chapter 9). In CPU builds, miso/include/miso/cuda_compat.hpp defines __host__ and __device__ as empty strings.
// excerpt from miso/include/miso/cuda_compat.hpp
#ifndef __device__
#define __device__
#endif
#ifndef __host__
#define __host__
#endif
For this reason, functions in MISO’s headers such as __host__ __device__ int size() const are also compiled as ordinary functions in CPU builds. When you write a function to be called from inside a lambda expression yourself, add __host__ __device__ if you plan to use it on the GPU later. For math functions, use miso::util::sqrt, util::exp, util::sin, and so on instead of std::sqrt (miso/include/miso/utility.hpp). The util functions work on both CPU and GPU, and accept both float and double.
We use MISO_LAMBDA so that the same source is compiled for both CPU and GPU. We use [=] and not [&] because, as seen in Section 2.3, lambda expressions that hold addresses on the CPU cannot be used on the GPU.
MISO’s unit test miso/tests/serial/test_execution.cpp uses [&] in some places in CPU-only tests. This also works on the CPU (backend::Host), but in code shared with the GPU, use MISO_LAMBDA.
4. Parallel Loops with for_each
Range1D and Range3D, which represent ranges
The range of indices passed to for_each is represented by Range1D and Range3D (miso/include/miso/execution.hpp).
// excerpt from miso/include/miso/execution.hpp (member functions omitted)
struct Range1D {
int begin, end;
};
struct Range3D {
Range1D range0, range1, range2;
};
Range1D is a struct with a start begin and an end end. It includes begin and does not include end (the same as for (int i = 0; i < n; ++i) in Chapter 2). A range that includes the start and excludes the end is called a half-open interval, and is sometimes written [begin, end). Range3D is the Range1Ds for the three directions lined up.
miso::Range1D r1{0, n}; // i = 0, 1, ..., n-1
miso::Range3D r3{{0, nx}, {0, ny}, {0, nz}}; // ranges in the 3 directions
Range1D{0, n} is the notation that initializes the struct’s members in declaration order, with begin = 0 and end = n. It corresponds to Fortran’s do i = 0, n - 1. Note that the upper limit is written differently by 1.
The Host version of for_each
for_each receives a range and the processing for each point, and calls the processing at every point in the range. The definition of the Host version is just a for statement.
// excerpt from miso/include/miso/execution.hpp
template <class F> inline void for_each(backend::Host, Range1D range, F f) {
for (int i = range.begin; i < range.end; ++i) {
f(i);
}
}
template <class F> inline void for_each(backend::Host, Range3D range, F f) {
for (int i0 = range.range0.begin; i0 < range.range0.end; ++i0) {
for (int i1 = range.range1.begin; i1 < range.range1.end; ++i1) {
for (int i2 = range.range2.begin; i2 < range.range2.end; ++i2) {
f(i0, i1, i2);
}
}
}
}
It has the same form as apply_to_all in Section 2.5. The first argument backend::Host is the tag type from Chapter 6 and has no argument name. Its contents are not used, and only its type selects between the Host and CUDA versions of for_each. The Host version runs in order on a single CPU thread. The three-dimensional version loops over the last index i2 innermost, to match the memory layout of Section 1.2.
The caller writes as follows.
#include <iostream>
#include <miso/array3d.hpp>
#include <miso/execution.hpp>
using namespace miso;
int main() {
const int nx = 2, ny = 3, nz = 4;
Array3D<double, backend::Host> ro(nx, ny, nz);
auto rov = ro.view(); // pass the View to the lambda
Range3D range{{0, nx}, {0, ny}, {0, nz}};
for_each(backend::Host{}, range, MISO_LAMBDA(int i, int j, int k) {
rov(i, j, k) = 100 * i + 10 * j + k;
});
std::cout << rov(1, 2, 3) << "\n"; // 123
}
Part |
Role |
|---|---|
|
The place to compute. Create a tag-type object with |
|
The ranges in the three directions. None of the directions includes its upper limit |
|
The processing to run at each point. For |
The lambda expression holds rov (a View) by value. ro (an Array3D) cannot be copied, so it is not used directly inside the lambda expression. Written in Fortran, it corresponds to the following computation.
do i = 0, nx - 1
do j = 0, ny - 1
do k = 0, nz - 1
ro(i,j,k) = 100*i + 10*j + k
enddo
enddo
enddo
Do not depend on the execution order
The Host version runs in index order, but the CUDA version runs on many threads at the same time, and the order is not fixed (Chapter 9). Write the processing passed to for_each so that the result does not change even when the points are executed simultaneously. Specifically, follow these two rules.
Do not read, in the processing of one point, an element written by the processing of another point. Separate the arrays you read from the arrays you write.
Do not write to the same element from the processing of different points.
For example, av[i] = av[i - 1] + 1.0 reads the result of i - 1 at i, so it works on the Host version but gives different results on the GPU.
A MISO example: gas_pressure
Let us read eos::IdealEOS::gas_pressure from Section 2.4 once more.
// excerpt from miso/include/miso/eos.hpp
const Real gm_ = gm; // device lambdas must not capture `this`
Range1D range{0, qq.size()};
for_each(
btag, range,
MISO_LAMBDA(int i) { pr[i] = (gm_ - Real(1)) * qq.ro[i] * qq.ei[i]; });
btagis a function argument and is eitherbackend::Host{}orbackend::CUDA{}.gas_pressureis a function template withBackendas a type argument, so the same code is used on both CPU and GPU.qqis anmhd::FieldsViewbundling the Views of the nine physical variables, andpris anArray3DViewto which the pressure is written.Although the arrays are three-dimensional,
Range1D{0, qq.size()}makes it a one-dimensional loop over the sequential indexi. The pressure at each point is determined only by the densityroand internal energyeiat that point, so there is no need to split it into three indices. The range is all elements, including ghost cells (Chapter 8).Real(1)is 1 of typeReal. It is a way of writing that avoids mixing in adoubleconstant whenRealisfloat.
The processing at each point writes to pr[i] and only reads qq.ro[i] and qq.ei[i], so the result is the same in whatever order it runs.
Exercise 7-4: Write the following Fortran loop using a and b (size n) of Array1D<double>, for_each, and MISO_LAMBDA. Assume the arrays are declared as a(0:n-1) and b(0:n-1).
do i = 1, n - 2
b(i) = 0.5d0 * (a(i-1) + a(i+1))
enddo
Also explain why, if you change this processing to write to a(i) instead of b, for_each cannot compute it correctly.
auto av = a.view();
auto bv = b.view();
for_each(backend::Host{}, Range1D{1, n - 1}, MISO_LAMBDA(int i) {
bv[i] = 0.5 * (av[i - 1] + av[i + 1]);
});
Fortran’s do i = 1, n - 2 includes n - 2, so the upper limit of Range1D is n - 1, which is larger by 1. What you pass to the lambda expression are the Views av and bv. If you write to the same array, as in a(i) = 0.5 * (a(i-1) + a(i+1)), then when the processing of point i reads a(i-1), whether the processing of point i-1 has already rewritten it depends on the execution order. The Host version runs in order, so it gives the same result as Fortran, but when run simultaneously, as on a GPU, the result changes. Separate the array you read the pre-update values from and the array you write the result to.
5. Reduction with reduce
What is a reduction
An operation that finds a single value, such as a sum, maximum, or minimum, from all the elements of an array is called a reduction. It corresponds to Fortran’s sum, maxval, and minval. MISO uses reduce (miso/include/miso/execution.hpp).
// excerpt from miso/include/miso/execution.hpp
template <class T, class F, class Op>
inline T reduce(backend::Host, Range1D range, T init, F map, Op op) {
T acc = init;
for (int i = range.begin; i < range.end; ++i)
acc = op(acc, map(i));
return acc;
}
Argument |
Role |
|---|---|
|
The place to compute |
|
The range. There are versions for |
|
The initial value. The result type |
|
A lambda expression that returns the value at each point. |
|
A lambda expression that combines two values into one. For a sum, |
Unlike for_each, map does not write to arrays, and returns the value at that point with return. The Host version just repeats acc = op(acc, map(i)) in order. You receive the result as the return value of reduce.
#include <iostream>
#include <miso/array1d.hpp>
#include <miso/execution.hpp>
using namespace miso;
int main() {
const int n = 5;
Array1D<double, backend::Host> a(n);
auto av = a.view();
for_each(backend::Host{}, Range1D{0, n}, MISO_LAMBDA(int i) {
av[i] = 0.5 * i; // 0, 0.5, 1, 1.5, 2
});
const auto value = MISO_LAMBDA(int i) { return av[i]; }; // value at each point
const auto plus = MISO_LAMBDA(double s, double t) { return s + t; }; // combine two values
const double total = reduce(backend::Host{}, Range1D{0, n}, 0.0, value, plus);
std::cout << total << "\n"; // 5
const double wrong = reduce(backend::Host{}, Range1D{0, n}, 100.0, value, plus);
std::cout << wrong << "\n"; // 105 (the initial value is also added)
}
The lambda expressions are put in the variables value and plus before being passed. You may also write them directly inside the call to reduce.
Make the initial value the identity element
In the second call in the example above, the initial value 100.0 was added to the result, giving 105. The Host version uses the initial value only once, but the CUDA version uses the initial value many times, as the values of threads outside the range (reduce_stage_kernel in execution.hpp). For that reason, make the initial value something that does not change the value when combined (the identity element).
What you want |
|
Initial value |
|---|---|---|
Sum |
|
|
Minimum |
The smaller one |
A sufficiently large value (such as |
Maximum |
The larger one |
A sufficiently small value (such as |
std::numeric_limits is a standard library feature, available with <limits>, that returns a type’s maximum and minimum values and so on. lowest() is the smallest value (negative, with the largest absolute value), and max() is the largest value.
The result type T is determined by the type of the initial value. If you write the initial value of a sum as 0 rather than 0.0, T becomes int, and each partial sum is truncated to an integer. Trying it with the array above returns 4 instead of 5. For reductions of floating-point numbers, write the initial value like 0.0.
For functions that combine maxima and minima, utility.hpp has util::max2 and util::min2. In MISO’s main as of October 2026, these round to float before comparing even when given double (MISO issue #189). For double arrays, use util::fmax and util::fmin, or write it yourself, such as a < b ? b : a.
The sum changes with the order of addition
The sum of floating-point numbers changes in the trailing digits depending on the order of addition. The Host version adds in index order, but the CUDA version makes partial sums for each group of threads and then adds them together. For that reason, the bit patterns of the sums on CPU and GPU do not necessarily match. When comparing results, distinguish whether you require an exact match or use a tolerance. Maxima and minima are chosen just by comparison, so they match if the arrays are the same.
A MISO example: Integrator::cfl
The time step dt is found from the CFL condition as the minimum over all grid points. The CUDA version of Integrator::cfl uses reduce. Even without a GPU, you can read it as an example of how to use reduce.
// excerpt from Integrator::cfl in miso/include/miso/mhd_integrator_cuda.cuh
auto qq_v = qq.const_view();
auto grid_v = grid.const_view();
auto cs_v = cs.view();
// device lambdas must not capture `this`
const Real cfl_number_ = cfl_number;
eos.sound_speed(backend::CUDA{}, qq_v, cs_v);
Range3D range{{grid.i_margin, grid.i_total - grid.i_margin},
{grid.j_margin, grid.j_total - grid.j_margin},
{grid.k_margin, grid.k_total - grid.k_margin}};
const auto f = MISO_LAMBDA(int i, int j, int k) {
Real vv = util::sqrt(qq_v.vx(i, j, k) * qq_v.vx(i, j, k) +
qq_v.vy(i, j, k) * qq_v.vy(i, j, k) +
qq_v.vz(i, j, k) * qq_v.vz(i, j, k));
Real ca = util::sqrt((qq_v.bx(i, j, k) * qq_v.bx(i, j, k) +
qq_v.by(i, j, k) * qq_v.by(i, j, k) +
qq_v.bz(i, j, k) * qq_v.bz(i, j, k)) /
qq_v.ro(i, j, k) * pii4<Real>);
Real total_vel = cs_v(i, j, k) + vv + ca;
Real dxyz = util::min3(grid_v.dx[i], grid_v.dy[j], grid_v.dz[k]);
return cfl_number_ * dxyz / total_vel;
};
const auto op = MISO_LAMBDA(Real a, Real b) { return util::min2(a, b); };
const Real dt_max = 1.e10;
const auto dt = reduce(backend::CUDA{}, range, dt_max, f, op, cfl_helper);
Before the lambda expression, it prepares the Views of the physical variables, grid, and sound speed (
qq_v,grid_v,cs_v) and the copycfl_number_of the member variablecfl_number(Section 2.4).The range is a
Range3Dcovering the interior grid points, excluding ghost cells (i_marginis explained in Chapter 8).freturns the “time step allowed at that point” for each point, andopreturns the smaller one. The initial value1.e10is used in place of the identity element, as a value sufficiently larger than any actual time step.The last
cfl_helperis a workspace used only by the CUDA version (Chapter 9).pii4<Real>is the variable template from Chapter 6.
The Host version of Integrator::cfl (miso/include/miso/mhd_integrator_host.hpp) does not use reduce, and writes the same computation with a triple for statement over i, j, and k and util::min2. In both, the minimum over all processes is finally taken with MPI_Allreduce (Chapter 8).
Exercise 7-5: Write code that finds the minimum of an Array1D<double> array a (size n) with reduce. Also explain what you choose for the initial value and why.
Since it uses std::numeric_limits, include <limits>.
auto av = a.view();
const double a_min = reduce(
backend::Host{}, Range1D{0, n}, std::numeric_limits<double>::max(),
MISO_LAMBDA(int i) { return av[i]; },
MISO_LAMBDA(double s, double t) { return s < t ? s : t; });
Make the initial value a value that is never chosen when compared with any element. If you use the largest double value, std::numeric_limits<double>::max(), the array’s value always remains when compared with the array’s values. If you use 0.0 as the initial value, the result becomes 0 when all the array’s values are positive. util::min2 rounds double to float (MISO #189), so here we pick the smaller one with the conditional operator.
6. Writing a Small Program
In Practice: Arrays and Parallel Loops, you work on vector addition, reduction, and one-dimensional advection, in that order. For the CMake setup, use “the minimal CMake setup for using MISO from your own program” in Chapter 1. At this stage, run with backend::Host. Running the same assignments on the GPU comes after reading Chapter 9.
Once you have checked array computations on the CPU, move on to domain decomposition with multiple processes in Chapter 8. If you would like to write the time evolution first, you can also work on the serial computation in Chapter 10.
::In Fortran code, OpenMP is often used, splitting loops among multiple CPU threads with !$omp parallel do. The Host version of MISO’s for_each is just a for statement, as seen in Section 4, and does not use OpenMP. MISO’s source (miso/ and demo/) has no #pragma omp, and CMake does not look for OpenMP either.
MISO’s parallelization is handled by domain decomposition with MPI on the CPU (Chapter 8) and by CUDA on the GPU (Chapter 9). To use multiple cores on the CPU, launch as many MPI processes as there are cores.
Even so, the advantage of writing the Host version with for_each is that you can switch to the CUDA version with the same source. If you added OpenMP to the Host version, you would change only the contents of for_each and reduce in execution.hpp, and the calling side would not change. This is an inference from MISO’s design, not a policy written in MISO’s documentation.