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 languages 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 (Fortran-convention) or row-major ordering (C-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.)
Basic Usage (multi_array)
MultiArray is a header-only library -- no special linking is necessary. You just need to include the header file:
#include<boost/multi_array.hpp>
The usage of the multi_array class is pretty intuitive. First we declare and construct an variable of type multi_array, specifying its dimensionality (2) and its dimensions (100x100).
boost::multi_array<double, 2> A(boost::extents[10][10]);
The number of arguments of extents has to agree with the dimensionality. Then you can use the matrix as usual,
for(int i=0; i!=A.shape()[0]; ++i)
for(int j=0; j!=A.shape()[1]; ++j)
A[i][j]=i+j;
For an N-dimensional array we would need N similar loops. For any dimensionality, the elements can be iterated in its underlying memory representation:
for(double* it=A.origin(); it!=A.origin()+A.num_elements(); ++it) (*it)*=2;
In this example, the array is multiplied by two, regardless of its dimensionality. Note that num_elements()==shape()[0]*shape()[1]...shape()[N].
The dimensionality of a *particular array* can be known by using A.num_dimensions(). The dimensionality of the *type* can be known by using the static constant "::dimensionality".
typedef boost::multi_array<double, 4> marray; marray A; assert( A.num_dimensions() == marray::dimensionality );
The official documentation and this paper have additional examples and some explanation of the library internals.
Using MultiArray on existing code together with C-arrays (multi_array_ref)
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_ref); ... 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).
Fortran ordering and a simple Fortran binding (fortran_storage_order)
Functions in libraries like Lapack expect arrays in column-major order format, this option 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. (By default, if not specified, boost::c_storage_order is used.)
boost::multi_array<double,2> A(boost::extents[10][10], boost::fortran_storage_order); ... // assign A value with Fortran index ordering
In order to imitate Fortran convention, there is even an option to declare the multi_array to be 1-based (instead of 0-based) (see below).
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 'wrapper_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.
Custom index ranges (reindex)
In previous examples we used some obscure entity called boost::extents, it is obviously used to pass the valid range of indices 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)].)
This is a good point to clarify the difference between A.data() and A.origin(). As we saw before origin() returns a pointer to the first element of the array, what ever this is. On the other hand, A.data() returns a pointer to A[0][0]..[0], regardless of whether that element is the first or not. The element A[0][0]..[0] may not even be part of the array, for example in 1-based arrays. In general we have to deal origin() and rarely with data(), but it is very easy to confuse one with the other.
It is also good to know that the memory allocated by an array is always located between A.origin() and A.origin()+A.num_elements() (non-inclusive).
Regular subarrays
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) row 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. The following examples will create arrays from the original 2D array A:
typedef std::complex<double> complex; typedef boost::multi_array<complex,2> array2d; using boost::extents; array2d A(extents[10][10]); // 2D, 100 allocated elements
Subarrays (sub_array)
typedef array2d::subarray<1>::type subarray1d; subarray1d A_row5 = A[5]; // 1D, 10 elements
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_row3_5to10 = A[3][indices[range(5,10)]]; // or A[indices[3][range(5,10)] view1d A_row3_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 arrays defined here make any copy of elements, remember that sub_arrays and array_view are just references, another way to view the original array.
Views are more general than subarrays. For example A[1] is an object of type subarray, which can be implicitly converted into a view. However A[1][indices[range(0,5)]] is an array_view object which can not be converted into a subarray.
A simple FFTW wrapper
Here I 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. Sometimes need more flexibility, like the one given by C++ classes. Naturally, this will need a deep understanding of the language and of the wrapped libraries. But once written they will make easier to program.
For example, in the case of FFTW, although written in C, it is obvious that 'fftw_plan' structure behaves almost as a C++ 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:
namespace fftw3{
using namespace boost;
class plan : private counted<plan>{
plan(fftw_plan const& p) : p_(p) {}
public:
template<class ConstMultiArray, class MultiArray>
plan(ConstMultiArray const& in, MultiArray const& out, int direction, unsigned flags=FFTW_ESTIMATE){
assert(boost::shape(in)==boost::shape(out));
if(in.num_dimensions()>1){
for(int i=in.num_dimensions()-1, dense_in=in.strides()[i], dense_out=out.strides()[i]; i!=-1; dense_in*=in.shape()[i], dense_out[i]*=out.shape()[i], --i){
assert(dense_in==in.strides()[i] and dense_out==out.strides()[i]); // check for dense representation of multidimensional
}
}
p_ = fftw_plan_many_dft(in.num_dimensions(), (const int *)in.shape(), 1 /*int howmany*/,
(double(*)[2])in.origin() , 0 /*const int *inembed, 0=in.shape()*/,
in.strides()[in.num_dimensions()-1] /*int istride*/, 0 /*int idist ignored for howmany=1*/,
(double(*)[2])out.origin(), 0 /*const int *onembed, 0=out.shape()*/,
out.strides()[out.num_dimensions()-1] /*int ostride*/, 0 /*int odist ignored for howmany=1*/,
direction, flags);
}
void operator()() const{return fftw_execute(p_);}
virtual ~plan(){fftw_destroy_plan(p_);}
private:
fftw_plan p_;
public:
static void cleanup(){
if(get_count()<=0) throw std::logic_error("trying to call fftw3::cleanup() with "+boost::lexical_cast<std::string>(get_count())+" active plans, should be zero");
fftw_cleanup();
}
};
void execute(plan const& p) {p.operator()();}
}
Then 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 execute(p); // 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 destruction needed (no risk to execute a destroyed plan).
The wrapper can be used with the condition that the input and output arrays are stored densely in memory. (This is not the case for some array_views.) This is a limitation of the library (which was designed to take advantage of contiguous memory) and this limitation must be reflected in the wrapper. Besides that, the wrapper can deal with arrays, sub_arrays and array_views (with strides==1 if multidimensional).
In-place transform can be achieved simply by using the same argument in the first and second placeholder. The wrapper will also check for size compatibility.
Using the original arrays, there are many possible ways to use the wrapper:
fftw3::plan p(A, B, FFTW_FORWARD); // 2D FFT 100x100 fftw3::plan p(A, A, FFTW_FORWARD); // 2D inplace-FFT fftw3::plan p(A[21], B[32], FFTW_FORWARD); // 1D FFT 100 elements, fftw3::plan p(A[10][indices[range(5,10)], B[20][indices[range(5,10)]], FFTW_FORWARD); // 1D FFT on 5 elements
If not specified, the option FFTW_ESTIMATE is used. If using this default option the point at which the plan is created is unimportant. Instead, when using the options FFTW_MEASURE or FFTW_PATIENT, the location of plan creation matters. This is so because in these cases the arrays can be overwritten when FFTW looks for the optimal algorithm to use. That means that the plans should be created before the matrix is assigned values and after the matrix has the desired size and shape.
multi_array<comple,2> A(extents[100][100]); fftw3::plan p(A,A, FFTW_FORWARD, FFTW_MEASURE); //overwrites values ... // assign values; execute(p); // performs the transform
Once the plans are created over existing arrays these arrays must be not reallocated! That means not resize operations please. Assignment (global or partial) is OK because in this case MultiArray guarranties no internal reallocation is done. (Global assignment is only allowed for arrays with same sizes). The is reflects a limitation of the FFTW in which case the array pointer doesn't have to change after the plan is created.
'Many' Transforms
The wrapper can also be asked to perform multiple individual transform, by assuming that the first index of the array is used to index smaller arrays and using the following syntax:
fftw3::plan p=fftw3::plan::many(A, B, FFTW_FORWARD); // 100 1D-FFT of 100 elements each
note that this is not equivalent to the first example.
Again, the hard work is done by the library, which uses certain optimizations to perform repeated transforms of the same size.
Special Memory Allocation
FFTW recommends to use its own version of malloc to allocate arrays. This is supposedly to improve speed taking advantage of some special allocation that is exploited by the numerical routines. For example they recommend to allocate the arrays to be used in the following way:
typedef std::complex<double> complex; complex* in=(complex*)fftw_malloc(sizeof(complex)*N); //instead of malloc(sizeof(complex)*N); fftw_free(in); //instead of free(in)
Up to this point we tried not to deal with memory allocation at all.
So, the question is: How can we take advantage of allocation of memory in a context of C++? (where manual allocation can --and should-- be avoided in most cases). The answer is 'custom allocators'. For STL containers and for multi_arrays, the allocator is passed as an extra template parameter.
In this example, I pass an special allocator that uses fftw_malloc internally, for the construction of an STL vector and also for a Boost.MultiArray multi_array.
std::vector<complex, fftw3::allocator<complex> > v(100); boost::multi_array<complex, 3, fftw3::allocator<complex> > A(extents[100][100][100]);
The memory is freed (automatically) by fftw_free. After this, the object can be used like any other standard object. If the allocator is not specified, std::allocator<T> is used.
Declaring the objects in this way is optional in the same way as the usage of fftw_malloc instead of malloc is optional.
Reshape (reshape)
Although MultiArray objects can be resized, the underlying memory layout is not specially good for resizing operations. In general, for resizing, the memory is reallocated *all* the elements copied. This is even true for operations reducing the size of the matrix (since elements have to be contiguous). So, resizing operations should be avoided as much as possible and they will not be covered here (otherwise see here).
Multidimensional arrays are different from one dimensional arrays because they can be reshaped. In the library, reshape operations are different from resize in the sense that they preserve the total number of elements. The contained elements are rearranged logically but no copies or reallocations are performed. The reshape operation is done as follows:
boost::multi_array<complex, 3> a(extents[5][4][3]); //assert(a.shape()[0]==5 and a.shape()[1]==4 and a.shape()[2]==3);
{
boost::array<int,3> new_shape={{5,6,2}}; //assert(5*6*2==5*4*3);
a.reshape(new_shape); //assert(a.shape()[0]==5 and a.shape()[1]==6 and a.shape()[2]==2);
}
In this example, the reshape operation is allowed because the number of elements is preserved since 5*6*2==5*4*3. If the number of elements is not preserved the operation will fail during execution.
The following three syntaxes are not allowed due to very strange technical limitations of C++:
a.reshape({{5,6,2}}); // does not compile
a.reshape({5,6,2}); // neither
a.reshape(boost::array<int, 3>({{5,6,2}})); // neither
What eventually works is the following syntax:
a.reshape((boost::array<int,3>){{5,6,2}}); // watch for parenthesis
or, depending on the context, it might be better to use:
using boost::array;
typedef array<int,3> array3;
a.reshape((array3){{5,6,2}});
Note 1: boost::array is not the same as boost::multi_array. boost::arrays objects are part of a much simpler library called Boost.Array. In this other library all the arrays are one-dimensional and the number after the element type indicates the (constant) number of elements and not the number of dimensions.
So far, the only constrain for reshaping is that it preserves the number of elements. The other requirement that is not very obvious is that the reshape operation must preserve the number of dimensions (in this case 3). This is not a big limitation since we can choose one of the dimensions to be 1 in size
a.reshape((array<int,3>){{5,1,12}}); // or {{5,12,1}} depending on your application
Note 2: Do not try to do something like a.reshape((array<int,2>){{5,12}}) for a 3D array to change dimensionality, it compiles but gives unexpected results. In fact, there is no way to have 'a' change its dimensionality during the program execution, since it is of type multi_array<double,3>. Another thing that is very confusing is that reshape does not work with extents; do not A.reshape(extents[5][4][3]). Again it compiles and gives the wrong result.
If you really need an object with smaller dimensionality with the same elements you can do
boost::multi_array<complex, 3> a(extents[5][4][3]); boost::multi_array_ref<complex, 2> a_ref(a.data(), extents[5][12]);
and use 'a_ref' instead of 'a'. (But in this case nobody will check that the sizes are compatible unless you add something like assert(a.num_elements()==a_ref.num_elements()); .)
Note 3: Of course you can always copy element by element from one type of array to another with different shape or different dimensionality. But this section is exactly all about avoiding copies.
Example of reshape with FFTW
The FFTW wrapper presented earlier honors the dimensionality of the arrays passed. I mean by this that the dimensionality of the transform follows the dimensionality of the passed array instead of for example interpreting the input as a one dimensional array which would of course give a completely different result. In this section I will illustrate how to force FFTW to interpret the data differently.
For simplicity let us assume that we have a two dimensional array, and that we want to do our transforms in-place on complex data:
typedef std::complex<double> complex; using boost::multi_array; using boost::extents; multi_array<complex, 2> A(extents[100][100]); ... // fill A with data
The simplest case is that we want a 2D FFT on the 2D data structure, I present it here as a review:
fftw3::plan pf_A2A_2D(A, A, FFTW_FORWARD); // implicitly 2D transform execute(pf_A2A_2D);
Now suppose that we need to do a different thing that is to transform each row of the the array independently, the approach in this case will be different, for example we could iterate over the rows and perform a 1D transform (this is allowed because each A row is contiguous data in memory):
for(int i=0; i!=A.shape()[0]; i++){
fftw3::plan pf_Arow2Arow_1D(A[i],A[i], FFTW_FORWARD); // implicitly 1D transform
execute(pf_Arow2Arow_1D);
}
Below there is a note against using FFTW plans inside loops. A different approach, *to do the same thing* would be to use the 'many' option of the library and the wrapper. This can be faster because the library can take advantage of how the overall data is distributed.
fftw3::plan pf_Arow2Arow_many_1D=fftw3::plan::many(A, A, FFTW_FORWARD); // many 1D transforms execute(pf_Arow2Arow_many_1D);
The last option is faster, there is no explicit loop, and only one plan is created, this is more elegant because the plan can be created just after the array is constructed, this is handy if want to create plans with the FFTW_MEASURE option. (But keep in mind that this is possible only because all the 1D transforms are performed inside a single multi_array.)
Note: There is also an insidious issue with creating many FFTW plans that I discovered accidentally. The library keeps track of all the created plans (even the 'destroyed' ones). Each created/destroyed plan leaves a small memory footprint. FFTW authors claim that this is for performance optimization purposes. But when creating, using and destroying plans repeatedly, let's say few million times, the used memory will grow up and the application crash. Authors of the library say this is not a memory leak, and in part they are right, just because there a function called fftw_cleanup() [called fftw3::plan::cleanup() in the wrapper] that makes sure all that memory is reclaimed. The problem is that the function must be called at some point where it is known that there are no plans active (the wrapper will checks for this condition). It is your responsibility to call this function and since there are only few points in the program at which we know that there are no plans active, where to call cleanup(), this *effectively* behaves like a memory leak. The way to avoid this altogether is to keep the number of created plans controlled, and specially avoid creating plans inside loops. The 'many' plan can save lots of plans to be created.
Now suppose that we have our data in a 2D multi_array but that we want to interpret the data differently and perform a long 1D transform over a contiguous sequence of rows in the array. This can be motivated by special boundary conditions in the problem but also illustrates the utility of reshape.
There are two options here, one using reshape and the other using multi_array_ref (see above). By using reshape we *temporarely* rearrange the array data (logically) and perform the transform:
using boost::array;
A.reshape((array<int,2>){{1,10000}}); //remember that original was 100x100
fftw3::plan pf_Aunwrappedrows_1D(A[0],A[0], FFTW_FORWARD); // implicitly 1D
A.reshape((array<int,2>){{100,100}});
execute(pf_Aunwrappedrows_1D);
First, we reshape the original matrix to be effectively a one dimensional array. Strictly speaking A is 2D, but it has only one row A[0] and we create a 1D FFT plan for that single row that contains all the data. After the plan is created we can reshape the array again, since FFTW will remember the shape of the array. Another variant is to simply declare a plan:
fftw3::plan pf_Aunwrappedrows_2D(A, A, FFTW_FORWARD);
Although in the case the arrays and the plans are 2D, the result should be the same as a 1D transform since there is only one row.
The second option is to create a multi_array reference with a different 'view' and work on that matrix:
using boost::multi_array_ref; multi_array_ref<complex, 1> Aunwrapped_ref(A.origin(), extents[10000]); assert(Aunwrapped_ref.num_elements()==A.num_elements()); fftw3::plan pf_Aunwrapped_1D(Aunwrapped_ref, Aunwrapped_ref, FFTW_FORWARD); // inplicitly 1D execute(pf_Aunwrapped_1D);
The FFTW will modify the reference matrix which is only a different way of changing the original array but the element access is logically different.
Boost.MPI Library
I will assume that it is the case that you 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 Install 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). It can be regarded as the parallel version of this serial example. 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 */
fftw_free(data);
// 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, because the splitting depends on the type of transform to be performed (e.g. in-place, out-of-place).
- 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). Even in this case, the allocated memory in that process can be different from zero! (in general, it is 1*sizeof(complex)).
- 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 results in the number 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* associative, 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).
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.)
Quick note on function objects (and Boost.Function)
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;
In fact it is very easy to create your own function object, they look almost like functions:
struct concat{ //concatenation
std::string operator()(std::string const& s1, std::string const& s2){
return s1 + s2;
}
};
Passing such function object to reduce will make all_reduce concatenate the strings produced by different processes. Note that this operation is not commutative and will depend on some standard convention used by MPI to collect the data (fortunately it is well defined).
The only difference between our function object an std::plus<std::string>, is that the later is already included in some standard library.
They differ from functions in two aspects, 1) you can parametrize them, 2) you can make copies of them, 3) you can compare them:
struct concat_with_sep{ // concatenation with separation
std::string separator_;
concat_with_sep(std::string const& separator) : separator_(separator){}
std::string operator()(std::string const& s1, std::string const& s2){
return s1 + separator_ + s2;
}
};
The beauty of this is that, although the separator string is arbitrary, 'concat_with_sep', viewed as function, still takes only two parameters.
Free functions (the normal C-functions) can be converted to function objects in many ways. The standard way is to use [1] library. Function objects can be combined with one another with the [2] and the (more advanced) [3].
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 dangling variables (and Boost.Optional)
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). Note that I refer this as "dangling variable", in the same sense we refer to "dangling pointers", both are equally dangerous. The following proposed solution uses the library Boost.Optional and it can be also used as an introduction to this other Boost library.
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.