BOOST Library
This document gives a practical introduction to the BOOST Libraries oriented to common tasks in scientific computing. For those using C++ as a programming language in general, BOOST is a very useful set of libraries that help programming in C++ by providing tools and commonly used idioms that solve recurring problems that emerge when developing in this language, even at a fairly low level. It also encourages certain style of programming, where generality and efficiency are priority. These tools include macro definitions, template functions and classes, functions, generic data containers, and common algorithms. It can be viewed as a natural extension of the C++ Standard Template Library (STL).
General Remarks
BOOST, like C++, has a very steep learning curve and the syntax can become sometimes very ugly and complicated, but it faster to learn BOOST than to reinvent the wheel and come with buggy home-brew solutions that at the end of the day will look even more ugly and won't be that good anyway. (By good I mean generic and efficient.) Note that I am not comparing BOOST and C++ with other tools like Matlab, or Fortran, or Python. Those language are likely to do better (and quick development) than C++ in many specific areas. Boost also closes the gap with those languages without introducing modifications to the core C++ language.
Libraries included in Boost range from very simple ones (like Boost.Conversion) to very complex ones (like Boost.MPI and Boost.BGL).
The aim of this document/tutorial is to incrementally document the utility and working practical examples. Extra documentation added by users is very welcome, specially in the area of non-trivial but concise examples, related to our common programming problems (i.e. scientific computing).
Topics covered in the first stages of this document will include: Boost building and installation, [Boost.MultiArray] and Boost.MPI. That are the libraries that I have to deal with at this point. Ideally I would like to add [Boost.uBlas] in the future and document examples of small libraries that are very usefull for numerical tasks. The focus will be to document the interaction of these libraries with numerical libraries such as FFTW and Lapack, and above all, the documentation of practical example code.
BOOST build and installation
It is a bless that most Linux distributions now include some versions of Boost. Many Boost libraries are header-only, which means that no linking is necessary (an #include<boost/ ... > line set us ready to use a particular library). For the rest of the libraries, linking is not hard either, it is just a matter of adding something like "-lboost_NAME" during compilationg (with certain naming convention). The current tutorial has been tested in wcr.stanford.edu.
The difficulty with linking is that most Boost distributions in the operating system are outdated; for example the current version of Boost is 1.37, while the most recent distributions come with version 1.34 at best. Once in a while we will need a library not included in older releases (for example Boost.MPI). For "header-only" libraries the situation is not too bad because the header files can be downloaded and they ready to use. For binary libraries we need to compile them. There are tools to compile specific libraries but it will be easier for the moment to compile the whole Boost set.
As usual, to be kind with our usage in administrated system, we will install Boost in our user space. for example in $HOME/user.
mkdir $HOME/usr
the directories ~/usr/include/boost and ~/usr/lib will be created and populated after the building of Boost. Be brave and install the last current version from the link bellow:
mkdir $HOME/soft cd $HOME/soft wget http://internap.dl.sourceforge.net/sourceforge/boost/boost_1_37_0.tar.gz
If the download fails, try to download the last version from Sourceforge Boost download page. Then decompress it:
tar -zxvf boost_1_37_0.tar.gz cd boost_1_37_0 ./configure --prefix=$HOME/usr
Append the following quoted line to user-config.jam to include Boost.MPI among the installed libraries.
echo "using mpi : mpicxx ;" >> user-config.jam
where mpicxx can be replaced by your mpi compiler wrapper, for example "/usr/lib/mpich/bin/mpiCC", or first selected interactively with "mpi-selector-menu" if is available. Then compile and install the package
make make install
Go lunch, that takes ~30 minutes. After successful compilation and installation, all the header files will be installed in ~/usr/include/boost-1-37/boost/*.hpp and the binary linkable files will be in ~/usr/lib/libboost_*:
Depending on how you want to compile your programs it might be necessary to add these directories to the include path of compilation for example, LD_LIBRARY_PATH.
Boost.MultiArray Library
The MultiArray library is the first to be described in this tutorial because it is the most simple of these libraries and yet it solves a very annoying problem with C/C++ for our area. C/C++ have a very limited support for built-in arrays and the situation is even worst for multidimensional arrays (2, 3 or more indices). Passing arrays to functions is very error prone, and functions with endless parameters with array dimensions must be continuously passed around. Moreover, writing algorithms for arbitrary array dimensionality is very hard.
We have to understand that the reason for this is that the language was designed to be very flexible, there are infinite ways to arrange multidimensional arrays in memory and infinite ways to resize them, so the language did not enforce any particular standard. Yet for scientific computing (which is only a subset of all computing areas), we can agree that there are a few multidimensional arrays that make sense in numerical problems, namely contiguous blocks memory with either column-major ordering (C-convention) or row-major ordering (Fortran-convention).
Among this constrain of the memory layout, MultiArray provides a very flexible, very small and nice interface. The official MultiArray manual provides rather simple but useful examples. Below there are more examples more suitable for existing numerical programs.
From the simple examples it can be noticed MultiArray usage is very similar to the default C-arrays with few very convenient differences 1) no explicit allocation (ie. malloc/new) is needed 2) no explict deallocation (free/delete[]) 3) indices ordering can be C-like or Fortran-like (useful for interfacing with Fortran routines) 4) Arrays keep track of their shape (no separate variables with size are needed).
The typing system of MultiArray is quite complicated and confusing, one way to alleviate this problem is to make extensive use of typedefs. (You may seldom benefit from using Template typedefs.)
Using MultiArray on existing code together with C-arrays
MultiArray can be adapted easily in programs that already use raw-arrays since the first element to the contiguous array in memory can be known for functions that expose the pointer to the array in memory; the shape of the array can be also accessed:
boost::multi_array<double, 2> A(boost::extents[10][10]); // allocates memory for 100 elements ... double* A_ptr = A.origin(); int A_N0 = array.shape()[0]; // A_N0 <- 10 int A_N1 = array.shape()[1]; // A_N1 <- 10 ... // use the array in a typical c-function somecfunction(A_N0, A_N1, A_ptr); // or as cfunction(A.shape()[0], A.shape()[1], A.origin() ); ... // .. no need to free memory, it is deallocated when A goes out of scope
Such function can even modify the elements of the array. In the example the array memory is managed by the multi_array A and therefore the memory is freed when A goes out of scope. Sometimes, for existing applications, this is not possible because the memory is managed in some other complicated way. In that case we can still take advantage of multi_array by using the "multi_array_ref" class, that is the same as an array but does not manage the underlaying memory:
int A_N0 = 10; int A_N1=10; double* A_arr = new double[N1*N2]; // or = malloc(sizeof(double)*A_N0*A_N1); ... boost::multi_array_ref<double, 2> A_ref(A_arr, boost::extents[A_N0][A_N1]); ... // access A_arr through A_ref or call a MultiArray aware function somemultiarrayfunction(A_view); ... delete[] A_arr; //or free(A_arr); (do not use A_view after this point!)
In this way we can go back and forth form C-arrays (raw-arrays) to boost::multi_array. Generalizations to higher dimensions are similarly easy. (One of the strong points of MultiArray is that dimensionality can be arbitrary). Let me stress that in the first case the memory is automatically managed for us, while in the second case we have to take care ourselves.
A simple FFTW example
In the context of numerical libraries the following is an example to call FFTW3 in combination with MultiArray. (Remember to Install FFTW3 first.)
boost::multi_array<std::complex<double>, 2> A(boost::extents[10][10]);
for(int i=0; i!=10; i++)
for(int j=0; j!=10; j++)
A[i][j]= std::complex<double>(i,j);
...
fftw_plan p=fftw_plan_dft_2d(A.shape()[0], A.shape()[1], A.origin(), A.origin(), FFTW_FORWARD, FFTW_ESTIMATE); // in place
fftw_execute(p);
fftw_destroy_plan(p);
Advanced interface of FFTW take stride arguments as well, this is not a limitation for MultiArray since, for example, array_view's (see below) also keep track of internal strides (e.g. A_center.stride()[n] is the stride in direction n of the referenced memory). The offset is given be A_center.origin()-A_center.data(). (origin() always refer to the first element of the (sub)array(view), while data() can return a pointer that does not belong to the data but is used for internal reference. For plain multi_array's (with 0-base) origin() and data() are the same.
Fortran ordering and a simple Lapack-type binding
Calling Lapack using array's makes sense only for one or two dimensional cases. Since Lapack expects arrays in column-major order format, it might be a good point to introduce this option, which is supported by Boost.MultiArray. (For linear algebra matrices it is more recommendable to use Boost.uBlas which is more oriented to linear algebra operations and one or two dimensional arrays, i.e. vectors and matrices.)
The ordering of the elements in memory is controlled by an option that can be set when declaring the array.
boost::multi_array<double,2> A(boost::extents[10][10], boost::fortran_storage_order); ... // assign A value with Fortran index ordering
There is even an option to declare the multi_array to be 1-based (instead as 0-based); details can be found here. Having said that, we can expand in the difference between .origin() and .data(). While .origin() always refer to the first element (whatever this one is, e.g. A[1][1]...[1]), .data() always refer to A[0][0]...[0] regardless of the existence of the element as part of the allocated memory.
We then can call a Fortran function as:
int A_N1=A.shape()[0]; int A_N2=A.shape()[1]; fortran_function(&A_N1, &A_N2, A.origin()); // Fortran arguments are passed by reference(pointer)
or write a wrapper to the Fortran function
wrapper_function(boost::multi_array<double, 2>& in){
int in_N1=in.shape()[0]; int in_N2=A.shape()[1];
fortran_function(&in_N1, &in_N2, in.origin()); // Fortran arguments are passed by reference(pointer)
}
and then call it simply by doing 'wraper_function(A);'.
Calling a Fortran function that takes stride arguments can be done by means of subarray's and this case is not covered here.
Subarrays and views
One of the useful properties of raw-arrays is that regular subarrays can be defined simply by defining offsets and strides. MultiArray provides this kind of creation of subarray in a quite elegant manner. For example given the array defined in the example above, the most simple subarray is the object A[5] which is automatically a one-dimensional array containing a reference to the sixth (indices start at 0) column of the array, so for example '(A[5])[1]=0' changes the value of element A[5][1]. Of course that there are more complicated subarray views that can be defined in terms the original array. The basic examples are given in the official MultiArray tutorial. The following are some digested examples, which at the same time introduces 'typedef' statements to reduce verbosity:
Subarrays
typedef std::complex<double> complex; typedef boost::multi_array<complex,2> array2d; using boost::extents; array2d A(extents[10][10]); // 2D, 100 allocated elements typename array2d::sub_array<1>::type A_row5 = A[5]; // 1D, 10 element (no copy)
Views (array_view)
typedef array2d::array_view<1>::type view1d; typedef boost::multi_array_types::index_range range; boost::multi_array_types::index_gen indices; view1d A_col3_5to10 = A[3][indices[range(5,10)]]; // or A[indices[3][range(5,10)] view1d A_col3_5to10odd = A[indices[3][range(5,10,2)]; typedef array2d::array_view<2>::type view2d; view2d A_center = A[indices[range(2,7)][range(2,7)]]; // 2D, 25 elements view2d A_upperleft_quadrant_oddonly = A[indices[range(0,5,2)][range(0,5,2)]]; // 2D, 4 (2x2 ) elements view2d A_lowerright_quadrant_evenonly = A(range(5,10,2),range(5,10,2)); // 2D, 4 (2x2 ) elements
Last two arrays refer to odd and even elements of the original array.
Also note that something like
view1d A_error = A[indices[range(2,7)][range(2,7)]];
is not allowed because the RHS is 2-dimensional and the LHS is declared to be 1-dimensional.
None of the array defined here make any copy of elements, remember that sub_arrays and array_view are just references, another way to view the original array.
Custom index ranges
In previous examples we used some obscure entity called boost::extents, it is obviously used to pass the valid range of index to the matrix declaration.
Technically extents is an (empty) global variable of type boost::extent_gen; this is just a trick to allow a nice syntax with arbitrary number of parameters, like passing extents[10][10][10][10] for a 4D-array.
The real power of this syntax is that we can really pass other types than just integers declaring the size in each dimension, namely we can declare ranges instead. For example, this will create a 1-based array of size 10x10.
using boost::extents; using boost::range; boost::multi_array<double, 2> A(extents[range(1,11)][range(1,11)]);
In this case the elements A[0][0], A[0][1], A[1][0] are out-of-bounds, just like A[11][11]. Furthermore, we can create an array that is behaves completely like a Fortran-array by also indicating the ordering.
boost::multi_array<double, 2> B_Flike(extents[range(1,11)][range(1,11)], boost::fortran_storage_order);
Let also mention that we can reindex an array after creation by calling the .reindex method:
A.reindex(0); // now all dimensions are 0-based, valid elements 0..9 x 0..9 ... // take advantage of 0-based array A.reindex(1); // 0-based index restored, valid elements 1..10 x 1..10
(Reindixing is just a new logical rearrangement of element referencing, not copies of the matrix or the elements is made.)
The indices can start an any arbitrary point (not just 0 or 1), even negative values -- in this example (taken from here), the indexing is totally arbitrary in each dimension:
typedef boost::multi_array<double, 3> array3D; array3D A(extents[8][range(1,4)][range(-1,3)][range(10,20)]);
(Note that, as an argument, [8] is the same as [range(0,8)].)
A simple FFTW wrapper
In my opinion, one of the problems of programming with C++ is that, since the core language allows only very basic constructs, we tend to repeat long code blocks even for doing basic things. What I propose here is to show that it is very easy to build a wrapper of known libraries that will make the syntax of the main program very easy to understand by packaging common tasks. The packaging of common tasks is traditionally achieved by writing simple C-like functions. The problem is that we sometimes need more flexibility such as C++ objects. Naturally this will need a more deep understanding of the language and of the wrapped libraries.
For example,in the case of FFTW, although written in pure-C, it is obvious that 'fftw_plan' structure behaves almost as an object: 1) It has a certain state [pointers and size of arrays], 2) lifetime [until fftw_destroy_plan] 3) and can be referred recurrently during it lifetime [can be executed several times during it lifetime]. If we combine this with what we already know of Boost.MultiArray we can design a robust interface. This is a draft of the wrapper classes:
using boost::multi_array;
typedef std::complex<double> complex;
typedef multi_array<complex,1> array1d;
typedef multi_array<complex,2> array2d;
namespace fftw3{
class plan{
public:
plan(array1d const& in, array1d& out, int direction, unsigned flags){
assert(in.size()==out.size()); // check size compatibility
p_=fftw_create_plan_1d(in.size(), in.origin(), out.origin(), direction, flags);
}
... // more complicated constructors here, for example array_view will need the advanced interface of fftw
void execute() const{fftw_execute(p_);}
virtual ~fftw3(){fftw_destroy_plan(p_);}
private:
fftw_plan p_;
};
}
The we can use it in the program:
boost::multi_array<std::complex<double>,2> A(boost::extents[100][100]); boost::multi_array<std::complex<double>,2> B(boost::extents[100][100]); // must have the same shape! fttw3::plan p(A,B, FFTW_FORWARD, FFTW_ESTIMATE); ... // assign values to A p.execute(); // repeated times ... // no need to explicitly destroy plan
The wrapper is elegant in some aspects, 1) it is uniform for all input/output dimensionalities and for in-place/out-of-place transform, 2) special cases where array_views are passed can be handled by the same syntax without having to call the advanced FFTW functions, 3) no explicit destroy needed (no risk to execute a destroyed plan).
Boost.MPI Library
I will assume that it is the case that know the very basics of C-MPI and you want to improve the readability and design of your code. If you don't know C-MPI, you may try to learn from this examples (the syntax is very nice to learn what mpi is supposed to do --in contrast with the ugly C version--) but be aware that Boost.MPI is not widely used.
I will try to work out an example that is not trivial by using FFTW in parallel. Hopefully this approach will be more related to our application area, while the official Boost.MPI pages are plenty of (deceptively) trivial examples.
Simple MPI FFTW example
In my opinion, using a real MPI numerical library is the only way to learn MPI and make it do useful work.
The example in this section uses the MPI version of FFTW3 which must be installed prior to begin. See detailed instructions to Intall FFTW3 first if necessary.
Before going to the first example let us review the command line to compile it
mpicxx -I$HOME/usr/include fft_mpi_test.cpp \
-L$HOME/usr/lib -lfftw3_mpi -lfftw3 \
-lboost_mpi-gcc43-mt-1_37 -lboost_serialization-gcc43-mt \
-o fft_mpi_test
besides the libfftw part, -L$HOME/usr/lib -lboost_mpi-gcc43-mt-1_37 -lboost_serialization-gcc43-mt do the job of linking against Boost.MPI and Boost.Serialization.
The reason we need Boost.Serialization (although we don't use it explicitly here) is that MPI is highly dependent on the serialization mechanism. Serialization is the capacity of converting any object or structure in memory into a linear stream of data. That is what we do over and over again when we save data to a file, in this case MPI communication is basically the same thing: one processes send streams of data to others. This integration allows to pass complex objects between processes without having to convert (deconstruct) the objects to numerical (or primitive) types first and then reconstructing the object on the other side as we would do it with C-MPI.
The following compilable example is based on the fftw3 simple mpi example (plain-C sources). In the following example we just add the Boost.MPI syntax (instead of the C-MPI):
#include <complex> /* must include before fftw3*.h */
#include <fftw3-mpi.h>
#include <boost/mpi.hpp>
#include <iostream>
namespace mpi = boost::mpi; /* if not defined have to use names as boost::mpi::... */
using std::clog;
using std::endl;
int main(int argc, char **argv){ /* initialize mpi world and fftw_mpi */
mpi::environment env(argc, argv); /* in C: MPI_Init(&argc, &argv); */
mpi::communicator world; /* this will work as MPI_COMM_WORLD */
fftw_mpi_init(); /* initialize fftw_mpi library */
const ptrdiff_t N0 = 18, N1 = 18; /* size of parallel array */
fftw_plan plan; /* this is the usual fftw plan */
std::complex<double> *data; /* in C: fftw_complex *data; this is the local (to process) data */
// given the matrix size (N0, N1) ask fftw_mpi how to split the matrix
ptrdiff_t alloc_local, local_n0, local_0_start, i, j;
alloc_local = fftw_mpi_local_size_2d(N0, N1, world,
&local_n0, &local_0_start); // alloc_local isn't always local_no*N1!!
//report allocation size, example of output from each process
if(world.rank()==0){ /* only process 0 prints this */
clog<<"Global size of array : [0:"<<N0<<")x[0:"<<N1<<"), total size "<<N0*N1<<endl;
}
//all processes print the following but each prints a different report (there are better ways to print)
clog<<"Process "<<world.rank()<<" has the array of size ["<<local_0_start<<","
<<local_0_start+local_n0<<")x[0,"<<N1<<"), total local size "<<(N0*N1)<<endl;
//now allocates data in the old fashioned way
data = (std::complex<double>*) fftw_malloc(sizeof(std::complex<double>)*alloc_local);
// in C: data = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * alloc_local);
//create fftw plan, note "_mpi_" and note the "_2d"
plan = fftw_mpi_plan_dft_2d(N0, N1, (double(*)[2])data, (double(*)[2])data, world,
FFTW_FORWARD, FFTW_ESTIMATE);
double pdata=0;
//initialize data and compute local power
for (i = 0; i < local_n0; ++i){
for (j = 0; j < N1; ++j){
std::complex<double> val=std::complex<double>(local_0_start + i, 0);
data[i*N1 + j] = val;
pdata+=norm(val);
}
}
clog<<"Process "<<world.rank() /* each process prints its local power*/
<<" has a power of "<<pdata<<" of the original data"<<endl;
//compute transform according to the plan
fftw_execute(plan);
double ptransform=0;
// normalize data and compute local power of the transform
for (i = 0; i < local_n0; ++i){
for (j = 0; j < N1; ++j){
data[i*N1+j]/=sqrt((double)(N0*N1)); // this is the normalization of the transform
std::complex<double> val=data[i*N1+j];
ptransform+=norm(val);
}
}
clog<<"Process "<<world.rank() //each process print its local power
<<" has a power of "<<ptransform<<" of the transformed data"<<endl;
fftw_destroy_plan(plan); /* destroy plan */
// boost::mpi::environment doesn't need finalization it is automatic, in C: MPI_Finalize();
}
The example is rather long, this is mainly because we are using the plain-C fftw interface, and manually allocated tables. Note that the only change relative to the original FFTW example is the usage of Boost.MPI and the use of C++ complex type. (Imagine how much cleaner would be the code if we could take advantage of boost.multiarray and of a C++ wrapper of FFTW. Eventually, this will be the topic of another section.)
Note that using C++ and Boost does not enforce us to leave the known C-interfaces, the idea is that we can reuse the old interfaces and improve them gradually as we need it. For example we used std::complex<double> instead of double[2], (this is possible because both types are bit-compatible); the former has a much support for treating complex numbers.
What the program does is to initialize parts of matrix in different processes, calculate local power (sum of squares) of the data, do a global Fourier transform of this matrix (in place) and the compute the local power of the transformed data. The global power (sum of locals powers) must be equal in the original data and in the transformed data. (See above for compilation command.)
Also note that although this is our first parallel MPI program the real MPI communication is performed by the FFTW routines and not by us. Rarely we will want to deal with passing very large amount of data ourselves between processes, libraries like FFTW and Lapack will do a faster/better job than us. The idea here is to set the environment for the usage of those well tunned numerical libraries.
Now let us download simple_boost_mpi_example.tar.gz and run the example:
wget http://micro.stanford.edu/mediawiki-1.11.0/images/Simple_boost_mpi_example.tar.gz -O simple_boost_mpi_example.tar.gz tar -zxvf simple_boost_mpi_example.tar.gz cd simple_boost_mpi_example make simple_boost_mpi_example
(Checkout the Makefile to see the compilation line). Run the example to see an output like the following:
$ mpirun -np 2 ./simple_boost_mpi_example Process 1 has the array of size [9,18)x[0,18), total local size = 162, local allocated size = 162 Global size of array : [0:18)x[0:18), total size 324 Process 0 has the array of size [0,9)x[0,18), total local size = 162, local allocated size = 162 Process 0 has a power of 3672 of the original data Process 0 has a power of 27729 of the transformed data Process 1 has a power of 28458 of the original data Process 1 has a power of 4401 of the transformed data
Two copies of the program are run at the same time, FFTW takes care of the communication between them. Before looking at any number, we can see that the output messages are not well ordered. we learned the first lesson of MPI in general: not synchronized instructions can be executed in any order. Besides this the program works as expected. The original 18x18 matrix is split among the two processors in two matrices of 9x18, and a memory for 162 elements is allocated in each process.
Also we can check that the power of the global data is the same before and after the FFT, 3672+28458 == 27729+4401 == 32130.
We can check running the program in one process.
$ mpirun -np 1 simple_boost_mpi_example Global size of array : [0:18)x[0:18), total size 324 Process 0 has the array of size [0,18)x[0,18), total local size = 324, local allocated size = 324 Process 0 has a power of 32130 of the original data Process 0 has a power of 32130 of the transformed data
Important notes about MPI-FFTW3 memory distribution (do not skip)
Although the distribution of the matrix and allocated memory seems to very reasonable in the previous example, this is only because of the simple dimensions and process number used in the example. There are some characteristics of the memory distribution for which the previous example can be misleading, and it is worth noting them:
- Only FFTW (via fftw_mpi_local_size_2d) knows how to split the data and allocated memory size among processes, so don't try to predict the splitting yourself.
- The splitting of the global matrix is only done over rows (first dimension, for N>2-dimensional arrays).
- 1D FFT's usage is totally different.
- The splitting of the matrix can be uneven (for example, #columns is not divisible by #processes).
- Some processes can be assigned zero columns (for example if there are more processes than columns).
- The logical size of the submatrix (local matrix) can be different from the local allocated memory, i.e. local_n0*N1 <= alloc_local. The outputs of local_size are not redundant. FFTW uses the (small) extra memory as a scratch space.
The next step in complexity could be to intercommunicate the processes in order to compute the total (global) power of the data.
Collective Reduce operation
There are several kinds of MPI operations, point-to-point and collective. Among the collective operations, the reduce operation is what we need to compute the global power of the data. Reduction means that we take *all* the objects/numbers from different processes and convert it into a *single* object/number by means of some function. More often than not, this function is simply a summation ().
In practice, a summation is a very simple function, because it is associative and commutative (sequence-product is as well), so many issues are simplified in this case because the result will not depend (mostly) on how the order in which the data is gathered. This kind of operations are defined only in terms of a binary function. That is, the summation is defined in terms of the binary addition alone. I said *mostly* because for machine numbers the usual operations are not *exactly* commutative, and the result can depend on how the data is split among processes. This is the second lesson of parallel programming, numeric results *can* depend on the number of processes. In fact, if the result completely changes with the number of processes then that may point to very deep conceptual errors in your numeric algorithm, since the result depend on how the data is summed (for example, when you sum small and large numbers together). This is very different from a 'bug' in the program (or in the MPI library).
In the C++ standard library, the binary addition is represented by a function object called std::plus<T> which represents the addition over type T, we are interested in the case in which T is the double type. This functions object have become the standard way of passing arbitrary functions to other functions. Function-object seem misterious sometimes but they are simply a representation of the operation, for example this code works as expected:
std::binary_function<double, double, double>* myop; myop=new std::plus<double>(); double plus_result=(*myop)(4,5); // or myop->operator()(4,5); assert(plus_result==9); delete myop; myop=new std::multiply<double>(); double multiply_result=(*myop)(3,5); assert(multiply_result==15); delete myop;
all_reduce
Going back to the global power example we can make all the processes receive the summation of values by inserting the lines:
double global_pdata=-1; mpi::all_reduce(world, pdata, global_pdata, std::plus<double>()); double global_ptransform=-1; mpi::all_reduce(world, ptransform, global_ptransform, std::plus<double>()); ... cout<<"Process "<<world.rank()<<" knows that global power of original data is "<< global_pdata << endl; cout<<"Process "<<world.rank()<<" knows that global power of transformed data is "<< global_ptransform << endl;
Now the program will compute the power of the data for us, which makes easier to verify that the Fourier transform is consistent. Note that the third argument is the variable which will be overwriten with the result and the second is the contribution from the current process to the reduction.
From this point on we can easily adapt the reduction to other needs. We can pass a different function, like std::multiply<double>() or, just as an example, selectively pass data if it fulfills certain condition:
mpi::all_reduce(world, pdata>0.?pdata:0, std::plus<double>()); // only pass positive values. mpi::all_reduce(world, world.rank()<4?pdata:0, std::plus<double>()); //only pass data of the first 4 processes.
Note that in both cases the condition of passing the data is *inside* the function call (via the ternary function condition?casetrue:casefalse). If you put the condition outside the function call, there exists the possibility of not fulfilling such condition call the MPI communication will hang because all the process will be waiting for a process that refused to send the data. This is the wrong code
if(pdata>5000){ // wrong way to try to selective pass data larger than 5000
mpi::all_reduce(world, pdata, std::plus<double>());
} // WRONG
What is worst: the program may or may not hang depending on the actual numeric value of the data. So, third lesson of parallel programming: it is too easy to create a deadlock. In particular all blocking-communication have to be executed in non-conditional branches of the program. (Another good reason to program without if's.)
reduce
Another variation of the example can be introduced if we realize that sometimes the result is only needed in one process. In this case the reduction result is sent to one process only, which is specified as the last argument in mpi::reduce:
double global_pdata=-1; mpi::reduce(world, pdata, global_pdata, std::plus<double>(),0); // process 0 receives the data cout<<"Process "<<world.rank()<<" *thinks* that global power of the original data is "<< global_pdata << endl;
In this example, only process 0 will have a meaningful value of the result; other processes will have the global_pdata variable unchanged. The process 0 will output the correct value, while the rest will output -1. Essentially, all process but one will have an uninitialized, invalid, or at best flagged to some magical value.
This selective data collection sounds reasonable idea but really it can cause trouble because we have to keep track of the process that has the right value by some conditional. There are two basic possibilities, either we "promise" to ignore the value in processes with the wrong value.
double global_pdata=-1;
mpi::reduce(world, pdata, global_pdata, std::plus<double>(),0); // process 0 receives (and send to itself) the data
if(world.rank()==0){ //always check for process number *before* using the variable.
cout<<"Process 0 knows that the global power of the origin data is "<<global_pdata<<", others do not really know"<<endl;
}
or we put the result variable inside the scope of an 'if' statement:
if(world.rank()==0){
double global_pdata=-1;
mpi::reduce(world, pdata, global_pdata, std::plus<double>(),0); // process 0 sends and receives the data
... // do other operations with result
cout<<"Process 0 knows that the global power of the origin data is "<<global_pdata<<", others do not know at all"<<endl;
}else{
mpi::reduce(world, pdata, std::plus<double>(),0); // others process only send the data
}
Note that the number of arguments is different in the two different calls to reduce. The first branch will receive in the result of all the processes (including itself) while the second will only send values.
In this way, only process 0 can ever access the value. Trying to use the the variable outside the first branch of the conditional will not even compile. In this case we do not have to promise about the usage of the variable but we have to promise that we will call similar versions of reduce in both branches of the conditional. If we forget the else branch (second branch), the program will hang in the same way as described above.
Next subsection will deal with this problem of 'dangling' variables. You can skip the following subsection in a first read, it will not improve your understanding of MPI. Just be aware the problem for the moment.
The issue of asymmetrical variables
A value that makes sense in only one process seems to be a natural situation, after all, each process is supposed to do different part of the total work. The problem is that, within MPI, all processes must share the same code and therefore will be forced to carry variables that may not make sense to them. This situation creates tension.
Note that in the first place, this issue arises because by collecting data in only one process we are imposing some asymmetry; it is nevertheless necessary sometimes. The situation is annoying for those used to serial programming, but we have to understand that this kind of problems are at the core of parallel programming difficulties. The first step to try to avoid the problem is to think whether it may be really useful --or at least not harming-- for all the processes to know the result of reduce/all_reduce.)
It seems that the second solution above (if/else block) is more elegant because there are no dangling undefined variables around, but is not free of problems. That is because if we need to compute something non-trivial with the result, the conditional statement can become several lines long and it can be difficult to handle cases in which more than one (but not all) process receive the value. We can put all these lines in a separate function, but that function will be forced to be serial code (because it is executed only for one process).
There is another possibility to consider that I discovered, that is the usage of Boost.Optional. (This solution is not suggested anywhere else to my knowledge.) In this case, the valid state of the variable is handled by the variable itself.
boost::optional<double> global_pdata(world.rank()==0, double()); //non empty at node 0 only.
if(global_pdata) mpi::reduce(world, pdata, *global_pdata, std::plus<double>(),0);
else mpi::reduce(world, pdata, std::plus<double>(),0);
...
if(global_pdata){ // only "use" if global_pdata is it is initialized
cout<<"Process "<<world.rank()
<<" knows that the global power of the origin data is "<<*global_pdata<<", others will die if they try to 'know'"<<endl;
}
// cout<<*global_pdata<<endl; // unchecked usage will abort the program in other processes
For those used to C, they will recognize that 'global_pdata' looks like a nullable pointer (whose pointee value is valid when !=0 by convention). The difference is that boost::optional behaves better than a pointer for this application, for example, boost::optional can be assigned values, or uninitialized with no memory allocation necessary, no dangling pointers can be created in this way. If we respect this check-then-use syntax, the variable becomes transparent to all the process where it is undefined.
Note that global_pdata is not associated with any particular process number, we can forget in which process the variable was valid, the validity test is done in the variable itself. We have to remember to test whether the global_pdata is locally valid. What if we do not check and use *global_pdata without test? Unless we compile in production mode (NDEBUG=0) the calling process will terminate after an 'uninitialized' fail assertion is called from Boost.Optional. This is still better than other options, like having a dangling pointer (nullable pointer solution), using an uninitialized or invalid value (first solution), or a deadlock (second solution).
So far I mentioned three (or four) possible solutions, none of them can save us at compilation time, i.e. before running our program. So, it is pretty much like saying 'choose your poison'. The situation will not improve until someone smart writes a compiler (and a language) capable of *really* interpreting parallel code and figuring out these issues before actually running the program. (Notice that MPI is really like low level programming.)
If you think that my suggested solution (using Boost.Optional) is in the right direction, there are at least two ways to improve it further. First, we can encapsulate the two calls to 'reduce' in a single function (a template function will be even more general):
boost::optional<double> optional_reduce(boost::mpi::communicator const& world, double const& in_value, int root){
boost::optional<double> ret(world.rank()==0, double());
if(ret) mpi::reduce(world, in_value, *ret, std::plus<double>(),0);
else mpi::reduce(world, in_value, std::plus<double>(),0);
return ret;
}
and then use it as
boost::optional<double> global_pdata = optional_reduce(world, pdata, 0); ... if(global_pdata) cout<<"global_pdata = "<<*global_pdata<<endl;
The second way to improve it is to wrap/extend optional<T> inside a class (that I will call likely<T>) that throws an exception (instead of failing after the assert) when accessing uninitialized variable, this gives the opportunity to recover from failed executions. For example:
try{
likely<double> global_pdata = likely_reduce(world, pdata,0);
... // use global_pdata
}catch(...){
clog<<"sorry I am a process that don't know about global_pdata"
}
The idea is that the invalid part of the try-block becomes transparent to processes. (The valid part should be exception safe though.) It would be a bad idea to use this syntax as a substitution for a conditional or for printing the result.