1.1 Rank polymorphism

Rank polymorphism is the ability to treat an array of rank \(r\) as an array of lower rank where the elements are themselves arrays.

For example, think of a matrix A, a 2-array with shape (n₀, n₁) where the elements A(i₀, i₁) are numbers. If we consider the subarrays A(0, ...), A(1, ...), ..., A(n₀-1, ...) as individual elements, then we have a new view of A as a 1-array of length n₀ with those rows as elements. We say that the rows A(i₀)≡A(i₀, ...) are the 1-cells of A, and the numbers A(i₀, i₁) are 0-cells of A. For an array of arbitrary rank \(r\) the (\(r\)-1)-cells of A are called its items. The prefix of the shape (n₀, n₁, ... nₙ₋₁₋ₖ) that is not taken up by the k-cell is called the k-frame.

An obvious way to store an array in linearly addressed memory is to place its items one after another. So we would store a 3-array as

A: [A(0), A(1), ...]

and the items of A(i₀), etc. are in turn stored in the same way, so

A: [A(0): [A(0, 0), A(0, 1) ...], ...]

and the same for the items of A(i₀, i₁), etc.

A: [[A(0, 0): [A(0, 0, 0), A(0, 0, 1) ...], A(0, 1): [A(0, 1, 0), A(0, 1, 1) ...]], ...]

This way to lay out an array in memory is called row-major order or C-order, since it’s the default order for built-in arrays in C (see Other array libraries). A row-major array A with shape (n₀, n₁, ... nᵣ₋₁) can be looked up like this:

A(i₀, i₁, ...) = (storage-of-A) [(((i₀n₁ + i₁)n₂ + i₂)n₃ + ...)+iᵣ₋₁] = (storage-of-A) [o + s₀i₀ + s₁i₁ + ...]

where the numbers (s₀, s₁, ...) are called the steps³. Note that the ‘linear’ or ‘raveled’ address [o + s₀i₀ + s₁i₁ + ...] is an affine function of (i₀, i₁, ...). If we represent an array as a tuple

A ≡ ((storage-of-A), o, (s₀, s₁, ...))

then any affine transformation of the indices can be achieved simply by modifying the numbers (o, (s₀, s₁, ...)), with no need to touch the storage. This includes common operations such as: transposing axes, reversing the order along an axis, most cases of slicing, and sometimes even reshaping or tiling the array.

A basic example is obtaining the i₀-th item of A:

A(i₀) ≡ ((storage-of-A), o+s₀i₀, (s₁, ...))

Note that we can iterate over these items by simply bumping the pointer o+s₀i₀. This means that iterating over (k>0)-cells doesn’t cost any more than iterating over 0-cells (see Cell iteration).

1.2 Drag along and beating

These two fundamental array optimizations are described in [Abr70] .

Drag-along is the process that delays evaluation of array operations. Expression templates can be seen as an implementation of drag-along. Drag-along isn’t an optimization in and of itself; it simply preserves the necessary information up to the point where the expression can be executed efficiently.

Beating is the implementation of certain array operations on the array view descriptor instead of on the array contents. For example, if A is a 1-array, one can implement reverse(A, 0) by negating the step and moving the offset to the other end of the array, without having to move any elements. More generally, beating applies to any function-of-indices (generator) that can take the place of an array in an array expression. For instance, an expression such as 1+iota(3, 0) can be beaten into iota(3, 1), and this can enable further optimizations.

1.3 Why C++

Of course the main reason is that (this being a personal project) I’m more familiar with C++ than with other languages to which the following might apply.

C++ supports the low level control that is necessary for interoperation with external libraries and languages, but still has the abstraction power to create the features we want even though the language has no native support for most of them.

The classic array languages, APL [FI73] and J [Ric08] , have array support baked in. The same is true for other languages with array facilities such as Fortran or Octave/Matlab. Array libraries for general purpose languages usually depend heavily on C extensions. In Numpy’s case [num17] this is both for reasons of flexibility (e.g. to obtain predictable memory layout and machine types) and of performance.

On the other extreme, an array library for C would be hampered by the limited means of abstraction in the language (no polymorphism, no metaprogramming, etc.) so the natural choice of C programmers is to resort to code generators, which eventually turn into new languages.

In C++, a library is enough.

1.4 Guidelines

ra:: attempts to be general, consistent, and transparent.

Generality is achieved by removing arbitrary restrictions and by adopting the rank extension mechanism of J. ra:: supports array operations with an arbitrary number of arguments. Any of the arguments in an array expression can be read from or written to. Arrays or array expressions can be of any rank. Slicing operations work for subscripts of any rank, as in APL. You can use your own types as array elements.

Consistency is achieved by having a clear set of concepts and having the realizations of those concepts adhere to the concept as closely as possible. ra:: offers a few different types of views and containers, but it should be possible to use them interchangeably whenever it makes sense. For example, it used to be the case that you couldn’t create a higher rank iterator on a ViewSmall, even though you could do it on a View; this was a bug.

Sometimes consistency requires a choice. For example, given array views A and B, A=B copies the contents of view B into view A. To change view A instead (to treat A as a pointer) would be the default meaning of A=B for C++ types, and result in better consistency with the rest of the language, but I have decided that having consistency between views and containers (which ‘are’ their contents in a sense that views aren’t) is more important.

Transparency is achieved by avoiding unnecessary abstraction. An array view consists of a pointer and a list of steps and I see no point in hiding it. Manipulating the steps directly is often useful. A container consists of storage and a view and that isn’t hidden either. Some of the types have an obscure implementation but I consider that a defect. Ideally you should be able to rewrite expressions on the fly, or plug in your own traversal methods or storage handling.

That isn’t to mean that you need to be in command of a lot of internal detail to be able to use the library. I hope to have provided a high level interface to most operations and a reasonably usable syntax. However, transparency is critical to achieve interoperation with external libraries and languages. When you need to, you’ll be able to guarantee that an array is stored in compact columns, or that the real parts are interleaved with the imaginary parts.

1.5 Other array libraries

Here I try to list the C++ array libraries that I know of, or libraries that I think deserve a mention for the way they deal with arrays. It is not an extensive review, since I have only used a few of these libraries myself. Please follow the links if you want to be properly informed.

Since the C++ standard library doesn’t offer a standard multidimensional array type, some libraries for specific tasks (linear algebra operations, finite elements, optimization) offer an accessory array library, which may be more or less general. Other libraries have generic array interfaces without needing to provide an array type. FFTW is a good example, maybe because it isn’t C++!

• Standard C++
• Blitz++
• Other C++ libraries
• Other languages

2.1 Using ra::

ra:: is a header only library with no dependencies other than the standard library, so you just need to place the ‘ra/’ folder somewhere in your include path and add #include "ra/ra.hh" at the top of your sources.

A compiler with C++20 support is required. For the current version this means at least gcc 11 with -std=c++20. Some C++23 features are available with -std=c++2b. Check the top README.md for more up-to-date information.

Here is a minimal program⁴:

#include "ra/ra.hh"
#include <iostream>

int main()
{
    ra::Big<char, 2> A({2, 5}, "helloworld");
    std::cout << ra::noshape << format_array(transpose<1, 0>(A), {.sep0="|"}) << std::endl;
}

-| h|w
   e|o
   l|r
   l|l
   d|d

The following headers are not included by default:

• "ra/dual.hh": A dual number type for simple uses of automatic differentiation.
• "ra/test.hh": Used by the test suite.

The header "ra/base.hh" can be used to configure Error handling. You don’t need to modify the header, but the configuration depends on including "ra/base.hh" before the rest of ra:: in order to override the default handler. All other headers are for internal use by ra::.

The following #defines affect the behavior of ra::.

• RA_DO_CHECK (default 1)

  If 1, check shape agreement (e.g. Big<int, 1> {2, 3} + Big<int, 1> {1, 2, 3}) and random array accesses (e.g. Small<int, 2> a = 0; int i = 10; a[i] = 0;). See Error handling.

• RA_DO_OPT_SMALLVECTOR (default 0)

  If 1, perform immediately certain operations on ra::Small objects, using small vector intrinsics. Currently this only works on gcc and doesn’t necessarily result in improved performance.

• RA_DO_FMA (default FP_FAST_FMA if defined, else 0)

  If 1, use fma in certain array reductions such as dot, gemm, etc.

2.2 Containers and views

ra:: offers two kinds of data objects. The first kind, the container, owns its data. Creating a container uses up memory and destroying it causes that memory to be freed.

There are three kinds of containers (ct: compile-time, rt: runtime): 1) ct size, 2) ct rank/rt shape, and 3) rt rank; rt rank implies rt shape. Some rt size arrays can be resized but rt rank arrays cannot normally have their rank changed. Instead, you create a new container or view with the rank you want.

For example:

{
    ra::Small<double, 2, 3> a(0.);     // a ct size 2x3 array
    ra::Big<double, 2> b({2, 3}, 0.);  // a rt size 2x3 array
    ra::Big<double> c({2, 3}, 0.);     // a rt rank 2x3 array
    // a, b, c destroyed at end of scope
}

Using the right kind of container can result in better performance. Ct shapes do not need to be stored in memory, which matters when you have many small arrays. Ct shape or ct rank arrays are also safer to use; sometimes ra:: will be able to detect errors in the shapes or ranks of array operands at compile time, if the appropriate types are used.

Container constructors come in two forms. The first form takes a single argument which is copied into the new container. This argument provides shape information if the container type requires it.⁵

using ra::Small, ra::Big;
Small<int, 2, 2> a = {{1, 2}, {3, 4}};  // explicit contents
Big<int, 2> a1 = {{1, 2}, {3, 4}};      // explicit contents
Small<int, 2, 2> a2 = {{1, 2}};         // error: bad shape
Small<int, 2, 2> b = 7;                 // 7 is copied into b
Small<int, 2, 2> c = a;                 // the contents of a are copied into c
Big<int> d = a;                         // d takes the shape of a and a is copied into d
Big<int> e = 0;                         // e is a 0-array with one element f()==0.

The second form takes two arguments, one giving the shape, the second the contents.

ra::Big<double, 2> a({2, 3}, 1.);     // a has shape [2 3], filled with 1.
ra::Big<double> b({2, 3}, ra::none);  // b has shape [2 3], default initialized
ra::Big<double> c({2, 3}, a);         // c has shape [2 3], a is copied into c

The last example may result in an error if the shape of a and (2, 3) don’t match. Here the shape of 1. [which is ()] matches (2, 3) by a mechanism of rank extension (see Rank extension). The special value ra::none can be used to request default initialization of the container’s elements.

The shape argument can have rank 0 only for rank 1 arrays.

ra::Big<int> c(3, 0);      // ok {0, 0, 0}, same as ra::Big<int> c({3}, 0)
ra::Big<int, 1> c(3, 0);   // ok {0, 0, 0}, same as ra::Big<int, 1> c({3}, 0)
ra::Big<int, 2> c({3}, 0); // error: bad length for shape
ra::Big<int, 2> c(3, 0);   // error: bad length for shape

When the content argument is a pointer or a 1D brace list, it’s handled especially, not for shape⁶, but only as the (row-major) ravel of the content. The pointer constructor is unsafe —use at your own risk!⁷

Small<int, 2, 2> aa = {1, 2, 3, 4}; // ravel of the content

ra::Big<double, 2> a({2, 3}, {1, 2, 3, 4, 5, 6}); // same as a = {{1, 2, 3}, {4, 5, 6}}

double bx[6] = {1, 2, 3, 4, 5, 6}
ra::Big<double, 2> b({3, 2}, bx); // {{1, 2}, {3, 4}, {5, 6}}

double cx[4] = {1, 2, 3, 4}
ra::Big<double, 2> c({3, 2}, cx); // *** WHO NOSE ***

using lens = mp::int_list<2, 3>;
using steps = mp::int_list<1, 2>;
ra::SmallArray<double, lens, steps> a {{1, 2, 3}, {4, 5, 6}}; // stored column-major: 1 4 2 5 3 6

These produce compile time errors:

Big<int, 2> b = {1, 2, 3, 4};           // error: shape cannot be deduced from ravel
Small<int, 2, 2> b = {1, 2, 3, 4 5};    // error: bad size
Small<int, 2, 2> b = {1, 2, 3};         // error: bad size

Sometimes the pointer constructor gets in the way (see scalar):

ra::Big<char const *, 1> A({3}, "hello"); // error: try to convert char to char const *
ra::Big<char const *, 1> A({3}, ra::scalar("hello")); // ok, "hello" is a single item
cout << ra::noshape << format_array(A, {.sep0="|"}) << endl;

-| hello|hello|hello

A view is similar to a container in that it points to actual data in memory. However, the view doesn’t own that data and destroying the view won’t affect it. For example:

ra::Big<double> c({2, 3}, 0.);     // a rt rank 2x3 array
{
    auto c1 = c(1);                // the second row of array c
    // c1 is destroyed here
}
cout << c(1, 1) << endl;           // ok

The data accessed through a view is the data of the ‘root’ container, so modifying the former will be reflected in the latter.

ra::Big<double> c({2, 3}, 0.);
auto c1 = c(1);
c1(2) = 9.;                        // c(1, 2) = 9.

Just as for containers, there are separate types of views depending on whether the shape is known at compile time, the rank is known at compile time but the shape is not, or neither the shape nor the rank are known at compile time. ra:: has functions to create the most common kinds of views:

ra::Big<double> c {{1, 2, 3}, {4, 5, 6}};
auto ct = transpose<1, 0>(c); // {{1, 4}, {2, 5}, {3, 6}}
auto cr = reverse(c, 0); // {{4, 5, 6}, {1, 2, 3}}

However, views can point to anywhere in memory and that memory doesn’t have to belong to an ra:: container. For example:

int raw[6] = {1, 2, 3, 4, 5, 6};
ra::View<int> v1({{2, 3}, {3, 1}}, raw); // view with shape [2, 3] steps [3, 1]
ra::View<int> v2({2, 3}, raw);           // same, default C (row-major) steps

Containers can be treated as views of the same kind (rt or ct) . If you declare a function

void f(ra::View<int, 3> & v);

you may pass it an object of type ra::Big<int, 3>.

2.3 Array operations

To apply an operation to each element of an array, use the function for_each. The array is traversed in an order that is decided by the library.

ra::Small<double, 2, 3> a = {{1, 2, 3}, {4, 5, 6}};
double s = 0.;
for_each([&s](auto && a) { s+=a; }, a);

⇒ s = 21.

To construct an array expression but stop short of traversing it, use the function map. The expression will be traversed when it’s assigned to a view, printed out, etc.

using T = ra::Small<double, 2, 2>;
T a = {{1, 2}, {3, 4}};
T b = {{10, 20}, {30, 40}};
T c = map([](auto && a, auto && b) { return a+b; }, a, b); // (1)

⇒ c = {{11, 22}, {33, 44}}

Expressions may take any number of arguments and be nested arbitrarily.

T d = 0;
for_each([](auto && a, auto && b, auto && d) { d = a+b; },
         a, b, d); // same as (1)
for_each([](auto && ab, auto && d) { d = ab; },
         map([](auto && a, auto && b) { return a+b; },
             a, b),
         d); // same as (1)

The operator of an expression may return a reference and you may assign to an expression in that case. ra:: will complain if the expression is somehow not assignable.

T d = 0;
map([](auto & d) -> decltype(auto) { return d; }, d) // just pass d along
  = map([](auto && a, auto && b) { return a+b; }, a, b); // same as (1)

ra:: defines many shortcuts for common array operations. You can of course just do:

T c = a+b; // same as (1)

2.4 Rank extension

Rank extension is the mechanism that allows R+S to be defined even when R, S may have different ranks. The idea is an interpolation of the following basic cases.

Suppose first that R and S have the same rank. We require that the shapes be the same. Then the shape of R+S will be the same as the shape of either R or S and the elements of R+S will be

(R+S)(i₀ i₁ ... i₍ᵣ₋₁₎) = R(i₀ i₁ ... i₍ᵣ₋₁₎) + S(i₀ i₁ ... i₍ᵣ₋₁₎)

where r is the rank of R.

Now suppose that S has rank 0. The shape of R+S is the same as the shape of R and the elements of R+S will be

(R+S)(i₀ i₁ ... i₍ᵣ₋₁₎) = R(i₀ i₁ ... i₍ᵣ₋₁₎) + S().

The two rules above are supported by all primitive array languages, e.g. Matlab [Mat] . But suppose that S has rank s, where 0<s<r. Looking at the expressions above, it seems natural to define R+S by

(R+S)(i₀ i₁ ... i₍ₛ₋₁₎ ... i₍ᵣ₋₁₎) = R(i₀ i₁ ... i₍ₛ₋₁₎ ... i₍ᵣ₋₁₎) + S(i₀ i₁ ... i₍ₛ₋₁₎).

That is, after we run out of indices in S, we simply repeat the elements. We have aligned the shapes so:

[n₀ n₁ ... n₍ₛ₋₁₎ ... n₍ᵣ₋₁₎]
[n₀ n₁ ... n₍ₛ₋₁₎]

This rank extension rule is used by the J language [J S] and is known as prefix agreement. The opposite rule of suffix agreement is used, for example, in Numpy [num17] ⁸.

As you can verify, the prefix agreement rule is distributive. Therefore it can be applied to nested expressions or to expressions with any number of arguments. It is applied systematically throughout ra::, even in assignments. For example,

ra::Small<int, 3> x {3, 5, 9};
ra::Small<int, 3, 2> a = x; // assign x(i) to each a(i, j)

⇒ a = {{3, 3}, {5, 5}, {9, 9}}

ra::Small<int, 3> x(0.);
ra::Small<int, 3, 2> a = {{1, 2}, {3, 4}, {5, 6}};
x += a; // sum the rows of a

⇒ x = {3, 7, 11}

ra::Big<double, 3> a({5, 3, 3}, ra::_0);
ra::Big<double, 1> b({5}, 0.);
b += transpose<0, 1, 1>(a); // b(i) = ∑ⱼ a(i, j, j)

⇒ b = {0, 3, 6, 9, 12}

An weakness of prefix agreement is that the axes you want to match aren’t always the prefix axes. Other array systems offer a feature similar to rank extension called ‘broadcasting’ that is a bit more flexible. For example, in the way it’s implemented in Numpy [num17] , an array of shape [A B 1 D] will match an array of shape [A B C D]. The process of broadcasting consists in inserting so-called ‘singleton dimensions’ (axes with length one) to align the axes that one wishes to match. You can think of prefix agreement as a particular case of broadcasting where the singleton dimensions are added to the end of the shorter shapes automatically.

A drawback of singleton broadcasting is that it muddles the distinction between a scalar and a vector of length 1. Sometimes, an axis of length 1 is no more than that, and if 2≠3 is a size mismatch, it isn’t obvious why 1≠2 shouldn’t be. To avoid this problem, ra:: supports broadcasting with undefined length axes (see insert).

ra::Big<double, 3> a({5, 3}, ra::_0);
ra::Big<double, 1> b({3}, 0.);
ra::Big<double, 3> c({1, 3}, ra::_0);

// b(?, i) += a(j, i) → b(i) = ∑ⱼ a(j, i) (sum columns)
b(ra::insert<1>) += a;

c = a; // 1 ≠ 5, still an agreement error

Still another way to align array axes is provided by the rank conjunction.

Even with axis insertion, it is still necessary that the axes one wishes to match are in the same order in all the arguments. Transposing the axes before extension is a possible workaround.

2.5 Cell iteration

map and for_each apply their operators to each element of their arguments; in other words, to the 0-cells of the arguments. But it is possible to specify directly the rank of the cells that one iterates over:

ra::Big<double, 3> a({5, 4, 3}, ra::_0);
for_each([](auto && b) { /* b has shape (5 4 3) */ }, iter<3>(a));
for_each([](auto && b) { /* b has shape (4 3) */ }, iter<2>(a));
for_each([](auto && b) { /* b has shape (3) */ }, iter<1>(a));
for_each([](auto && b) { /* b has shape () */ }, iter<0>(a)); // elements
for_each([](auto && b) { /* b has shape () */ }, a); // same as iter<0>(a); default

One may specify the frame rank instead:

for_each([](auto && b) { /* b has shape () */ }, iter<-3>(a)); // same as iter<0>(a)
for_each([](auto && b) { /* b has shape (3) */ }, iter<-2>(a)); // same as iter<1>(a)
for_each([](auto && b) { /* b has shape (4 3) */ }, iter<-1>(a)); // same as iter<2>(a)

In this way it is possible to match shapes in various ways. Compare

ra::Big<double, 2> a = {{1, 2, 3}, {4, 5, 6}};
ra::Big<double, 1> b = {10, 20};
ra::Big<double, 2> c = a * b; // multiply (each item of a) by (each item of b)

⇒ a = {{10, 20, 30}, {80, 100, 120}}

with

ra::Big<double, 2> a = {{1, 2, 3}, {4, 5, 6}};
ra::Big<double, 1> b = {10, 20, 30};
ra::Big<double, 2> c({2, 3}, 0.);
iter<1>(c) = iter<1>(a) * iter<1>(b); // multiply (each item of a) by (b)

⇒ a = {{10, 40, 90}, {40, 100, 180}}

Note that in this case we cannot construct c directly from iter<1>(a) * iter<1>(b), since the constructor for ra::Big matches its argument using (the equivalent of) iter<0>(*this). See iter for more examples.

Cell iteration is appropriate when the operations take naturally operands of rank > 0; for instance, the operation ‘determinant of a matrix’ is naturally of rank 2. When the operation is of rank 0, such as * above, there may be faster ways to rearrange shapes for matching (see The rank conjunction).

FIXME More examples.

2.6 Slicing

Slicing is an array extension of the subscripting operation. However, tradition and convenience have given it a special status in most array languages, together with some peculiar semantics that ra:: supports.

The form of the scripting operator A(i₀, i₁, ...) makes it plain that A is a function of rank(A) integer arguments⁹. An array extension is immediately available through map. For example:

ra::Big<double, 1> a = {1., 2., 3., 4.};
ra::Big<int, 1> i = {1, 3};
map(a, i) = 77.;

⇒ a = {1., 77., 3, 77.}

Just as with any use of map, array arguments are subject to the prefix agreement rule.

ra::Big<double, 2> a({2, 2}, {1., 2., 3., 4.});
ra::Big<int, 1> i = {1, 0};
ra::Big<double, 1> b = map(a, i, 0);

⇒ b = {3., 1.} // {a(1, 0), a(0, 0)}

ra::Big<int, 1> j = {0, 1};
b = map(a, i, j);

⇒ b = {3., 2.} // {a(1, 0), a(0, 1)}

The latter is a form of sparse subscripting.

Most array operations (e.g. +) are defined through map in this way. For example, A+B+C is defined as map(+, A, B, C) (or the equivalent map(+, map(+, A, B), C)). Not so for the subscripting operation:

ra::Big<double, 2> A {{1., 2.}, {3., 4.}};
ra::Big<int, 1> i = {1, 0};
ra::Big<int, 1> j = {0, 1};
// {{A(i₀, j₀), A(i₀, j₁)}, {A(i₁, j₀), A(i₁, j₁)}}
ra::Big<double, 2> b = A(i, j);

⇒ b = {{3., 4.}, {1., 2.}}

A(i, j, ...) is defined as the outer product of the indices (i, j, ...) with operator A, because this operation sees much more use in practice than map(A, i, j ...).

You may give fewer subscripts than the rank of the array. The full extent is assumed for the missing subscripts (cf all below):

ra::Big<int, 3> a({2, 2, 2}, {1, 2, 3, 4, 5, 6, 7, 8});
auto a0 = a(0); // same as a(0, ra::all, ra::all)
auto a10 = a(1, 0); // same as a(1, 0, ra::all)

⇒ a0 = {{1, 2}, {3, 4}}
⇒ a10 = {5, 6}

This supports the notion (see Rank polymorphism) that a 3-array is also an 2-array where the elements are 1-arrays themselves, or a 1-array where the elements are 2-arrays. This important property is directly related to the mechanism of rank extension (see Rank extension).

Besides, when the subscripts i, j, ... are scalars or integer sequences of the form (o, o+s, ..., o+s*(n-1)) (linear ranges), the subscripting can be performed inmediately at constant cost, and without needing to construct an expression object. This optimization is called beating.

ra:: isn’t smart enough to know when an arbitrary expression might be a linear range, so the following special objects are provided:

Special object: iota count [start:0 [step:1]] ¶

Create a linear range start, start+step, ... start+step*(count-1).

This can used anywhere an array expression is expected.

ra::Big<int, 1> a = ra::iota(4, 3 -2);

⇒ a = {3, 1, -1, -3}

Here, b and c are Views (see Containers and views).

ra::Big<int, 1> a = {1, 2, 3, 4, 5, 6};
auto b = a(iota(3));
auto c = a(iota(3, 3));

⇒ a = {1, 2, 3}
⇒ a = {4, 5, 6}

iota() by itself is an expression of rank 1 and undefined length. It must be used with other terms whose lengths are defined, so that the overall shape of the array expression can be determined. In general, iota<n>() is an array expression of rank n+1 that represents the n-th index of an array expression. This is similar to Blitz++’s TensorIndex.

ra:: offers the shortcut ra::_0 for ra::iota<0>(), etc.

ra::Big<int, 1> v = {1, 2, 3};
cout << (v - ra::_0) << endl; // { 1-0, 2-1, 3-2 }
// cout << (ra::_0) << endl; // error: undefined length
// cout << (v - ra::_1) << endl; // error: undefined length on axis 1
ra::Big<int, 2> a({3, 2}, 0);
cout << (a + ra::_0 - ra::_1) << endl; // {{0, -1, -2}, {1, 0, -1}, {2, 1, 0}}

When undefined length iota() is used as a subscript by itself, the result isn’t a View. This allows view(iota()) to match with expressions of different lengths, as in the following example.

ra::Big<int, 1> a = {1, 2, 3, 4, 5, 6};
ra::Big<int, 1> b = {1, 2, 3};
cout << (b + a(iota())) << endl; // a(iota()) is not a View

-| 3
2 4 6

Note the difference between

• ra::iota<3>() — an expression of rank 4 and undefined length, representing a linear sequence over the tensor index of axis 3
• ra::iota(3) ≡ ra::iota<0>(3) — an expression of rank 1, representing the sequence 0, 1, 2.

Special object: all ¶

Create a linear range 0, 1, ... (nᵢ-1) when used as a subscript at the i-th place of a subscripting expression. This might not be the i-th argument; see insert, dots.

This object cannot stand alone as an array expression. All the examples below result in View objects:

ra::Big<int, 2> a({3, 2}, {1, 2, 3, 4, 5, 6});
auto b = a(ra::all, ra::all); // (1) a view of the whole of a
auto c = a(iota(3), iota(2)); // same as (1)
auto d = a(iota(3), ra::all); // same as (1)
auto e = a(ra:all, iota(2)); // same as (1)
auto f = a(0, ra::all); // first row of a
auto g = a(ra::all, 1); // second column of a
auto g = a(ra::all, ra::dots<0>, 1); // same

all is a special case (dots<1>) of the more general object dots.

Special object: dots<n> ¶

Equivalent to as many instances of ra::all as indicated by n, which must not be negative. Each instance takes the place of one argument to the subscripting operation.

If n is defaulted (dots<>), all available places will be used; this can only be done once in a given subscript list.

This object cannot stand alone as an array expression. All the examples below result in View objects:

ra::Big<int, 3> a({3, 2, 4}, ...);
auto h = a(); // all of a

auto b = a(ra::all, ra::all, ra::all); // (1) all of a
auto c = a(ra::dots<3>); // same as (1)
auto d = a(ra::all, ra::dots<2>); // same as (1)
auto e = a(ra::dots<2>, ra::all); // same as (1)
auto f = a(ra::dots<>); // same as (1)

auto j0 = a(0, ra::dots<2>); // first page of a
auto j1 = a(0); // same
auto j2 = a(0, ra::dots<>); // same

auto k0 = a(ra::all, 0); // first row of a
auto k1 = a(ra::all, 0, ra::all); // same
auto k2 = a(ra::all, 0, ra::dots<>); // same
auto k3 = a(ra::dots<>, 0, ra::all); // same
// auto k = a(ra::dots<>, 0, ra::dots<>); // error

auto l0 = a(ra::all, ra::all, 0); // first column of a
auto l1 = a(ra::dots<2>, 0); // same
auto l2 = a(ra::dots<>, 0); // same

This is useful when writing rank-generic code, see examples/maxwell.cc in the distribution for an example.

The following special objects aren’t related to linear ranges, but they are meant to be used in a subscript context. Using them in other contexts will result in a compile time error.

Special object: len ¶

Represents the length of the i-th axis of a subscripted expression, when used at the i-th place of a subscripting expression.

This works like end in Octave/Matlab, but note that ra:: indices begin at 0, so the last element of a vector a is a(ra::len-1).

ra::Big<int, 2> a({10, 10}, 100 + ra::_0 - ra::_1);
auto a0 = a(ra::len-1); // last row of a; ra::len is a.len(0)
auto a1 = a(ra::all, ra::len-1); // last column a; ra::len is a.len(1)
auto a2 = a(ra::len-1, ra::len-1); // last element of last row; the first ra::len is a.len(0) and the second one is a.len(1)
auto a3 = a(ra::all, ra::iota(2, ra::len-2)); // last two columns of a
auto a4 = a(ra::iota(ra::len/2, 1, 2)); // odd rows of a
a(ra::len - std::array {1, 3, 4}) = 0; // set to 0 the 1st, 3rd and 4th rows of a, counting from the end

ra::Big<int, 3> b({2, 3, 4}, ...);
auto b0 = b(ra::dots<2>, ra::len-1); // ra::len is a.len(2)
auto b1 = b(ra::insert<1>, ra::len-1); // ra::len is a.len(0)

Special object: insert<n> ¶

Inserts n new axes at the subscript position. n must not be negative.

The new axes have step 0 and undefined length, so they will match any length on those axes by repeating items. insert objects cannot stand alone as an array expression. The examples below result in View objects:

auto h = a(insert<0>); // same as (1)
auto k = a(insert<1>); // shape [undefined, 3, 2]

insert<n> main use is to prepare arguments for broadcasting. In other array systems (e.g. Numpy) broadcasting is done with singleton dimensions, that is, dimensions of length one match dimensions of any length. In ra:: singleton dimensions aren’t special, so broadcasting requires the use of insert. For example:

ra::Big<int, 1> x = {1, 10};
// match shapes [2, U, U] with [U, 3, 2] to produce [2, 3, 2]
cout << x(ra::all, ra::insert<2>) * a(insert<1>) << endl;

-| 2 3 2
1 2
3 4
5 6

10 20
30 40
50 60

We were speaking earlier of the outer product of the subscripts with operator A. Here’s a way to perform the actual outer product (with operator *) of two Views, through broadcasting. All three lines are equivalent. See from for a more general way to compute outer products.

cout << (A(ra::dots<A.rank()>, ra::insert<B.rank()>) * B(ra::insert<A.rank()>, ra::dots<B.rank()>)) << endl;
cout << (A(ra::dots<>, ra::insert<B.rank()>) * B(ra::insert<A.rank()>, ra::dots<>)) << endl; // default dots<>
cout << (A * B(ra::insert<A.rank()>)) << endl; // index elision + prefix matching

• Subscripting and rank-0 views

using T = std::complex<double>;
int f(T &);
Big<T, 2> a({2, 3}, 0); // ct rank
Big<T> b({2, 3}, 0); // rt rank

cout << a(0, 0).real_part() << endl; // ok, a(0, 0) returns complex &
// cout << b(0, 0).real_part() << endl; // error, View<T> has no member real_part
cout << ((T &)(b(0, 0))).real_part() << endl; // ok, manual conversion to T &

cout << f(b(0, 0)) << endl; // ok, automatic conversion from View<T> to T &
// cout << f(a(0)) << endl; // compile time error, conversion failed since ct rank of a(0) is not 0
// cout << f(b(0)) << endl; // runtime error, conversion failed since rt rank of b(0) is not 0

2.7 Functions

You don’t need to use map every time you want to do something with arrays in ra::. A number of array functions are already defined.

• Standard scalar operations
• Conditional operations
• Special operations
• Elementwise reductions
• Special reductions
• Shortcut reductions

cout << exp(ra::Small<double, 3> {4, 5, 6}) << endl;

  -| 54.5982 148.413 403.429

ra::Big<double> x ...
ra::Big<double> y = where(x>0, expensive_expr_1(x), expensive_expr_2(x));

ra::Big<double> y = map([](auto && w, auto && t, auto && f) -> decltype(auto) { return w ? t : f; }
                        x>0, expensive_expr_1(x), expensive_function_2(x));

ra::Big<int, 1> o = {};
ra::Big<int, 1> e = {};
ra::Big<int, 1> n = {1, 2, 7, 9, 12};
ply(where(odd(n), map([&o](auto && x) { o.push_back(x); }, n), map([&e](auto && x) { e.push_back(x); }, n)));
cout << "o: " << ra::noshape << o << ", e: " << ra::noshape << e << endl;

-| o: 1 7 9, e: 2 12

template <class A>
inline auto op_reduce(A && a)
{
    T c = op_default;
    for_each([&c](auto && a) { c = op(c, a); }, a);
    return c;
}

    ra::Big<double, 2> a({m, n}, ...);

    ra::Big<double, 1> sum_rows({n}, 0.);
    iter<1>(sum_rows) += iter<1>(a);
    // for_each(ra::wrank<1, 1>([](auto & c, auto && a) { c += a; }), sum_rows, a) // alternative
    // sum_rows += transpose<1, 0>(a); // another

    ra::Big<double, 1> sum_cols({m}, 0.);
    sum_cols += a;

2.8 The rank conjunction

We have seen in Cell iteration that it is possible to treat an r-array as an array of lower rank with subarrays as its elements. With the iter<cell rank> construction, this ‘exploding’ is performed (notionally) on the argument; the operation of the array expression is applied blindly to these cells, whatever they turn out to be.

for_each(my_sort, iter<1>(A)); // (in ra::) my_sort is a regular function, cell rank must be given
for_each(my_sort, iter<0>(A)); // (in ra::) error, bad cell rank

The array language J has instead the concept of verb rank. Every function (or verb) has an associated cell rank for each of its arguments. Therefore iter<cell rank> is not needed.

for_each(sort_rows, A); // (not in ra::) will iterate over 1-cells of A, sort_rows knows

ra:: doesn’t have ‘verb ranks’ yet. In practice one can think of ra::’s operations as having a verb rank of 0. However, ra:: supports a limited form of J’s rank conjunction with the function wrank.

This is an operator that takes one verb (such operators are known as adverbs in J) and produces another verb with different ranks. These ranks are used for rank extension through prefix agreement, but then the original verb is used on the cells that result. The rank conjunction can be nested, and this allows repeated rank extension before the innermost operation is applied.

A standard example is ‘outer product’.

ra::Big<int, 1> a = {1, 2, 3};
ra::Big<int, 1> b = {40, 50};
ra::Big<int, 2> axb = map(ra::wrank<0, 1>([](auto && a, auto && b) { return a*b; }),
                            a, b)

⇒ axb = {{40, 80, 120}, {50, 100, 150}}

It works like this. The verb ra::wrank<0, 1>([](auto && a, auto && b) { return a*b; }) has verb ranks (0, 1), so the 0-cells of a are paired with the 1-cells of b. In this case b has a single 1-cell. The frames and the cell shapes of each operand are:

a: 3 |
b:   | 2

Now the frames are rank-extended through prefix agreement.

a: 3 |
b: 3 | 2

Now we need to perform the operation on each cell. The verb [](auto && a, auto && b) { return a*b; } has verb ranks (0, 0). This results in the 0-cells of a (which have shape ()) being rank-extended to the shape of the 1-cells of b (which is (2)).

a: 3 | 2
b: 3 | 2

This use of the rank conjunction is packaged in ra:: as the from operator. It supports any number of arguments, not only two.

ra::Big<int, 1> a = {1, 2, 3};
ra::Big<int, 1> b = {40, 50};
ra::Big<int, 2> axb = from([](auto && a, auto && b) { return a*b; }), a, b)

⇒ axb = {{40, 80, 120}, {50, 100, 150}}

Another example is matrix multiplication. For 2-array arguments C, A and B with shapes C: (m, n) A: (m, p) and B: (p, n), we want to perform the operation C(i, j) += A(i, k)*B(k, j). The axis alignment gives us the ranks we need to use.

C: m |   | n
A: m | p |
B:   | p | n

First we’ll align the first axes of C and A using the cell ranks (1, 1, 2). The cell shapes are:

C: m | n
A: m | p
B:   | p n

Then we’ll use the ranks (1, 0, 1) on the cells:

C: m |   | n
A: m | p |
B:   | p | n

The final operation is a standard operation on arrays of scalars. In actual ra:: syntax:

ra::Big A({3, 2}, {1, 2, 3, 4, 5, 6});
ra::Big B({2, 3}, {7, 8, 9, 10, 11, 12});
ra::Big C({3, 3}, 0.);
for_each(ra::wrank<1, 1, 2>(ra::wrank<1, 0, 1>([](auto && c, auto && a, auto && b) { c += a*b; })), C, A, B);

⇒ C = {{27, 30, 33}, {61, 68, 75}, {95, 106, 117}}

Note that wrank cannot be used to transpose axes in general.

I hope that in the future something like C(i, j) += A(i, k)*B(k, j), where i, j, k are special objects, can be automatically translated to the requisite combination of wrank and perhaps also transpose. For the time being, you have to align or transpose the axes yourself.

2.9 Compatibility

• Using other C and C++ types with ra::
• Using ra:: types with the STL
• Using ra:: types with other libraries

std::vector<int> x = {1, 2, 3};
ra::Small<int, 3> y = {6, 5, 4};
cout << (x-y) << endl;

-| -5 -3 -1

std::vector<int> x = {1, 2, 3};
// cout << sum(x) << endl; // error, sum not found
cout << sum(ra::start(x)) << endl;
cout << ra::sum(x) << endl;

-| 6

int p[] = {1, 2, 3};
int * z = p;
ra::Big<int, 1> q {1, 2, 3};
q += p; // ok, q is ra::, p is foreign object with shape info
ra::start(p) += q; // can't redefine operator+=(int[]), foreign needs ra::start()
// z += q; // error: raw pointer needs ra::ptr()
ra::ptr(z) += p; // ok, shape is determined by foreign object p

ra::Big<int, 2> A = {{3, 0, 0}, {4, 5, 6}, {0, 5, 6}};
std::accumulate(ra::begin(A), ra::end(A), 0); // or just sum(A)

⇒ 29

cout << std::ranges::count(range(A>3), true) << endl; // or sum(cast<int>(A>3))

⇒ 5

// count rows with 0s in them
cout << std::ranges::count_if(range(iter<1>(A)), [](auto const & x) { return any(x==0); }) << endl;

⇒ 2

ra::Big<int> x {3, 2, 1}; // x is a Container
auto y = x(); // y is a View on x
// std::sort(ra::begin(y), ra::end(y)); // error: begin(y) is not std::random_access_iterator
std::sort(ra::begin(x), ra::end(x)); // ok, we know x is stored in row-major order

⇒ x = {1, 2, 3}

2.10 Extension

• New scalar types
• New array operations

struct W { int x; }
ra::Big<W, 2> w = {{ {4}, {2} }, { {1}, {3} }};
cout << W(1, 1).x << endl;
cout << amin(map([](auto && x) { return w.x; }, w)) << endl;

-| 3
   1

template <> constexpr bool ra::is_scalar_def<W> = true;
...
W ww {11};
for_each([](auto && x, auto && y) { cout << (x.x + y.y) << " "; }, w, ww); // ok

-| 15 13 12 14

struct U { int x; }
U uu {11};
for_each([](auto && x, auto && y) { cout << (x.x + y.y) << " "; }, w, uu); // error: can't find ra::start(U)

return_type my_fun(...) { };
...
namespace ra {
template <class ... A> inline auto
my_fun(A && ... a)
{
    return map(::my_fun, std::forward<A>(a) ...);
}
} // namespace ra

namespace ra {
template <class ... A> inline auto
my_fun(A && ... a)
{
    return map([](auto && ... a) { return ::my_fun(a ...); }, std::forward<A>(a) ...);
}
} // namespace ra

2.11 Error handling

ra:: tries to detect a many errors as possible at compile time, but some errors, such as mismatch between expressions of runtime shape or out of range indices, aren’t apparent until runtime.

Runtime error handling in ra:: is controlled by two macros.

• RA_ASSERT(cond, ...) — check that cond evaluates to true in the ra:: namespace. The other arguments are informative.
• RA_DO_CHECK — must have one of the values 0, 1, or 2.

They work as follows:

• If RA_DO_CHECK is 0, runtime checks are skipped.
• If RA_DO_CHECK is not 0, runtime checks are done.
  • If RA_ASSERT is defined, using RA_ASSERT.
  • If RA_ASSERT isn’t defined, the method depends on the value of RA_DO_CHECK. The two options are 1 (plain assert) and 2 (prints the informative arguments and aborts). Other values are an error.

ra:: contains uses of assert for checking invariants or for sanity checks that are separate from uses of RA_ASSERT. Those can be disabled in the usual way with -DNDEBUG, but note that -DNDEBUG will also disable any asserts that are a result of RA_DO_CHECK=1.

The performance cost of the runtime checks depends on the program. Without custom RA_ASSERT, RA_DO_CHECK=1 usually has an acceptable cost, but RA_DO_CHECK=2 may be more expensive. The default is RA_DO_CHECK=1.

The following example shows how errors might be reported depending on RA_DO_CHECK.

ra::Big<int, 2> a({10, 3}, 1);
ra::Big<int, 2> b({40, 3}, 2);
cout << sum(a+b) << endl;

• RA_DO_CHECK=2

  *** ra::./ra/expr.hh:389,17 (check()) Mismatched shapes [10 3] [40 3]. ***}
  

• RA_DO_CHECK=1

  ./ra/expr.hh:389: constexpr ra::Match<checkp, std::tuple<_UTypes ...>, std::tuple<std::integral_constant<int, I>...> >::Match(P ...)
  [with bool checkp = true; P = {ra::Cell<int, const ra::SmallArray<ra::Dim, std::tuple<std::integral_constant<int, 2> >,
  std::tuple<std::integral_constant<int, 1> >, std::tuple<ra::Dim, ra::Dim> >&, std::integral_constant<int, 0> >, ra::Cell<int, const
  ra::SmallArray<ra::Dim, std::tuple<std::integral_constant<int, 2> >, std::tuple<std::integral_constant<int, 1> >, std::tuple<ra::Dim,
  ra::Dim> >&, std::integral_constant<int, 0> >}; int ...I = {0, 1}]: Assertion `check()' failed.
  

You can redefine RA_ASSERT to something that is more appropriate for your program. For example, if you run ra:: code under a shell, an abort may not be acceptable. examples/throw.cc shows how to throw a user-defined exception instead.

4.1 Reuse of expression objects

Expression objects are meant to be used once. This applies to anything produced with ra::map, ra::iter, ra::start, or ra::ptr. Reuse errors are not checked. For example:

ra::Big<int, 2> B({3, 3}, ra::_1 + ra::_0*3); // {{0 1 2} {3 4 5} {6 7 8}}
std::array<int, 2> l = { 1, 2 };
cout << B(ra::ptr(l), ra::ptr(l)) << endl; // ok => {{4 5} {7 8}}
auto ll = ra::ptr(l);
cout << B(ll, ll) << endl; // ??

4.2 Assignment to views

FIXME With rt-shape containers (e.g. Big), operator= replaces the left hand side instead of writing over its contents. This behavior is inconsistent with View::operator= and is there only so that istream >> container may work; do not rely on it.

4.3 View of const vs const view

FIXME Passing view arguments by reference

4.4 Rank extension in assignments

Assignment of an expression onto another expression of lower rank may not do what you expect. This example matches a and 3 [both of shape ()] with a vector of shape (3). This is equivalent to {a=3+4; a=3+5; a=3+6;}. You may get a different result depending on traversal order.

int a = 0;
ra::scalar(a) = 3 + ra::Small<int, 3> {4, 5, 6}; // ?

  ⇒ a = 9

Compare with

int a = 0;
ra::scalar(a) += 3 + ra::Small<int, 3> {4, 5, 6}; // 0 + 3 + 4 + 5 + 6

  ⇒ a = 18

4.5 Performance pitfalls of rank extension

In the following example where b has its shape extended from (3) to (3, 4), f is called 12 times, even though only 3 calls are needed if f doesn’t have side effects. In such cases it might be preferrable to write the outer loop explicitly, or to do some precomputation.

ra::Big<int, 2> a = {{1, 2, 3, 4}, {5, 6, 7, 8} {9, 10, 11, 12}};
ra::Big<int, 1> b = {1, 2, 3};
ra::Big<int, 2> c = map(f, b) + a;

4.6 Initialization of nested types

This doesn’t do what you probably think it does.

ra::Big<int2, 0> b({}, int2 {1, 3}); // b = {int2 {3, 3}} (*)

What happens here is that rank-0 b is prefix matched to rank-1 int2 {1, 3} which results in two assignments, first b ← 1 and then b ← 3. Compare with

ra::Big<int2, 0> b({}, ra::scalar(int2 {1, 3})); // b = {int2 {1, 3}} @c [ma116]

The use of scalar forces b to match the whole of int2 {1, 3} as a rank-0 object. After that, the single element in b is matched to int2 {1, 3}.

Assignment or initialization from higher rank to lower rank is something of a footgun. Assignments have some legitimate uses, but initializations like (*) may become errors in future versions of ra::.

4.7 Chained assignment

FIXME When a=b=c works, it operates as b=c; a=b; and not as an array expression.

4.8 Unregistered scalar types

FIXME View<T, N> x; x = T() fails if T isn’t registered as is_scalar.

1. Item 0
2. Item 1
3. Item 2

4.9 The cost of RA_DO_CHECK

FIXME

5.1 Headers

The header structure of ra:: is as follows.

• tuples.hh - Generic macros and tuple library.
• base.hh - Basic types and concepts, introspection.
• expr.hh - Expression template nodes.
• ply.hh - Traversal, including I/O.
• small.hh - Array type with compile time dimensions.
• big.hh - Array type with run time dimensions.
• ra.hh - Functions and operators, main header.
• test.hh - (accessory) Testing/benchmarking library.
• dual.hh - (accessory) Dual number type and operations.

5.2 Type hierarchy

Some of the categories below are C++20 ‘concepts’, some are still informal.

• Container — Big, Shared, Unique, Small

  These are array types that own their data in one way or another.

• View — ViewBig, ViewSmall

  These are array views into data in memory, which may be writable. Any of the Container types can be treated as a View, but one may also create Views into memory that has been allocated independently.

• Iterator — CellBig, CellSmall, Iota, Ptr, Scalar, Expr, Pick

  This is a traversable object. Iterators are accepted by all the array functions such as map, for_each, etc. map produces an Iterator itself, so most array expressions are Iterators. Iterators are created from Views and from certain foreign array-like types primarily through the function start. This is done automatically when those types are used in array expressions.

  Iterators have two traversal functions: .adv(k, d), moves the iterator along any dimension k, and .mov(d), is used on linearized views of the array. The methods len(), step(), keep() are used to determine the extent of these linearized views. In this way, a loop involving Iterators can have its inner loop unfolded, which is faster than a nested loop, especially if the inner dimensions of the loop are small.

  Iterators also provide an at(i ...) method for random access to any element.

5.3 Term agreement

The execution of an expression template begins with the determination of its shape — the length of each of its dimensions. This is done recursively by traversing the terms of the expression. For a given dimension k≥0, terms that have rank less or equal than k are ignored, following the prefix matching principle. Likewise terms where dimension k has undefined length (such as iota() or dimensions created with insert) are ignored. All the other terms must match.

Then we select a order of traversal. ra:: supports ‘array’ orders, meaning that the dimensions are sorted in a certain way from outermost to innermost and a full dimension is traversed before one advances on the dimension outside. However, currently (v28) there is no heuristic to choose a dimension order, so traversal always happens in row-major order (which shouldn’t be relied upon). ply_ravel will unroll as many innermost dimensions as it can, and in some cases traversal will be executed as a flat loop.

Finally we select a traversal method. ra:: has two traversal methods: ply_fixed can be used when the rank and the traversal order are known at compile time, and ply_ravel can be used in the general case.

5.4 Traversal

5.5 Introspection

The following functions are available to query the properties of ra:: objects.

Function: rank_t rank e ¶

Return the rank of expression e.

Function: array shape e ¶

Function: expr shape e k ¶

The first form returns the shape of expression e as an array. The second form returns the length of axis k, i.e. shape(e)[k] ≡ shape(e, k). It is possible to use len in k to mean the rank of e.

    ra::Small<int, 2, 3, 4> A;
    cout << shape(A, ra::len-1) << endl; // length of last axis

-| 4

If k is not a scalar, shape(e, k) returns a map of k over shape(e, ·).

    cout << shape(A, ra::iota(2, ra::len-2)) << endl; // last two lengths

-| 3 4

Note that the return type array of shape(e) might be a foreign type (such as std::array or std::vector) instead of a ra:: type. In that case, shape(A)(ra::iota(2, ra::len-2)) will not work.

5.6 Building the test suite

ra:: uses its own test and benchmark suites and it comes with three kinds of tests: proper tests, benchmarks, and examples. To build and run them all, run the following from the top directory:

CXXFLAGS=-O3 scons

or

CXXFLAGS=-O3 cmake . && make && make test

ra:: is highly dependent on the optimization level and the test suite may run much slower with -O0 or -O1.

The SCons/CMake scripts accept the following options:

• scons --no-sanitize / cmake -DSANITIZE=0

  Remove the sanitizer flags. This is necessary because appending -fno-sanitize=all to CXXFLAGS, as in -fsanitize=... -fno-sanitize=all, doesn’t completely remove the sanitizers. The default is to build with sanitizers.

• scons --blas

  Try to use BLAS in the benchmarks. The CMake script tries to detect BLAS automatically, but the SCons script doesn’t, hence this option.

6.1 Error messages

FIXME

6.2 Reductions

FIXME

6.3 Etc

FIXME