Interoperability

The fundamental goal of the library is that the arrays and iterators can be used out-of-the-box with C++ Standard Library, in particular the Standard Template Library (STL) and the algorithms it provides. The most dramatic example of this is that std::sort works with array as it is shown in a previous section.

Along with STL itself, the library tries to interact with other existing quality C++ libraries listed below.

The library works out of the box with Eric Niebler’s Range-v3 library, a precursor to the standard Ranges library (see above). The library helps to remove explicit iterators (e.g. begin, end) from the code when possible.

Every Multi array object can be regarded as range. Every subarray references (and array values) are interpreted as range views.

For example for a 2D array d2D, d2D itself is interpreted as a range of rows. Each row, in turn, is interpreted as a range of elements. In this way, d2D.transposed() is interpreted as a range of columns (of the original array), and each column a range of elements (arranged vertically in the original array).

In this other example, a 2D Multi array (or subarray) is modified such that each element of a column is subtracted the mean value of such column.

The ability to serialize arrays is essential for storing data in a persistent medium (files on disk) and communicating values via streams or networks (e.g., MPI). Unfortunately, the C++ language does not provide facilities for serialization, and the standard library doesn’t either.

However, there are a few libraries that offer a certain common protocol for serialization, such as Boost.Serialization and Cereal. The Multi library is compatible with both (and doesn’t depend on any of them). The user can choose one or the other, or none, if serialization is not needed. The generic protocol is such that variables are (de)serialized using the (>>)<< operator with the archive; operator & can be used to have a single code for both. Serialization can be binary (efficient) or text-based (human-readable).

Here, it is a small implementation of save and load functions for an array to JSON format with the Cereal library. The example can be easily adapted to other formats or libraries. (An alternative for XML with Boost.Serialization is commented on the right.)

(online)

These templated functions work for any dimension and element type (as long as the element type is serializable in itself; all basic types are serializable by default). However, note that the user must ensure that data is serialized and deserialized into the same type; the underlying serialization libraries only do minimal consistency checks for efficiency reasons and don’t try to second-guess file formats or contained types. Serialization is a relatively low-level feature for which efficiency and economy of bytes are a priority. Cryptic errors and crashes can occur if serialization libraries, file formats, or C++ types are mixed between writes and reads. Some formats are human-readable but still not particularly pretty for showing as output (see the section on Formatting on how to print to the screen).

References to subarrays (views) can also be serialized; however, size information is not saved in such cases. The reasoning is that references to subarrays cannot be resized in their number of elements if there is a size mismatch during deserialization. Therefore, array views should be deserialized as other array views with matching sizes.

The output JSON file created by Cereal in the previous example looks like this.

(The Cereal XML and Boost XML output would have a similar structure.)

Large datasets tend to be serialized slowly for archives with heavy formatting. Here it is a comparison of speeds when (de)serializing a 134 MB 4-dimensional array of with random double.

Using Interprocess allows for shared memory and for persistent mapped memory.

The library works out-of-the-box in combination with the Thrust library.

(live)

which uses the default Thrust device backend (i.e. CUDA when compiling with nvcc, HIP/ROCm when compiling with a HIP/ROCm compiler, or OpenMP or TBB in other cases). Universal memory (accessible from normal CPU code) can be used with thrust::universal_allocator (from <thrust/universal_allocator.h>) instead.

More specific allocators can be used ensure CUDA backends, for example CUDA managed memory:

In the same way, to ensure HIP backends please replace the cuda namespace by the hip namespace, and in the directory name <thrust/system/hip/memory.h>. <thrust/system/hip/memory.h> is provided by rocThrust in the ROCm distribution (in /opt/rocm/include/thrust/system/hip/, and not by the NVIDIA distribution.)

Multi doesn’t have a dependency on Thrust (or vice versa); they just work well together, both in terms of semantics and efficiency. Certain "patches" (to improve Thrust behavior) can be applied to Thrust to gain extra efficiency and achieve near native speed by adding the #include<multi/adaptors/thrust.hpp>.

Multi can be used on existing memory in a non-invasive way via (non-owning) reference arrays:

Finally, the element type of the device array has to be device-friendly to work correctly; this includes all build in types, and classes with basic device operations, such as construction, destruction, and assigment. They notably do not include std::complex<T>, in which can be replaced by the device-friendly thrust::complex<T> can be used as replacement.

CUDA is a dialect of C++ that allows writing pieces of code for GPU execution, known as "CUDA kernels". CUDA code is generally "low level" (less abstracted), but it can be used in combination with CUDA Thrust or the CUDA runtime library, specially to implement custom algorithms. Although code inside kernels has certain restrictions, most Multi features can be used. (Most functions in Multi, except those involving memory allocations, are marked device to allow this.)

Calling kernels involves a special syntax (<<< …​ >>>), and they cannot take arguments by reference (or by values that are not trivial). Since arrays are usually passed by reference (e.g. multi::array<double, 2>& or Array&&), a different idiom needs to be used. (Large arrays are not passed by value to avoid copies, but even if a copy would be fine, kernel arguments cannot allocate memory themselves.) Iterators (e.g. .begin()/.end()) and "cursors" (e.g. .home()) are "trivial to copy" and can be passed by value and represent a "proxy" to an array, including allowing the normal index syntax and other transformations.

Cursors are a generalization of iterators for multiple dimensions. They are cheaply copied (like iterators) and they allow indexing. Also, they have no associated .size() or .extensions(), but this is generally fine for kernels. (Since cursors have minimal information for indexing, they can save stack/register space in individual kernels.)

Here it is an example implementation for matrix multiplication, in combination with Thrust and Multi,

(live)

Expressions such as A.begin() (iterators) can also be passed to kernels, but they could unnecessarely occupy more kernel "stack space" when size information is not needed (e.g. A.begin()→size()).

The SYCL library promises the unify CPU, GPU and FPGA code. At the moment, the array containers can use the Unified Shared Memory (USM) allocator, but no other tests have been investigated.

(live)

Algorithms are expected to work with oneAPI execution policies as well (not tested)

The library doesn’t have a "pretty" printing facility to display arrays. Although it is not ideal, arrays can be printed and formated by looping over elements and dimension, as shown in other examples (using standard streams).

Fortunately, the library automatically works with the external library {fmt}, both for arrays and subarrays. The fmt library is not a dependency of the Multi library; they simply work well together using the "ranges" part of the formatting library. fmt allows a high degree of configurability.

This example prints a 2-dimensional subblock of a larger array:

yielding the "flat" output A2D subblock = [[3, 4], [6, 7]]. (A similar effect can be achieved with experimental C23 `std::print` in libc.)

For 2 or more dimensions the output can be conveniently structured in different lines using the fmt::join facility:

with the output:

In a similar way, the size of the array can be printed simply by passing the sizes output; fmt::print("{}", A2(…​).sizes() ), which will print (2, 2).

((live)).

When saving arrays to files, please consider using serialization (see section) instead of formatting facilities.

The library works out-of-the-box via de automatic cppyy bindings, providing a syntax similar to that of C. Additionally, the memory layouts are compatible with Numpy which makes possible to interoperate C and python code without memory copies.

Here it is a complete Python session:

Take into account that assignment in Python differs semantically from that of C++. This can be particularly confusing when assigning elements from a list.

The last line reassign the variable into a Python list, changing its type. To reassign the elements without changing the type, use .assign instead:

Pointers from arrays can be extracted normally and reused in other libraries, to access the elements through other Python libraries; for example, an array that uses contiguous memory can be converted into a Numpy array, without making copies

The other way around,

Note that the elements fpr these Python objects (a2d, npa2d and a2d_ref) all refer to the same memory; so modifying one, also nodifies the rest. Special care must be taken if the original element is resized, or assigned to a different Python object, in which case the memory will not be valid. For this reason alone, it would be an error to name the last variable as the original (a2d) instead of the new name (a2d_ref).

Multidimensional array data structures exist in all languages, whether implicitly defined by its strides structure or explicitly at the language level. Functions written in C tend to receive arrays by pointer arguments (e.g., to the "first" element) and memory layout (sizes and strides).

A C-function taking a 2D array with a concrete type might look like this in the general case:

such a function can be called from C++ on Multi array (arr), by extracting the size and layout information,

or

Although the recipe can be applied straightforwardly, different libraries make various assumptions about memory layouts (e.g., 2D arrays assume that the second stride is 1), and some might take stride information in a different way (e.g., FFTW doesn’t use strides but stride products). Furthermore, some arguments may need to be permuted if the function expects arrays in column-major (Fortran) ordering.

For these reasons, the library is accompanied by a series of adaptor libraries to popular C-based libraries, which can be found in the include/multi/adaptors/ subdirectory:

Interface for BLAS-like linear algebra libraries, such as openblas, Apple’s Accelerate, MKL and hipBLAS/cuBLAS (GPUs). Simply #include "multi/adaptors/blas.hpp" (and link your program with -lblas for example).

Interface for Lapack linear solver libraries. Simply #include "multi/adaptors/lapack.hpp" (and link your program with -llapack for example).

Interface for FFTW libraries, including FFTW 3, MKL, cuFFT/hipFFT (for GPU). Simply #include "multi/adaptors/fftw.hpp" (and link your program with -lfftw3 for example).

Use arrays (and subarrays) as messages for distributed interprocess communication (GPU and CPU) that can be passed to MPI functions through datatypes. Simply #include "multi/adaptors/mpi.hpp".

Popular in HPC environments, can display arrays in human-readable form (for simple types, like double or std::complex). Simply #include "multi/adaptors/totalview.hpp" and link to the TotalView libraries, compile and run the code with the TotalView debugger.