# New mathematics for computer performance: array algebra and cost functions

Gaétan Hains, Lenore Mullin

► To cite this version:

Gaétan Hains, Lenore Mullin. New mathematics for computer performance: array algebra and cost functions. 2026. hal-05442652

HAL Id: hal-05442652

https://hal.science/hal-05442652v1

Preprint submitted on 5 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Article

# New mathematics for computer performance: array algebra and cost functions

Gaétan Hains¹,†,‡ and Lenore Mullin²,*

1 LACL, Université Paris-Est Créteil (UPEC) France; gaetan.hains@u-pec.fr

2 University at Albany (SUNY) NY, USA; lmullin@albany.edu

* Correspondence: gaetan.hains@u-pec.fr;

† Current address: LACL, UPEC, 61 avenue du Général de Gaulle 94000 Créteil, France.

‡ These authors contributed equally to this work.

**Abstract:** MoA (mathematics of arrays) is a theory of parallel operations on arrays that can describe all known algorithms in linear algebra, signal processing, HPC because they are based on primitive-recursion and array shapes. Mapping parallel algorithms to computer architectures remains more of an art than a science and specific mathematical techniques are needed to provide a basis for performance evaluation at a level abstract enough to constitute an experimental science.

In this paper we present a methodology for implicit parallel programming with MoA expressions. Then we formalize the Moa operators in an operational semantics that relates them to the linear space of memory elements in computer architecture.

Then we define a theory of execution costs that can generalize all known concepts of this kind, is based on classical operations research and is formally related to MoA-based parallel programming. It constitutes the basis for a formalized and mechanizable approach to performance prediction, portability and maximization.

**Keywords:** Multidimensional arrays, tensors, algebra, mathematics of arrays, high-performance programming, automata, semantics, time- and energy measures

## 1. Introduction and related work

The productivity and reliability of high-performance computing (HPC) is very low due to complex and ever-evolving hardware architectures, compilers and the general lack of performance predictability.

Languages and libraries used for producing parallel code are based on diverse and historically unstable concepts: APL[1] was declarative and array-centered but is no more in use, FORTRAN[2] has been relatively stable but requires the addition of explicitly-parallel libraries like MPI; high-performance Fortran (HPF) [3] was an attempt at removing this limitation by adding a process-mapping concept to Fortran, C [4] remains a stable alternative to Fortran for this purpose but shares the same defects as the libraries it is extended with (MPI etc.). The Message passing interface (MPI) [5] proposes a clear concept of local process and global operations. But its semantics is non-deterministic and its syntax scoping is completely orthogonal to the process structures, a defect that obscures the parallel meaning of even the simplest codes. OpenMP [6] only applies to shared-memory computing so does not cover the scopes of MPI or CUDA. CUDA [7] is GPU-specific and proprietary. Other libraries emerge regularly that are more or less general purpose, for example OpenACC [8] uses annotations to target hybrid CPU-GPU architectures, while Pytorch [9] is specific to machine learning.

Received:

Revised:

Accepted:

Published:

**Citation:** Hains, G.; Mullin, L.; New mathematical tools for computer performance and energy management. *Mathematics* **2025**, *1*, 0.

https://doi.org/

**Copyright:** © 2026 by the authors.
Submitted to *Mathematics* for possible
open access publication under the
terms and conditions of the Creative
Commons Attribution (CC BY)
license (https://creativecommons.
org/licenses/by/4.0/).
One could argue that sequential programming languages are also constantly evolving: Fortran, Algol, C, Java, Python have all been popular at different periods and propose various programming concepts. But they are all based on the well-proven concepts of sequential computing. The panorama of *parallel* programming languages and libraries is not converging on a single concept that can combine simplicity of expression and predictable performance. In our view this is why high-performance computing remains a niche specialty, without serious software engineering. Today, even with the enormous popularity of machine learning, few data scientists understand the operators that built their models, or the effects of hardware on their execution speed.

The increasing variety of hardware platforms and hardware structures has worsened this overall parallel software crisis. Vector units, multicore CPUs, GPUs, TPUs, multinode systems etc.: all have different compute-to-communication speed ratios, different sizes so that the best mapping of a given algorithm to the target hardware is a kind of moving target for the developer and for the compiler. Yet all hardware structures have one fundamental feature in common: they are built from *arrays* of repeated elements: vector units, cores, nodes, execution units etc. This observation is the basis of our proposal for array-based parallel programming.

Whatever the language or library that is used, parallel programming's only purpose is to accelerate computation. But that can only happen if computation time is reasonably predictable. That requirement leads to the concept of so-called *cost models*. A cost model is theory that allows the production of a run-time estimate from a description of the algorithm combined with a summary of the target architecture. Two famous families of cost models are the roofline model and the bulk-synchronous parallel (BSP) performance model. We now outline their diverging advantages and shortcomings.

* BSP computing [10,11] assumes that execution is a sequence of asynchronous parallel sets of traces, interleaved with global exchange and synchronization phases. It benchmarks the target architecture for speed, bandwidth and global synchronization delay, and combines this information with the algorithm's computation graph to produce a reliable and extensible cost i.e. computation time estimate. The flat parallel structure of BSP made it unsuitable for mapping to hierarchical systems so that it could only represent one-level multicore / multinode or a single GPU's elements. But that restriction has been lifted with several multi-level versions of BSP [12–14]. The work of Louergue, Matsuzaki et al [15] has derived a class of BSP algorithms mathematically, but BSP's main drawback remains the lack of an algebra of programs so that programs must be written in an ad hoc manner. Moreover it assumes a fixed series-parallel (or nested series-parallel) form of computation which may not fit all future architectures.

* The roofline model [16,17] is a performance analysis technique for parallel executions. It relates global memory-access intensity with computation Flop rate. It supports performance tuning by giving hints at what part of the source code is responsible for not reaching maximal Flop/s rate. But that is precisely its main defect: it is blind to any algorithm optimization, such as blocking or Strassen-like matrix multiplication, that would reduce the total number of operations. Like BSP, the roofline model has been given a multi-level extension [18]. In summary the roofline model is relevant and useful but without a concept of source program so it is useless for parallel software engineering.

* Another approach to parallel programming is polyhedral compilation [19] with its well-studied theory and practice tiling and mapping source code to virtual processors. Its properties are dual to those of the BSP model. The source code is modeled as a
set of loop nests, so the model has an obvious link to software engineering, but the mapping to virtual processors is done without any cost model.

To go beyond the deficiencies of existing theories and techniques, we are building a declarative programming approach to HPC whereby a pure mathematical formulation of the computation is converted to executable code while maximizing performance and making it predictable. The source description uses an algebra of multi-dimensional arrays called *A Mathematics of Arrays (MoA)* [20–26] that provides a space of transformations allowing the automatic exploration of equivalent programs. Reliable code generation uses formal semantics for MoA and the simplest possible language fragments: pure loop nests in C language, and annotations or instructions for vectorization-, multicore- or multinode execution. Last and most importantly, performance prediction requires execution cost models that produce good execution time and energy consumption estimates, but ones that are mathematically related to the original mathematical formulation i.e. MoA expression. The success of this approach to HPC can produce massive improvements and enormous economic gains: low-coding to no-coding for the scientist or engineer, portable and reliable high-performance on a wide variety or computing hardware and lack of programming errors by design.

In [27] we proposed a proof of concept and software architecture for the declarative programming approach we advocate. If coupled in a logically reliable way to a generative AI system [28], it will become a *no-programming* method for HPC. But the cost prediction that was used in [27], despite being faithful, scalable and so far portable, was not mathematically connected to the MoA algebra.

In this paper we fill this gap by providing mathematical results that connect MoA to execution costs in a manner that will allow full and simple automation: algebraic correspondences between arrays of data and operations on the one hand, and hierarchies of arrays of hardware elements on the other. This work unites concepts from algorithm design, parallel programming, and parallel architectures whose disconnected use led to the productivity and reliability crisis in HPC.

The next sections describe: the algebra of arrays, its use for declarative high-performance computing, its semantics and equational transforms, a framework for execution cost models and theorems that preserve cost information while transforming MoA expressions.

## 2. Mathematics of arrays (MoA)

### 2.1. *Principles of MoA*

The idea of Mathematics of Arrays (MoA) considers the processing of data stored in multi-dimensional arrays (n-d arrays) by modelling the data array as a vector (1-d array). The mathematical operations (add/subtract, multiply/divide) then manipulate the data by considering indexing of the elements within the vector as depicted in Figure 1. Here, the n-d array (a tensor) is flattened so that the data is stored in contiguous memory that is then accessed by the suitable PU that performs the memory access and computation operations.
<figure><img src="image_2.png" /><figcaption>Figure 1. Array processing using MoA: Data and hardware perspective.</figcaption></figure>

Consider a small example as shown in Figure 2 to illustrate how MoA considers multi-dimensional arrays, where the dimension can be any dimension required:

<figure><img src="image_3.png" /><figcaption>Figure 2. MoA: Array flattening and placed in contiguous memory.</figcaption></figure>

Here a 2x3 matrix holds integer values 1, 2, 3, 4, 5, 6. This will be defined in software, for example in C as:

int example[2][3] = {{1, 2, 3}, {4, 5, 6}};

## 2.2. MoA: The Basics

Mathematically, the idea of MoA and Psi Calculus, starts with a scalar, $\sigma$, and its shape $\rho\sigma$, an empty vector, denoted by $<>$ or $\Theta$. Since a vector is an array, it has a shape, $(\rho <> ) \equiv < 0 >$, the one element vector containing the integer 0, denoting a scalar is a 0-dimensional array. Thus, from the shape, we can determine dimensionality, $\delta\sigma \equiv (\rho(\rho\sigma))[0] \equiv < 0 > [0] \equiv 0$ and the total number of components, $\pi(\rho\sigma) \equiv \pi <> \equiv 1$, in the array. Algorithms on shapes describe how to index arrays with Psi and are defined, such that, Psi takes a valid index vector, or an array of valid index vectors, and an array as arguments. For example, for a scalar, $\sigma$, $<>$, the empty vector, is the only valid index. Next, 2 scalars (0-d) and operations between them can be considered, with an extension to operations with n-dimensional (n-d) arguments.

Scalar operations are at the heart of computation, $\sigma_l f \sigma_r$, and in general for n-d arrays, $\xi_l f \xi_r$, where f is an arbitrary scalar function <fnref n="1" />.

Thus, in general,

$$ \begin{gathered} \forall \vec{i} \ni 0 \le \vec{i} < \rho \xi_l \\ \vec{i} \psi(\xi_l f \xi_r) \equiv (\vec{i} \psi \xi_l) f (\vec{i} \psi \xi_r) \end{gathered} \qquad (1) $$

<footnote n="1">¹ When f is commutative, $\vec{i}\psi(\xi_l f \xi_r) \equiv \vec{i}\psi(\xi_r f \xi_l)$. In addition, if g is a scalar operation, $\vec{i}\psi(\xi_l f(\xi_l g \xi_r)) \equiv \vec{i}\psi((\xi_l f \xi_l) g \xi_r)$ and associative.</footnote>
Here, $\vec{i}$ is a valid full index vector. That is, $\tau\vec{i} = \delta\xi$. The * notation denotes that all components of the index vector are within the bounds of its corresponding shape component.

That is, indexing distributes over scalar operations, or in compiler optimization terms, loop fusion.

For a scalar there is only one valid index, $\vec{i}$ is $<>$, the empty vector.

$$<> \psi(\sigma_l f \sigma_r) \equiv (<> \psi \sigma_l) f (<> \psi \sigma_r)$$

The equations above denote their Denotational Normal Form (DNF) and are in terms of cartesian coordinates.

$$\equiv (\text{rav } \sigma_l)[\gamma_{row}(<>;<>)] f (\text{rav } \sigma_r)[\gamma_{row}(<>;<>)]$$

Here we transform the DNF to its Operational Normal Form (ONF), or how to build the design, by flattening the array and relating that to offsets from the start of the array. $\gamma_{row}$ gets the offset using index and shape.

$$\equiv (\text{rav } \sigma_l)[0] f (\text{rav } \sigma_r)[0]$$

Here, the definition of $\gamma_{row}$ is applied in order to get 0.

$$\equiv (@\sigma_l + 0) f (@\sigma_r + 0)$$

This ONF now represents a generic form that uses bracket notation to denote how to build the code, through $\gamma$ and pseudo-code, i.e., $\text{rav}(\vec{i}\psi\xi) \equiv (\text{rav } \xi)[\gamma(\vec{i};\rho\xi)]$. With a family of gamma, $\gamma$, functions, e.g., row major($\gamma_{row}$), column major($\gamma_{column}$), a Cartesian index is related to an offset of the array in memory laid out contiguously. **rav** flattens an array based on its layout. Subscripts relate to left or right arguments and superscripts specify dimensionality. With this, introduce the important identify that is true in general for n-d arrays, denoted by $\xi$:

$$(\iota(\rho\xi))\psi\xi \equiv \xi$$

This means, with an array's shape, $\rho\xi$, generate an array of indices, $\iota(\rho\xi)$. Then, using that array as an argument to $\text{Psi}$, the original array $\xi$, is returned.

In the scalar case, where $\sigma$ denotes a scalar,

$$((\iota(\rho\sigma))\psi\sigma) \equiv \sigma$$

The shape of a scalar is the empty vector $<>$ and the only valid index of a scalar.

$$(\iota<>)(\psi\sigma) \equiv \sigma$$

$$<> \psi\sigma \equiv \sigma$$

Thus, as long as we can get shapes, e.g., $\rho\xi$, we can get all indices from shapes, e.g., $\iota(\rho\xi)$. This is true in general, no matter what the dimensionality of the array, including scalars, a 0-dimensional array.

Now building on scalar operations, and extending to scalar vector operations, scalar extension is introduced, later generalized as the outer product. This generalization is completed by adding reductions/contractions and, the inner product (defined using reduction and outer product). Mullin's research connects these four ideas mathematically in MoA; scalar operations: Matrix Multiplication (MM), Hadamard Product (HP), and the Kronecker Product (KP) using one algorithm/circuit [29]. These designs were extended to
Blocked matrix matrix multiplication by first formulating in MoA, then derived, built, and tested in C, proving its performance exceeds modern DGEMM libraries [30].

MoA defines all operations using shapes and the Psi indexing function. Hence, when operations are formulated in the MoA algebra, they can be reduced through Psi Reduction to a semantic/denotational normal form (DNF), i.e., Cartesian indices of all operations are composed. This is called the Psi Calculus. Then, one of the mapping functions, e.g., $\gamma_{row}$, transforms an index to an offset in memory, the ONF, or how to build it with knowledge of data layout that capitalizes on starts, stops. and strides. Next, formulating dimension lifting, i.e., partitioning shapes, leads to mapping data to hardware [31]. Combined with Lambda Calculus, iteration, sequence, and control [32], a Turing complete paradigm is realized to reason about array computation in general. These fundamentals allow one to define and optimize programs in any domain that use the algebra and could be enhanced subsets of any language with arrays, e.g., Fortran[2], Python[33], or Julia. That is, when there is semantic equivalence across programming languages, soft or hard, automatic linear and multi-linear transformations are easily applied, proving correctness by construction. MoA is a Universal Algebra (UA) that can optimize all domain specific languages that use arrays/tensors as their primary data structure.

### 2.3. Dense Array Processing

It is important to understand why MoA's design of Matrix Multiplication is different than other designs. Although the arithmetic complexity is the same and the arithmetic is the same, it is the fact that all arrays are accessed contiguously that give MoA's MM an edge in performance. That is, the IO (Input-Output) complexity is better than other algorithms. Even when other algorithms have a better arithmetic complexity, e.g., Strassen, the fact that contiguous access does not occur makes for poorer performance. Note, this is why the MoA MM outperforms standard libraries [30,34].

Moreover, MoA's design can treat any dimensional operation in the same way by first flattening the arrays, in row or column major ordering, followed by the use of starts, stops, and strides. The MoA MM is also very general, it is a 2-dimensional version of the inner product which can also take any two scalar operations as functional inputs, e.g., and with or. In addition, the MoA circuit not only performs the matrix multiplication, it can perform the Kronecker and Hadamard Product [35]. The only thing that differentiates the execution of multiple algorithms is that different starts, stops, and strides are fed into the circuit. This has important power and space attributes. That is, one circuit can perform four different algorithms.

Contiguous memory access easily allows the breaking up of loops, or what MoA refers to as **dimension-lifting** or reordering of the array by shape refactorization. When one loop is broken into two, it may be partitioned to more than one resource, e.g., more than one MoA circuit or vector registers versus scalar registers. This concept can easily be done for as many resources as desired in an implementation based on the size of the problem and/or the number of resources available. This concept was used to block matrices in [30,34] where each block was accessed contiguously in the L1 cache.

### 2.4. MoA Matrix Multiplication

The general matrix-matrix multiplication (GEMM) in MoA is a special case of the inner product for 2-D arrays (matrices). Define **A** as an $m \times n$ matrix, **B** as $n \times p$, and **C** as $m \times p$. Let the following notation denote the 2-d Matrix Multiplication.

$$C = A \bullet B \qquad (2)$$
In MoA notation, the shapes are:

$$ \rho \mathbf{A} = \langle m, n \rangle \quad \rho \mathbf{B} = \langle n, p \rangle \quad \rho \mathbf{C} = \langle m, p \rangle \qquad (3) $$

Next, define the valid indices of the matrices:

$$ \forall i, j, k \ni \{0 \le i < m \quad 0 \le j < p \quad 0 \le k < n\} \qquad (4) $$

$\mathbf{C} = \mathbf{A} \bullet \mathbf{B}$ is defined by the MoA psi expression, given the shapes above.

$$ \langle i \rangle \psi \mathbf{C} \equiv +_{\text{red}} (\langle i, k \rangle \psi \mathbf{A} \times (\langle k \rangle \psi \mathbf{B})) \qquad (5) $$

Equation 5 uses the general $\psi$ indexing operator and is the DNF of our GEMM expression. Once reduced to its DNF, $\gamma_{row}$ is applied thus producing the MoA Operational Normal Form (ONF) for GEMM and is given by the following "generic" program.

$$ \mathbf{C}[(i \times p) + j] := \sum_{k=0}^{n-1} \mathbf{A}[(i \times n) + k] \times \mathbf{B}[(k \times p) + j] \qquad (6) $$

Equation (6) is the generic code for a sequential program in MoA. In each of the ith rows of the resultant array C, each scalar-vector operation involving the column index j is independent of each other. The ith row of C is contiguously filled in by the summation of scalar-vector multiplications involving each matrix element at the ith row and kth column of A (the scalar) with each kth row of B (the vector) obtained by accessing the arguments contiguously. A row-wise sum reduction is then applied over the k index to yield the final answer as the ith row in C. The design is parallel at every step:

1. Parallel row operations.

2. Each n parallel row operation is performing m parallel scalar vector multiplies.

3. All n operations can be summed in log n steps (ideally).

4. Everything is accessed contiguously.

2.4.1. Example: Consider performing a 3 by 4 with a 4 by 3, 2-d matrix multiply

Consider the following example with two matrices A and B.

$$ \text{Let: } A \equiv \begin{bmatrix} 0 & 1 & 2 & 3 \\ 4 & 5 & 6 & 7 \\ 8 & 9 & 10 & 11 \end{bmatrix} \quad \text{and let } B \equiv \begin{bmatrix} 12 & 13 & 14 \\ 15 & 16 & 17 \\ 18 & 19 & 20 \\ 21 & 22 & 23 \end{bmatrix} $$

Then, the matrix multiplication is $C \equiv A \bullet B$ is performed in the following way. The first row of the result is filled in contiguously as follows while noticing how the left and right matrices are accessed contiguously.

$$ \langle 0 \rangle \psi \mathbf{C} \equiv (\sum_{n=0}^{3} (\langle 0, n \rangle \psi \mathbf{A}) \times (\langle n \rangle \psi \mathbf{B})) \qquad (7) $$

(8)

That is, the entire 0th row of C is:

$$ (0 \times \langle 12, 13, 14 \rangle) + (1 \times \langle 15, 16, 17 \rangle) + (2 \times \langle 18, 19, 20 \rangle) + (3 \times \langle 21, 22, 23 \rangle) $$
All other rows are processed as follows. Note the parallelism in all operations. The resultant C array is a 3 by 3 array.

$$ (<1> \psi C) \equiv (\sum_{n=0}^{3} (<1, n> \psi A) \times (<n> \psi B)) \quad (9) $$

$$ (<2> \psi C) \equiv (\sum_{n=0}^{3} (<2, n> \psi A) \times (<n> \psi B)) \quad (10) $$

The result is the expected 3 by 3 matrix:

$$ C \equiv \begin{bmatrix} 114 & 120 & 126 \\ 378 & 400 & 422 \\ 642 & 680 & 718 \end{bmatrix} $$

**Reduce:**

In Mullin's dissertation, operator $\pi$ was introduced and related to the mathematical notation $\Pi$. This is again, mentioned in the definition of reduction, which is now explained here to higher dimensional arrays.

$$ \forall i \text{ s.t. } 0 \le i < (\rho\xi)[0] $$

Product reduction (or any scalar operation), over the primary axis, is defined as

$$ ((\rho\xi)[0])-1 \prod_{i=0}^{((\rho\xi)[0])-1} <i> \psi\xi $$

This can be defined this way due to definitions of scalar operations on whole arrays, thus extending the example of

$$ \times \text{red} \vec{v} \equiv \prod_{i=0}^{((\tau\vec{v})-1)} \vec{v}[i] $$

This is important to relate a recursive definition to iteration, essential for an optimal implementation of theory.

**Reshape**

In Mullin's dissertation, reshape explained how to reshape an arbitrary array. The example given uses a row major ordering. However, $\gamma$ can reshape given any ordering, e.g. column major, or any other defined ordering. When implemented, it is only the shape portion of the array that is changed, leaving the linear representation of the array, e.g. row or column major ordering, unchanged. It is important to note that reshape maps an initial shape to a desired shape based on the layout of the data in memory. Ideally, contiguity of memory access (costs are modeled in later sections) is best in row major order. But that mapping can be arbitrary, as long as there is a 1-1 correspondence from one shape to another and it can be defined by a gamma $\gamma$ function.

Reshape defines dimension-lifting, a phrase defined to explain how each component of the original array's shape, after psi-reduction to it's DNF, is partitioned. This partition is based on the sizes, speeds, etc. or whatever attributes the DNF is mapped to given a set of architectural components. As each component is broken into equal parts, the original shape becomes reshaped.

### 2.4.2. Example:

Once the new shape is discovered, i.e. the lifted shape that maps with costs to the chosen architectural components. So, what was an 8 by 16 array, becomes a 2 by 4 by 2 by 8 array given that the rows are broken up over 2 processors and the columns are fed to 2 cycles of vector registers with length 8. That is, $\forall 0 \le * < i, j \le * < 8, 16 \le i, j \le \psi\xi$
becomes $\forall 0 \le * < i', i'', j', j'' > < * < 2,4,2,8 >$: $< i', i'', j', j'' > \psi_\xi$. Now casts can be calculated given the array size and datatype with how the array is mapped.

279

280

## 3. Semantics

281

In this section we use the array domains of [36] to give precise and complete denotational semantics to core MoA. In that paper we defined arrays as finite the finite Scott domains (algebraic, bounded-complete and directed-complete partial orders) of partial functions from a boxed-shaped set of indices to a given base type. The information partial order is the increasing amount of defined base-type values from an array that would have all of its values undefined. The base type is considered to be a flat domain: the undefined "bottom" value is less than all well-defined values, and those are relatively incomparable. For example value 3.0 is not comparable to 2.0 because both contain maximal information but both are greater than the undefined number. However *index* 3 is greater than *index* 2 because shape < 3 > defines a box that contains the box of shape < 2 > etc. In this order an array is greater than another one if its shape is greater or equal to that of the other one, and it is pointwise greater for its values i.e. it contains coherently more information.

282

283

284

285

286

287

288

289

290

291

292

293

### 3.1. *Denotational semantic domains*

294

The purpose of a denotational semantics is to define cleanly and completely the values and operations on arrays, including border cases i.e. null values and infinite limits. The former are well covered by MoA definitions and that is one of the theory's major advantages: it fixes the classical confusion between scalar values (arrays of empty shape), arrays of size one and base-type elements that are not arrays. Array domain theory clearly distinguishes all of them and connects pure mathematics with computing practice in a way that linear algebra ignores: array shapes are (virtually folded) parts of computer memory that has been *allocated*, undefined array elements are undefined points in that memory and defined array values are base-type elements that have been assigned. In summary: array domains are coherent with all MoA definitions and extend them to partly-defined arrays and infinite or even infinitely-dimensional arrays. Infinite values are a natural consequence of domain theory, but find no use for them in HPC or embedded systems for which we model and maximize the efficiency of *algorithms* i.e. finite computations.

295

296

297

298

299

300

301

302

303

304

305

306

307

308

Array domain theory brings credit to the following aspects of MoA definitions:

* The shape of an array is part of its value and for a given base-type *A* there is an array type (domain) *A*<sup>□</sup> that contains arrays of all shapes with elements from *A*. If shapes were part of the array's *type* there would be no notion of partial order on sub-arrays or sub-shapes which would isolate array-based programming from standard Computer Science theory where information orders are pervasive.

309

310

311

312

313

* Arrays have box shapes and not random patterns of index positions. Given a random subset of index positions where a base-type value is defined, the least array defined by those values has the shape that is the least enclosing box.

314

315

316

* Shapes containing the value 1 are allowed, e.g. shape < 3, 1, 3 > is a plane in dimensions x,z and shape < 3, 3 > is a plane in dimensions x,y. Therefore as in Pytorch [9] a dimension of length 1, allows the alignment of a tensor of fewer dimensions with another one of more dimensions before combining them.

317

318

319

320

* Shapes containing 0 values are allowed, they generate the empty ideal as input set for the content map i.e. an empty array. The shape is part of the array value (not its type) so empty < 2, 0 > is not the same as empty < 0 >: the first one is an empty plane in dimensions x,y while the other is an empty line in dimension x: can't concatenate < 0 > with < 2, 1 > but you can concatenate < 2, 0 > with < 2, 1 >.

321

322

323

324

325
In summary, the denotational semantic (array) domains are the extension to all dimensions and to all limits of the usual notion that an arrays is a map from a box-shaped integer-indexed domain to a domain of values.

## 3.2. Operational semantics

We now provide a more concrete basis semantic space for arrays and their operators. It is based on the notion that a linear ordering of the array values can serve as a relay between the logical space of indices (the denotational array domains) and the physical space of data layout in computer architectures. Some operational semantics definitions are less intuitive than their denotational counterpart. But they have the enormous advantage of spelling out the data layout. All operator definitions are based on the following.

**Definition 1. Arrays (operational definition)**

Assume a non-empty set of values A. An array $\zeta \in A^\square$ is a pair $(\zeta.\rho, \zeta.\kappa)$ where $\zeta.\rho : \mathbf{N}^*$, the **shape** is a possibly empty list of non-negative integers and $\zeta.\kappa : A^*$, the **content** is a list of values from A such that the product of the integers in the shape is the length of the content list $\Pi(\zeta.\rho) = |\zeta.\kappa|$.

Arrays of the form $([ ], [a])$ are called scalars and are in one-to-one correspondence with elements of A. When any element of $\rho$ is 0, the array is said to be empty being of the form $(\rho, [])$. An empty array is denoted $\Theta$ or $\Theta_\rho$ if its form matters (this can happen when defining a concatenation along a certain axis).

**Remark 1.** There is an infinite variety of empty arrays: any array whose shape contains a 0 is empty and is distinct from empty arrays having different shapes.

**Remark 2.** A scalar array $([], [a])$ is distinct from an array of size one with non-empty shape e.g. $([1, 1], [a])$. This is because the shape is part of the array's value and not its type. A scalar cannot be concatenated, while a non-scalar array of size 1 can be.

An array's shape defines a box-shaped set of indices and its content is given in lexicographic order w.r.t. those indices. For example $((2, 3), [(0, 0); (0, 1); (0, 2); (1, 0); (1, 1); (1, 2)])$ is the array of type $(\mathbf{N})$ shape (2, 3) whose content is the list of the indices. The complete denotational meaning of arrays, including partial order, infinite limits and information continuous functions is defined in [36]. For the purpose of understanding and using MoA for HPC applications it is sufficient to assume the usual correspondence between definition 1 and the usual "logical" meaning of arrays as maps from index boxes to values. We call the above "flat" definition the operational semantics because it resembles computer memory and is convenient for defining mappings from data and operations to hardware elements. We call the "logical" meaning the denotational semantics because it uses standard domain theory [36] and better matches the application domains for which arrays are defined.

**Definition 2. Flat indexing**

Let $\zeta$ be an array of shape $\rho = [\rho_0; \dots; \rho_{d-1}]$ and content $\kappa$. The index domain of $\zeta$ is the cartesian product $\text{dom}(\zeta) = \Pi_{i=0}^{d-1} \{0; \dots; \rho_i - 1\}$. An index of $\zeta$ is a vector $\vec{i} \in \text{dom}(\zeta)$. The flat indexing of $\zeta$ by $\vec{i} = [i_0; \dots; i_{d-1}]$ is the content element

$$ \zeta[i_0; \dots; i_{d-1}] = \kappa[i_0 \times (\rho_1 \times \dots \times \rho_{d-1}) + \dots + i_{d-1}] = \kappa[\sum_{a=0}^{d-1} (i_a \times \Pi_{b=a+1}^{d-1} \rho_b)] $$

* If $\zeta = \Theta$ then $\text{dom}(\zeta) = \emptyset$ and there exists no index $\vec{i}$ of $\zeta$. In this case $\Theta[\vec{i}]$ is not defined.

365
• Another example is $[|\xi|] = ([2,3,4], \kappa)$ and $\xi[1,2,2] = \kappa[1 \times (3 \times 4) + 2 \times 4 + 2] = \kappa[22]$. 366

The MoA algebra generalizes this indexing with the Ψ operator whose value is itself an array. The index argument of Ψ is not restricted to be of the same length as the shape of the array argument. 367 368 369

### 3.3. The array algebra 370

We now give a formal development of the core of MoA following the dependencies shown in figure 3. Our definitions take the form of an operational semantics $[|\_|]$ whose input is an expression of the MoA algebra, and whose value is an operational semantics that uses the flat definition. 371 372 373 374

**Definition 3.** MoA algebra 375

MoA is a one-sorted concrete algebraic structure: its base-type of elements is abstract but the array operators themselves are not axiomatized but rather defined. All constants and variables are of array type written $X^\square$ where X is either $\mathbf{N}$ or a pure set$^2$ of base elements A. All other function symbols represent transformers, so their types are of the form $A^\square \times \dots \times A^\square \to A^\square$. Its constants are the empty arrays $\Theta, \Theta_\rho$ and empty vector $<>= \Theta_{<0}>$. There are scalar variables $\sigma, a, \dots, z$ and general array variables $\xi, \xi^n, \xi_\rho$ etc. The function symbols of arity 1 or more are those of figure 3. Symbol $\rho$ is overloaded many times. As a variables it denotes a shape vector (one-dimensional array), or a shape list in the operational semantics. As a unary function symbol it denotes the operation that takes an array and returns this shape vector (or shape list meaning). As a binary function symbol it is called reshape and denotes an operation that takes a new shape, an array and produces a new array modified by this new shape. The operational semantics $[|\_|]$ takes a syntactic element of MoA and maps it to an array or array transformation in the flat representation. It is defined below for each operator. The algebraic laws of MoA are equivalence theorems proven from the denotational semantics. 376 377 378 379 380 381 382 383 384 385 386 387 388 389

**Remark 3.** Array constructors and destructors are implicit. In practice this means that a software tool or API implementing MoA has predefined array constants, input/output operations that import/export arrays to or from other data types and channels. The action of such operations is not part of MoA theory as it is context dependent and not necessary for defining parallel algorithms. 390 391 392 393

<footnote n="2">² For applications to linear algebra it will be assumed that $A$ is a field $F$ or at least a semi-ring. As we will see, the only elementary operator that requires an associative law is the reduction.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 3. MoA definition dependencies</figcaption></figure>

**Definition 4.** *The shape operator*

The shape operator $\rho_-$ is unary, of type $A^\square \to N^\square$ and its operational meaning is given by

$$ \frac{[|\xi|] = (\xi.\rho', \xi.\kappa) \quad |\xi.\rho'| = n}{[|\rho(\xi)|] = ([n], \xi.\rho')} $$

In other words $\rho$ returns a vector whose content is the shape of its argument. For example, if $[|\xi|] = ([2;2], [a;b;c;d])$ then $|[2;2]| = 2$, therefore $[|\rho(\xi)|] = ([2], \xi.\rho') = ([2], [2;2])$ which is usually written $< 2;2 >$ in vector (1-dim arrays) MoA notation.

Operator $\delta$ is called "dimensionality" to avoid confusing with dimension sizes. It counts the number of dimensions i.e. length of the shape of an array.

**Definition 5.** *The dimensionality operator*

The dimensionality operator $\delta_-$ is unary, of type $A^\square \to N^\square$ and has meaning

$$ \frac{[|\xi|] = (\xi.\rho, \xi.\kappa) \quad |\xi.\rho| = m}{[|\delta(\xi)|] = ([ ], [m])} $$

So $\delta$ returns a scalar whose value is the length of the shape, or number of dimensions of its argument. For example, if $[|\xi|] = ([2;2], [a;b;c;d])$ then $|\xi.\rho| = 2$ therefore $[|\delta(\xi)|] = ([ ], [2])$ which, in MoA notation is the scalar $\sigma = 2$.

**Definition 6.** *The (total) size operator*

The size operator $\tau_-$ is unary, of type $A^\square \to N^\square$ and has meaning

$$ \frac{[|\xi|] = (\xi.\rho, \xi.\kappa) \quad \Pi(\xi.\rho) = t}{[|\tau(\xi)|] = ([ ], [t])} $$
So $\tau$ returns a scalar whose value is the product of the shape elements, or size of the index domain, of its argument.

For example, if $[|\xi|] = ([2; 2], [a; b; c; d])$ then $\Pi(\xi.\rho) = 4$, therefore $[|\tau(\xi)|] = ([], [4])$ which is the scalar $\sigma = 4$.

**Remark 4.** By definition of arrays, the total size value is the length of the content list.

We give here the well-known definition of concatenation and map on finite lists.

**Definition 7.** List concatenation.

List concatenation _ ++_ has type $A^* \times A^* \rightarrow A^*$, neutral element $_-_$, and is defined by
$$[x_0; \dots; x_{a-1}] ++ [y_0; \dots; y_{b-1}] = [x_0; \dots; x_{a-1}; y_0; \dots; y_{b-1}].$$

List concatenation is associative so whenever $A$ is non-empty, $(A^*, ++, [\_])$ is a monoid,
and list-length is a homomorphism $A^* \rightarrow (\mathbf{N}, +, 0)$.

List concatenation has an analog for vectors (1-dimensional arrays).

**Definition 8.** Vector concatenation.

Vector concatenation _ ++_ : $A^\square \times A^\square \rightarrow A^\square$ is defined when both arguments are vectors.

$$ \frac{[|\vec{u}|] = ([m], \vec{u}.\kappa) \quad [|\vec{v}|] = ([n], \vec{v}.\kappa)}{[|\vec{u}++\vec{v}|] = ([m+n], \vec{u}.\kappa++\vec{v}.\kappa)} $$

In other words, vector-concatenation concatenates contents and adds lengths to produce a vector.

Vector concatenation is associative with the empty vector as left/right neutral element.

**Definition 9.** List map.

List map _* has type $(A \rightarrow B) \rightarrow (A^* \rightarrow B^*)$ and is defined by $f^*[\quad] = [\quad]$ and $f^*([a] ++l) = [f(a)] ++l$. $f^*$ is an endomorphism on the list monoid.

The semantics of the take (array) operation uses the following list prefix-suffix operation that will be used for taking content.

**Definition 10.** List take, or "prefix-suffix".

We define the list-take operation _↑_ : $\mathbb{Z} \times A^* \rightarrow A^*$. Assume $\kappa$ is a list of length $n$ and $s : \mathbb{Z}$ such that $|s| \le n$. When $s$ equals 0, then $s \uparrow \kappa = []$. When $s$ is positive, $s \uparrow \kappa$ equals $[\kappa[0]; \dots; \kappa[s-1]]$ but if $s$ is negative, then it equals $[\kappa[n+s]; \dots; \kappa[n-1]]$.

When the right argument of take is an array $\xi$, operation $2 \uparrow \xi$ is well-defined if $\xi$'s first dimension $^3$ is at least of length 2. Let us suppose that $\xi$'s shape list is $[3; 4; 5]$. Then $\xi$ is abstractly equivalent to a 3-list of arrays of shape $[4; 5]$. Then producing the value of $2 \uparrow \xi$ is equivalent to taking the first two such sub-arrays to produce an array of shape list $[2; 4; 5]$. The operational semantics of the MoA take operation is defined below. It may seem less intuitive or complex than the above explanation using nested array. But it has enormous logical and computational advantages: (a) absence of nesting removes the use of pointers, shape verifications and multiple types (b) moreover an implementation can translate the flat semantics directly and statically into hardware mapping information.

**Definition 11.** The take operator

The take operator _↑_ is binary, of type $\mathbb{Z}^\square \times A^\square \rightarrow A^\square$. It is defined when its left argument is a

<footnote>$^3$ The technical definition of array take (definition 11) allows for a scalar right argument without "first" dimension.</footnote>
scalar integer, its right argument is an array of dimensionality at least 1, and it has the following meaning.

$$
\begin{array}{l}
\sigma = (\llbracket \cdot \rrbracket, \llbracket s \rrbracket); \llbracket |\xi| \rrbracket = (\xi.\rho, \xi.\kappa); |\xi.\rho| \ge 1; \\
\qquad 1 \uparrow \xi.\rho = [n]; |s| \le n; (1 \downarrow \xi.\rho) = \rho'; \Pi(\rho') = m \\
\qquad \llbracket |\sigma \uparrow \xi| \rrbracket = (\llbracket |s| \rrbracket + +\rho', (s \times m) \uparrow (\xi.\kappa))
\end{array}
$$

In other words a take expression is well defined if the left argument σ is an integer scalar (s) whose absolute value is no greater than the product of the first dimension of the array. Its result is an array whose content is a list prefix (or suffix if s < 0) of the argument's content and the shape begins with s and ends with the tail of the original shape.

* For example if $\sigma = 2$ is the left argument of take, $\zeta$ the right argument is an array ($[3;4;5]$, $[\kappa[0];...;\kappa[59]]$). Then the definition's parameters are $s = 2$, $n = \Pi([3]) = 3$, $s \le n$ is verified, $\rho' = [4;5]$, $m = \Pi([4;5]) = 20$. Therefore the result $[|\sigma \uparrow \zeta|]$ has shape $[|s|] ++\rho' = [2] ++[4;5] = [2;4;5]$, its content is $(s \times m) \uparrow (\zeta.\kappa) = (2 \times 20) \uparrow (\zeta.\kappa) = 40 \uparrow (\zeta.\kappa) = [\kappa[0];...;\kappa[39]]$, the content of the first two "virtual sub-arrays" of shape $[4;5]$.

* If the left scalar argument had been $\sigma = -2$, the shape of the result would be the same, and the content would have been $(-40) \uparrow (\xi.\kappa)$ i.e. the last two "virtual sub-arrays" of shape $[4;5]$.

* If the left scalar argument has been $\sigma = 0$, the result's shape would be $[0;4;5]$ and the content would be $0 \uparrow (\xi.\kappa) = []$ and that defines $\Theta_{[0;4;5]}$.

The semantics of the drop (array) operation is symmetrically based on the following list suffix-prefix operation that will be used for taking content.

**Definition 12.** *List drop, or "suffix-prefix".*

Assume $\kappa$ is a list of length $n$ and $s: \mathbb{Z}$ such that $|s| \le n$.

When $s \ge 0$, $s \downarrow \kappa = [\kappa[s]; \dots; \kappa[n-1]]$, while if $s < 0$, then $s \downarrow \kappa = [\kappa[0]; \dots; \kappa[n-1+s]]$.

The meaning of the drop operator is a mirror image of that of the take operator.

**Definition 13.** *The drop operator*

The drop operator $\_\downarrow$ is binary, of type $\mathbb{Z}^\square \times \mathbb{A}^\square \to \mathbb{A}^\square$. It is defined when its left argument is a scalar, its right argument is an array of dimensionality 1 or more, and it has the following meaning.

$$
\begin{array}{l}
\sigma = (\llbracket \cdot \rrbracket, \llbracket s \rrbracket); \llbracket |\xi| \rrbracket = (\xi.\rho, \xi.\kappa); |\xi.\rho| \ge 1; \\
\qquad 1 \uparrow \xi.\rho = [n]; |s| \le n; (1 \downarrow \xi.\rho) = \rho'; \Pi(\rho') = m \\
\qquad \llbracket |\sigma \downarrow \xi| \rrbracket = (\llbracket n - |s| \rrbracket + +\rho', (s \times m) \downarrow (\xi.\kappa))
\end{array}
$$

In other words a drop expression is well defined if the left argument $\sigma$ is an integer scalar (s) whose absolute value is no greater than the first dimension of the array argument. Its result is an array whose content is a list suffix (or prefix if $s < 0$) of the argument's content and the shape begins with $n-s$ and ends with the tail of the original shape.

The next definition is an index constructor $\iota$. We begin with its list version and use it to define the multi-dimensional MoA meaning in a manner similar to how list take was used to define array take.

**Definition 14.** *The list index generator*

Assume $n: \mathbb{N}$, then $\iota: \mathbb{N} \to \mathbb{N}^*$ is defined as follows: $\iota 0 = []$ and $\iota(n+1) = \iota(n) ++[n]$.

For example $\iota(3) = \iota(2) ++[2] = \iota(1) ++[1] ++[2] = \dots = [0;1;2]$.

The MoA index generator produces an array of indices.
**Definition 15.** The index generator
The index generator $\iota$ is unary of type $\mathbf{N}^\square \to \mathbf{N}^\square$. It is defined when its argument is a scalar or a vector (of integers), and has the following meaning.

$$ \frac{[|\sigma|] = ([], [n])}{[|\iota(\sigma)|] = ([n], \iota(n))} \quad \frac{[|\xi|] = ([n], v) (\prod v) = 0}{[|\iota(\xi)|] = ([0], [])} \quad \frac{[|\tilde{\xi}|] = ([n], v) (\prod v) \ge 1}{[|\iota(\tilde{\xi})|] = (v ++[n], \kappa)} = \kappa $$

where $\iota^*(\rho) \in (\mathbf{N}^*)^*$ is the mapping of the list-index generator on the list of integers $\rho$, $\Pi(\iota^*(\rho))$ is the cartesian product of those index lists, and finally $++\Pi(\iota^*(\rho))$ is the concatenation of all those indexing lists of equal length. The condition $(\prod v) \ge 1$ ensures none of the elements of $v$ is zero.

* For example, the rule for scalar arguments applies to $\iota(\sigma = 3)$. The semantics of $\sigma$ is $([], [3])$ so $[|\iota(\sigma)|] = ([3], \iota(3))$ which is the vector $< 0; 1; 2 >$ in MoA notation. In summary $\iota(\sigma = 0) = < >$, $\iota(\sigma = 1) = < 0 >$, $\iota(\sigma = 2) = < 0; 1 >$, etc.

* If the argument vector contains a 0, then the second rule applies, $\prod v = 0$ and $\iota(\xi)$ is then the empty vector $<>$.

* If the argument is a vector, the third rule applies. Let's apply $\iota$ to the empty vector $\Theta_{[0]} = < >$ whose meaning is $([0], [])$. In the left part of the inference rule's denominator, $([n], v) = ([0], [])$ so $n = 0$ and $v = []$. The right part is then $\kappa = ++\Pi(\iota^*(v)) = ++\Pi(\iota^*([])) = ++\Pi(’) = ++([’]) = [’]$. The denominator then yields $\iota(<>) = (v ++[n], \kappa) = ([] ++[0], []) = ([0], [])$ which is again the empty vector. So $\iota(<>) = < >$.

* As an example of the most general case (without zeros) assume $\xi = < 3; 4 >$ and let us compute $\iota(\xi)$. The numerator of the second rule produces $n = 2$, $v = [3; 4]$ and $\kappa = ++\Pi([0; 1; 2]; [0; 1; 2; 3]) = [0; 0; 0; 1; 0; 2; 0, 3; 1; 0; \dots; 2; 3]$ where every block of length 2 is a possible index for an array of shape $[3; 4]$. The meaning of $\iota(< 3; 4>)$ is thus $(v ++[n], \kappa) = ([3; 4; 2], [0; 0; 0; 1; 0; 2; 0, 3; 1; 0; \dots; 2; 3])$. Intuitively it corresponds to a 3 array of indices of length 2.

Next we consider indexing. Flat indexing (definition 3.3) $\xi[i_0; \dots; i_{d-1}]$ has two defects. 1° First it only applies to index lists that are exactly $d$ long, where $d$ is the dimensionality of $\xi$. We generalize this by (a) giving the result an array type, not the base type, and (b) allowing shorter index vectors. 2° Another, less obvious defect is that its definition uses enumeration lists $[i_0; \dots; i_{d-1}]$ specific to a certain dimensionality $d$. We replace this semi-formal description by a binary indexing operator $\Psi$ whose arguments are arrays. Its left argument provides indices, possibly produced by $\iota$, and the right one provides content. General array indexing is built from the following auxiliary function.

**Definition 16.** The one-integer indexing operator.
One-integer indexing is binary, of type $\_[-] : A^\square \times \mathbf{N} \to A^\square$ and is well-defined when the array argument is non-empty, the integer argument is no greater than the length of the first dimension of its array argument.

$$ \frac{\xi \neq \Theta \quad [[\xi]] = (\xi.\rho, \xi.\kappa) \quad 1 \uparrow \xi.\rho = [n] \quad 0 \le i < n \quad \Pi(1 \downarrow \xi.\rho) = m}{[[\xi[i]]] = (1 \downarrow \xi.\rho, m \uparrow ((i * m) \downarrow (\xi.\kappa)))} $$

* Scalar array arguments are not allowed, they are excluded by the $1 \uparrow \xi.\rho = [n]$ hypothesis in the rule. Indeed for a scalar $\xi = (\sigma = s), \xi.\rho = []$ and so $1 \uparrow \xi.\rho$ is undefined. As a result, the variable $n$ is also undefined etc.

* If the array is 1-dimensional i.e. a vector, then by hypothesis $\xi \neq \Theta$ so it is a non-empty vector $\xi = < x_0; \dots; x_{n-1} >$. Then $\xi.\rho = [n]$ so $1 \uparrow (\xi.\rho) = [n]$, so the possible indices are $i = 0; \dots; i = n-1$. Also, $\Pi(1 \downarrow \xi.\rho) = \Pi(’) = 1 = m$, and the denominator of the inference rule yields an array of empty shape $1 \downarrow \xi.\rho = []$ and content $m \uparrow ((i * m) \downarrow (\xi.\kappa))) = 519$
1 ↑ ((i * 1) ↓ (ξ.κ))) = 1 ↑ (i ↓ ([x₀; ... ; xₙ₋₁])) = xᵢ. So in this case [[ξ[i]]| is the scalar (σ = xᵢ). 520

* Assume $\xi.\rho = [2;256;2]$, $\xi.\kappa = [x_0;...;x_{1023}]$ and let us compute $\xi[i]$ for $i = 0$. In the numerator of the inference rule we have $1 \uparrow \xi.\rho = [2]$ so $n = 2$, and $0 \le i < n$ is verified. Also $\Pi (1 \downarrow \xi.\rho) = \Pi [256;2] = 512 = m$. The denominator of the rule yields for value of $\xi[i]$ the following shape $1 \downarrow \xi.\rho = [256;2]$ and content $m \uparrow ((i * m) \downarrow (\xi.\kappa)) = 512 \uparrow ((0 * 512) \downarrow [x_0;...;x_{1023}]) = 512 \uparrow ([x_0;...;x_{1023}]) = [x_0;...;x_{511}]$) = which is the first (0th) sub-array of shape [256;2] in $\xi$. 521
522
523
524
525
526
527

* With the same array $\xi$, if the integer-indexing was $i = 1$, the numerator would give $n = 2$, and $0 \le i < n$ still verified. $m = 512$ would not change, and the denominator would yield the same shape of [256;2] but content $512 \uparrow ((1 * 512) \downarrow (\xi.\kappa)) = 512 \uparrow ([x_{512};...;x_{1023}]) = [x_{512};...;x_{1023}]$ i.e. the second (rank 1) sub-array of shape [256;2] in $\xi$. 528
529
530
531

**Definition 17.** The indexing operator
The indexing operator _Ψ_ is binary of type N^□ × A^□ → A^□. It is defined when its left argument
i̇ is a vector no longer than the shape of its right argument ξ. If the left argument is empty, the
operation is the identity.

$$
\frac{\overline{[i]}}{\overline{[i\Psi\xi]}} = \xi
$$

$$
\frac{\overline{[i]}}{\overline{[i\Psi\xi]}} = \xi
$$

In other words (< > Ψ_) is the identity function on arrays.
If the left argument is of length 1 or more, indexing produces a sub-array of ξ whose content is
defined recursively by one-integer indexing.

536
537
538

$$
\frac{\xi \neq \Theta \quad 1 \uparrow \vec{i} = <i> \quad 1 \downarrow \vec{i} = \vec{j} \quad \xi[i] = \xi'}{[\vec{i}\Psi\xi] = \vec{j}\Psi\xi'}
$$

* An empty array argument to Ψ requires a special rule because one-integer indexing Θ[i] is not defined. Operation i̇Ψ is the identity in this case.
539
540

* Suppose the indexing vector is of length 1 as in < i > Ψξ. Then 1 ↑< i > =< i >, 1 ↓ i̇ = < >, ξ[i] = ξ' and the inference rule gives [[ < i > Ψξ]] = < > Ψ(ξ[i]) which equals ξ[i] by the first inference rule. So < i > Ψξ is the one-integer indexing ξ[i]. 541
542
543

* The general case for i̇Ψξ is when i̇ = < i > ++ j̇. The second inference rule is then i̇Ψξ = j̇Ψ(ξ[i]). 544
545

* For example let $\xi$ be a 3-dimensional array of shape $[5;5;2]$, content $[x_0;x_1;...;x_{49}]$ and let us compute $< 1;0 > \Psi\xi$. The above rephrasing of rule 2 produces $< 1;0 > \Psi\xi = < 0 > \Psi(\xi[1]) = < > \Psi(\xi[1])[0]$ which, by the first rule is equal to $(\xi[1])[0]$. Two applications of definition 16 then give a result of shape [2] and content
  2 ↑ ((0 * 2) ↓ (10 ↑ ((1 * 10) ↓ [x_0;...;x_{49}])))) =
  2 ↑ ((10 ↑ (10 ↓ [x_0;...;x_{49}])))) =
  2 ↑ ((10 ↑ [x_{10};...;x_{49}])) = 2 ↑ ([x_{10};...;x_{19}]) = [x_{10};x_{11}]
  i.e. the first sub-array of shape [2] in the second sub-array of shape [5;2] of the initial array.
546
547
548
549
550
551
552
553

A direct consequence of definition 17 is that it distributes over index concatenation.
554

**Proposition 1.** If $\delta(\vec{u}) = \delta(\vec{v}) = 1$ and $\vec{u}, \vec{v} : \mathbf{N}^{\square}$ and $(\tau\vec{u}) + (\tau\vec{v}) \leq \delta\xi$ then
555

$$
(\vec{u} + \vec{v}) \Psi \xi = \vec{v} \Psi (\vec{u} \Psi \xi)
$$

In other words, when all indexings are well defined (ū++v̄)Ψ_ = (v̄Ψ_) ◦ (ūΨ_). Whenū and v̄
are both of length one, this property can be written <i,j> ΨA = <j> Ψ(<i> > ΨA).
556
557
Reshaping transforms an array's shape while retaining its content. It lies at the core of our methodology for mapping MoA expressions to parallel architectures.

**Definition 21.** The reshape operator

The reshape operator `_ρ_` is binary, of type $N^□ \times A^□ \rightarrow A^□$ and is well-defined when its left argument is a vector.

$$ \frac{\zeta \neq \Theta \ \mathrm{rav}(\zeta)[0] = x}{[\ll \hat{\rho}\zeta] = ([], x)} $$

$$ \frac{\zeta \neq \Theta \quad [|\zeta|] = (\zeta.\rho, \zeta.\kappa) \quad \delta(\rho') = (\sigma = 1) \quad \Pi(\zeta, \rho) = n \quad \Pi(\rho') = n' \\ \forall i. (0 \le i < n'). x_i = \zeta.\kappa[i \bmod n] \quad \kappa' = [x_0; \dots; x_{n'-1}] \\ [|\rho' \hat{\rho}\zeta|] = (\rho', \kappa')} $$

In other words, the result of reshaping by an empty vector produces a scalar whose value is the first element in the content of the array argument. This non-intuitive property, is a necessary consequence of the semantics of MoA and the property that reshaping by $ρ'$ produces an array of shape $ρ'$. That property also holds for the general case, and the content of the reshaped array is either a prefix or an extension of the original content by cyclic repetition.

The next operator we define is core to all linear-algebra applications of arrays, namely reduction of the array for a given binary operations. Reduction depends on the following auxiliary definition of scalar extension, namely pointwise extension to arrays of a given scalar operation.

**Definition 22.** Scalar (pointwise) extension.

The extension of a given operation on type A, applies to arrays of $A^□$ and operates on its elements pointwise.

$$ \frac{\mathrm{uop}_- : A \rightarrow A \quad [|\xi|] = (\xi.\rho, \xi.\kappa) \quad \xi.\kappa = [x_0; \dots; x_{n-1}] \quad \kappa' = [\mathrm{uop}(x_0); \dots; \mathrm{uop}(x_{n-1})] \\ [|\mathrm{uop}(\xi)|] = (\xi.\rho, \kappa')} $$

$$ \frac{\mathrm{_op}_- : A \times A \rightarrow A \quad [|\xi|] = (\xi.\rho, \xi.\kappa) \quad [|\xi'|] = (\xi'.\rho, \xi'.\kappa) \quad \xi.\rho = \xi'.\rho \\ \xi.\kappa = [x_0; \dots; x_{n-1}] \quad \xi.\kappa = [y_0; \dots; y_{n-1}] \quad \kappa'' = [(x_0)\mathrm{op}(y_0); \dots; (x_{n-1})\mathrm{op}(y_{n-1})] \\ [|\xi)\mathrm{op}(\xi')|] = (\xi.\rho, \kappa'')} $$

Both definitions apply to empty arrays and produce empty arrays: $\mathrm{uop}(\Theta_\rho) = (\Theta_\rho)$ and $(\Theta_\rho)\mathrm{op}(\Theta_\rho) = (\Theta_\rho)$.

**Definition 23.** The reduction operator

Suppose an associative binary operation $\oplus$ on type A, having neutral element $id(\oplus)$ i.e. a monoid $(A, \oplus, id(\oplus))$. The reduction operator for arrays of $A^□$, is of type

$$ \mathrm{red}_{\oplus}: (A \times A \rightarrow A) \times A^\square \rightarrow A^\square $$

and is defined when its argument has dimensionality 1 or more.

$$ \overline{\mathrm{red}_{\oplus}(\Theta_{\rho})} = (1 \downarrow \rho)\hat{\rho}(\mathrm{id}(\oplus)) $$

$$ \frac{\tau(\rho) \ge 1 \quad \zeta \neq \Theta}{\mathrm{red}_{\oplus}(\zeta) = (<0> \Psi\zeta) \oplus (\mathrm{red}_{\oplus}(1 \downarrow \zeta))} $$

Notice: the central $\mathrm{red}_{\oplus}$ operation in the last rule's denominator is a scalar extension to arrays.
Neutral elements are needed to apply the first rule. For example if $(A, \oplus, \mathrm{id}(\oplus))$ is $(\mathbb{N} \cup \{\infty\}, \max, 0)$ then $\mathrm{id}(\oplus) = 0$. In $(\mathbb{N} \cup \{\infty\}, \min, \infty)$ the neutral element is $\mathrm{id}(\oplus) = \infty$. In $(\mathbb{R}, \times, 1)$ we have $\mathrm{id}(\oplus) = \Pi(\emptyset) = 1$, In $(\mathbb{R}, +, 0)$ we have $\mathrm{id}(\oplus) = \Sigma(\emptyset) = 0$, etc. But the classical theory of monoids, semi-rings (for graph algorithms) and rings/fields (for general linear algebra) does not handle shapes in a systematic and multi-dimensional manner as MoA does. This is why the first rule defines the reduction of an empty as a very specific reshaping of the scalar $\mathrm{id}(\oplus)$. In fact the shape of the result of reduction is always 1 dropped from the shape of the argument.

**Proposition 3.** $\rho(\mathrm{red}_{\oplus}(\xi)) = 1 \downarrow (\rho(\xi))$.

We now illustrate definition 23 by an exhaustive set of examples.

The first rule applies to any empty array $\Theta_{\rho}$. By definition of emptiness, that array's shape $\rho$ cannot be empty, it must contain one 0 value. Suppose that $\rho = [0] ++\rho'$ where $\rho'$ could be empty. Then $\mathrm{red}_{\oplus}(\Theta_{\rho}) = (1 \downarrow \rho)\hat{\rho}(\mathrm{id}(\oplus)) = \rho'\hat{\rho}(\mathrm{id}(\oplus))$. If $\rho' = []$, the result is the scalar $\mathrm{id}(\oplus)$. If $\rho'$ is not empty, for example $\rho' = [4, 5]$ then the result is a $4 \times 5$ array filled with the neutral element $\mathrm{id}(\oplus)$ (consequence of definition 21). An even less intuitive but coherent special case would be $\rho' = [4, 0]$ i.e. the result is again an empty array $\Theta_{[4,0]}$.

Consider now the second rule of definition 23. It is a head-recursive definition of reduction over the shape of any array argument. Assume for example a vector argument. It states that for a non-empty vector, $\vec{v}$, the reduction is the abstract sum of its first element with the reduction of its tail (remainder $1 \downarrow \vec{v}$). For example in $(\mathbb{R}, +, 0)$, if $\vec{v} = \langle x_0, x_1 \rangle$ then $\mathrm{red}_+(\vec{v}) = x_0 + \mathrm{red}_+(1 \downarrow \vec{v}) = x_0 + \mathrm{red}_+(\langle x_1 \rangle) = x_0 + (x_1 + \mathrm{red}_+(\langle \rangle)) = x_0 + (x_1 + \mathrm{red}_+(\Theta_{[0]})) = x_0 + (x_1 + \mathrm{id}(+)) = x_0 + (x_1 + 0)$. In general $\mathrm{red}_+(\langle x_0, \dots, x_{n-1} \rangle) = \sum_{i=0}^{n-1} x_i$. In summary, if $\zeta = \vec{v}$ is a non-empty vector, then the denominator of the third rule reads

$$ \mathrm{red}_{\oplus}(\vec{v}) = (< 0 > \Psi \vec{v}) \oplus (\mathrm{red}_{\oplus}(1 \downarrow \vec{v})) = \vec{v}[0] \oplus (\mathrm{red}_{\oplus}(1 \downarrow \vec{v})). $$

Take now a 2-dimensional example. Let $\zeta_{1,5}$ be $[1,5]\hat{\rho}(\sigma = 1)$, the array of shape $[1,5]$ whose content list is $[1, 1, 1, 1, 1]$.

$$ \begin{array}{l} < 0 > \Psi(\zeta_{1,5}) = \text{the first (only) row i.e. the vector } \langle 1, 1, 1, 1, 1 \rangle \\ 1 \downarrow (\zeta_{1,5}) = \Theta_{[0,5]} \text{ the empty array of shape } [0,5] \text{ so} \\ \mathrm{red}_+(1 \downarrow (\zeta_{1,5})) = \mathrm{red}_+(\Theta_{[0,5]}) = [5]\hat{\rho}(\sigma = 0) = \langle 0, 0, 0, 0, 0 \rangle \text{ so finally} \\ \mathrm{red}_+(\zeta_{1,5}) = (< 0 > \Psi(\zeta_{1,5})) + (\mathrm{red}_+(1 \downarrow (\zeta_{1,5}))) = \\ < 1, 1, 1, 1, 1 > + < 0, 0, 0, 0, 0 > = \langle 1, 1, 1, 1, 1 \rangle. \end{array} $$

More interestingly, let $\zeta_{2,5}$ be $[2,5]\hat{\rho}(\sigma = 1)$, of shape $[2,5]$ filled with values 1, then

$$ \begin{array}{l} < 0 > \Psi(\zeta_{2,5}) \text{ is the first row i.e. the vector } \langle 1, 1, 1, 1, 1 \rangle \\ 1 \downarrow (\zeta_{2,5}) = \zeta_{1,5} \text{ so the reduction is } \mathrm{red}_+(\zeta_{2,5}) = \langle 1, 1, 1, 1, 1 \rangle + \mathrm{red}_+(\zeta_{1,5}) = \\ < 1, 1, 1, 1, 1 > + < 1, 1, 1, 1, 1 > = \langle 2, 2, 2, 2, 2 \rangle. \end{array} $$

In summary, $+-$reduction of $\zeta_{2,5}$ amounts to a $+-$reduction of each line of the matrix. Similarly, $+-$reduction of a 3-dimensional array would reduce lines orthogonal to the first-dimension's planes, etc. To compute a single scalar from the whole content of the array one applies reduction as many times as there are dimensions (associativity guarantees the unique value).

**Remark 5.** Our definition 23 requires that the binary operation be associative. This is logically independent from the formal definition of MoA reduction as written above, and hence appears unnecessary. Yet there are two compelling reasons to avoid non-associative reductions (like arbitrary
accumulators, arbitrary loop bodies in programming languages, or inverses of associative operations like subtraction):

1. Applications of arrays to numerical algebra where the base set $A$ is a semi-ring (or often a field) impose binary operations that are both associative.

2. Parallel algorithms are all based on either reduction, scan or both operations [38]. Even so-called "embarassingly parallel" applications require reductions for their final statistics.

The reason for the deep relationship between parallel algorithms and the associative law is the following: reduction with an associative binary operation can be parenthesized in any order, and this allows for parallel execution of sub-trees in the reduction's syntax. For example $\text{red}_+(< x_0, \dots, x_{1023} >) = \sum_{p=0}^{15} (\sum_{i=0}^{63} x_{i+16 \times p})$, and this syntax gives the structure of a 16-wide parallel algorithm for the sum: 16 sequential sums of length 64 etc. Its binary form is given by the following property in one dimension.

**Proposition 4.** For any $k < \delta(\vec{v})$, $\text{red}_{\oplus}(\vec{v}) = (\text{red}_{\oplus}(k \uparrow \vec{v})) \oplus (\text{red}_{\oplus}(k \downarrow \vec{v}))$.

Recursive application of this property (for example in binary form when $k = \delta(\vec{v})/2$) naturally produces a parallel algorithm for reduction.

**Example 1.** The only core MoA operators that we omit from this presentation are the symmetries, including the general binary transpose. It applies an arbitrary permutation of the array index "axes". Assume $T: A^{\square} \to A^{\square}$ is a transpose that exchanges the i-th and j-th axes of an array. Assume also that $\zeta$ has 5 dimensions numbered (0,...4). Then reduction of $\zeta$ along its dimension 3 can be obtained by $T' \circ \text{red} \circ T_{0,3}$ where $T'$ is the appropriate transposition to return the remaining dimensions to their initial order. We have thus shown how reduction can apply to any axis.

# 4. Cost models, mapping arrays to processors

## 4.1. MoA expressions as high-level programs

The MoA algebra gives rise to a large number of equivalences among expressions on arrays. We have only given a glimpse of them but the reader can imagine their combinatorial space as all possible combinations of transformers and the only reducer (red) with axis permutations (transpositions) that we have not defined here for lack of space. That space contains the algebraic definition of all known (and probably all possible) operations in linear algebra: matrix multiplication, outer products, Hadamard product, FFT, etc.

As a result one can use MoA to give a very high-level definition of those algorithms. Once that is done, it is possible to transform the expression into the MoA operational normal form (ONF, subsection 2.4). In the remainder we will explain our methodology and thus introduce the key concept of cost function:

1. represent the specification of an algorithm as an MoA expression

2. reduce it to normal form and derive a corresponding concrete algorithm from it (i.e. an execution schedule)

3. compute a *cost* for that concrete algorithm that is an abstraction of execution time (or space, or energy)

4. and then if needed algebraically transform the normal form to seek a lesser cost for this equivalent concrete algorithm.

Let us now use this methodology and the generic example of matrix multiplication.

We have seen how matrix multiplication $C = A \bullet B$ for arrays of shapes $\rho A = \langle m, n \rangle$, $\rho B = \langle n, p \rangle$, $\rho C = \langle m, p \rangle$ and indices $0 \le i < m, 0 \le j < p, 0 \le k < n$ is defined by an MoA expression that we call the **denotational normal form** (DNF) because it is pure MoA and its indexing operations are at the leaves of the expression. Function $\gamma$ is then applied
architectures. Related work [18] calls it a performance diagnostic method and estimates execution time by the sum of one term per memory-hierarchy level caches L1, L2, L3. Each such term is a ratio derived from published information or measurement of assembly operations width and bandwidth for that part of the hardware. A later work [39] refines the same idea to account for time variations of the processor clock frequency, and application to a power-energy formula. Here power is a polynomial formula that is quadratic over f (frequency) and depends linearly on execution time.

Such work is of great importance because it gradually covers all the practical parame-
ters that determine the execution time of a given program on a given architecture. From
this we draw several observations:

a techniques such as the roofline model are necessary for realistic diagnostic of execution time, space or energy.

b but they only refer to measurements, not to source code like our MoA ONF or fragment-C;

c they depend on the architecture's parameters: transfer speeds relative to arithmetic operation speeds, number of computation units, hierarchies such as in caches etc.

d so they can require a great variety of execution time formulas, yet those formulas follow a certain structure. For example the refinement cited above generalizes clock frequency to be variable and not constant. In this formula, execution time may be an integral or weighted sum rather than a simple count of the arithmetic operations and transfer times.

Those observations lead us to propose an axiomatic notion of execution cost, directly inspired by operations research and graph theory, but also connected to our structured notion of a source program, namely MoA ONF.

It is folklore in the HPC community that execution time prediction is not a science but an art. We wish to demonstrate that it is rather skilfull engineering that can grow into an experimental science if an appropriate theory is provided.

4.3. *The cost of an MoA program: theory to explain the experiments*

The fragment-C code of example 2 does not specify any parallel execution. A compiler must produce parallel (vector, multithread etc.) instructions from such loops and that is usually hidden to the developer who is left guessing what source code and *directives* (requests to the compiler that have no fixed semantics) must be produced to obtain the best performance.

We introduce here a simple and general mathematical basis for this compilation process. First, we make the compiler's choice for parallel execution explicit. Loops are either named seq for sequential execution or par for parallel execution. Consider for example the inner loop of our previous example, put in explicit sequential mode and assume that *p* = 4.

**Example 3.**

```c
{ seq (j=0;j<4;j++)
   { C[(4*i)+j]+=A[(i*n)+k]*B[(4*k)+j]; }
}
```

The execution graph is then a linear chain of such blocks (figure 4)<fnref n="4" />

<footnote n="4">⁴ The inner loop is proportional to the size of the original MoA expression so it is reasonable and customary to treat it as an atomic block of sequential instructions.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 4. Execution graph of the seq(4) algorithm variant, seq j:[0..4)</figcaption></figure>

If on the contrary the compiler choice was parallel execution, the explicit is:

Example 4.

$$
\begin{array}{l}
\{ \text{par}(j=0; j<4; j++) \\
\quad \{ C[(4*i)+j]+=A[(i*n)+k]*B[(4*k)+j]; \} \}
\end{array}
$$

and the execution graph would be the one shown in (figure 5).

<figure><img src="image_6.png" /><figcaption>Figure 5. Execution graph of the par(4) algorithm variant, par j:[0..4)</figcaption></figure>

In complexity theory, those two variants are considered equivalent for work i.e. arith-
metic operation count, but of course parallel time is O(1) for the parallel variant, while it
is O(p) in the sequential variant, assuming p = 4 processing units. But concrete parallel
algorithms either designed by humans [38] at static time, by compilers at compile time or
by the runtime, are never fully parallel but rather block-parallel. There are good reasons for
that: physical data transfers can never occur as fast as local state changes <fnref n="5" />. A block-parallel
version of our example comes by splitting the loop into 2 blocks of 2 sequential lines each.

Example 5.

$$
\begin{array}{l}
\{ \text{par}(\text{proc}=0; j<2; j++) \\
\quad \{ \text{seq}(j\text{seq}=0; j\text{seq}<2; j++) \\
\quad \quad \{ \text{C}[(4*i)+(2*\text{proc})+j\text{seq}]+=\text{A}[(i*n)+k]*\text{B}[(4*k)+(2*\text{proc})+j\text{seq}]; \} \}
\end{array}
$$

The execution graph becomes (figure 5).

<footnote n="5">⁵ This may also be the case for future quantum computers, but this question is beyond the scope of this paper.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Execution graph of the parseq(2,2) algorithm variant, par proc:[0..2) seq jseq:[0..2)</figcaption></figure>

This type of block decomposition can happen in a variety of combinations depending on block sizes, which loop is being split, whether the seq-par alternation is nested etc. In [27] we have exhaustively explored this combinatorial space on example 2 and target architectures with the two lowest levels of parallelism, namely vector registers and multiple cores. Analysis of the timings reveals that optimal choices are made when algorithm block sizes are exactly or close to the architecture's vector register size and/or number of cores. This behavior is known to apply to all levels of parallel architectures: multiple nodes, GPU/TPU processing units etc. As a result a mathematical theory can be built that relates algorithm shape with estimated processing time. That theory has two parts:

1. Producing explicit seq-par loops from MoA expressions. That is called *dimension lifting* because it increases the number of dimensions in one of the arrays [27], and corresponds to the transformation from a *par* to a *par-seq* loop seen above. Dimension lifting is a purely algebraic transformation applied to the DNF, using the reshape $\rho$ operator and its effect on the ONF or fragment-C code is to produce blocks that are executed in parallel. The space of possible dimension liftings is essentially that of possible reshaping operations on the given algorithm. The laws of linear algebra ensure that all such transformations preserve the algorithm's result.

2. Estimating the *cost* i.e. estimated execution time of each algorithm variant. This builds on algorithmic theories like BSP [40] and its hierarchical variant the scatter-gather language SGL [12]: the cost of an execution graph combines instruction counts *and* transfer (effects of bandwidth, or synchronization i.e. latency) times. In practice the cost can be either normalized in units of arithmetic operations (assuming that their rate per second is mostly stable) or purely in (fractions of) seconds. The appropriate theory for defining costs comes from operations research, namely weights from a *semiring* that label the edges/nodes of a directed graph [41].

Let us apply this concept to our execution graphs above, in a manner that is a pure application of the standard concept of *semiring costs* (e.g. max+ algebras) and a realistic and flexible tool for actual execution speeds on computers.

Counting arithmetic operations poses no conceptual problem. The basic step in our algorithm is of the form $C[(i*p)+j]+= A[(i*n)+k]*B[(k*p)+j]$ so it requires 4 multiplications and 4 additions (3 for indices and one for accumulation $+=$). Let us assume that this amounts to 8 units of time.

The *semiring* is $(\mathbf{N}, \max, +, 0, 0)$ so that we associate a unit of time to every edge or node being traversed, accumulate them by addition on any path and compute the maximal cost for a path in the algorithm directed acyclic graph (dag). But applying this theory to pure arithmetic costs would not satisfy the physical requirement for transfer times.

That mathematical requirement is trivial: to associate costs not only to the nodes being traversed, but also to the edges. The scientific problem however is extremely difficult: to
associate such costs in a manner that it models actual measurements on computer architectures. This problem is compounded by the great variety and non-scientific design culture of processors and computers. But the decades of experiments with multinode systems and the BSP cost model [11], our recent experiments with GPU [42] and vector-multicore CPU systems [27] confirm all of the above remarks: linear-algebraic combinations of transfer and arithmetic costs, structured in hierarchical blocks that match the hardware structures define a tractable and realistic search space for optimal and predictable performance.

Let us therefore estimate all costs for our algorithm variant examples, and generalize them to a realistic but very general theory. This is a concrete thought experiment using simplistic constants but being based on standard hardware behavior it brings out one possible general algebraic form for the costs. Assume that the target architecture has two processing units. We have seen that every basic execution block requires 8 arithmetic time units. It also requires one read (load) from array B, one write (store) to array C, and *at most* one read (load) from array A. Assume that the cost of an arithmetic operation is one unit of time, and that of a transfer operation is *t*: N. Realistic values of *t* can range from a dozen to thousands depending on the specific hardware level.

[seq] When executing the seq algorithm variant, only one processing unit is active and it may "prefetch" the A element $A[(i*n)+k]$ i.e. load it only once because it is independent of loop index *j*. Each one of the four sequential steps requires a read from B and write to C. So the overall arithmetic cost is that of read(A) followed by four time (read(B);write(C)) i.e. $9t$. The arithmetic cost is the sum of the cost of one step $4 \times 8$.

[par] When executing the par algorithm variant, we exceed the capacity for parallel execution by creating four parallel processes or threads. As a result, the two processing units will interleave the four steps. The degree of parallelism implies that the arithmetic cost is the sum of the costs of two sequential steps. In other words the arithmetic cost of one thread which is $2 \times 8$. For transfer costs our model is as follows: each thread may load(A) a single time, but there is no guarantee that the hardware can process this concurrent read in one step. Our modeling is that it could cost $2t$ if that operation is serialized. Each thread reads from B twice, and that is fully independent between the threads because the array elements are distinct. Each thread also writes to C twice, also a fully parallel operation. As a result the transfer cost is $2t + 2t + 2t = 6t$.

[parseq] When executing the parseq algorithm, the loops create exactly two threads of two sequential steps each, which matches hardware size. The arithmetic cost is the same as the (two-thread) execution of par: $2 \times 8$. However, each thread is able to amortize the loading of $A[(i*n)+k]$ because the sequence of two operations is known at compile time, and not the effect of dynamic thread scheduling. So each processing unit reads from A only once, then reads twice from B and writes twice to C. The transfer cost is then $5t$. The costs are summarized in the following table.

<table><thead><tr><td>algorithm variant</td><td>arithmetic</td><td>transfer</td><td>cost</td></tr></thead><tbody><tr><td>seq(4)</td><td>4 × 8</td><td>9t</td><td>9t + 32</td></tr><tr><td>par(4)</td><td>2 × 8</td><td>6t</td><td>6t + 16</td></tr><tr><td>parseq(2,2)</td><td>2 × 8</td><td>5t</td><td>5t + 16</td></tr></tbody></table>

Two conclusions can be drawn from our thought experiment. They are specific to this example but, because of the physical properties of hardware they are experimentally verified:

1. Sequential execution is slowest. This holds when $t > 1$, and this is true except in very unusual contexts.
2. Block-conscious parallel algorithm variant is fastest because it maximizes the use of static information about data- and architecture *shapes*. In our tiny example this led to the economy if one load operation while retaining gain.

MoA dimension lifting can thus explore block-parallel variants of a given algorithm with the assurance of finding a cost that is optimal for the standard arithmetic+transfer measurements. There is a huge side-effect of this methodology which is not explicit here: some parallel compilation methods attempt to match a random computation graph to a random architecture structure, which leads to an intractable problem similar to job scheduling. By contrast, MoA-based exploration of algorithm variants by reshaping leads to a tractable and usually optimal cost.

### 4.4. General costs

We can now generalize from the above observations and summarize our applied mathematical basis for MoA-based parallel computing.

A given mathematical expression on arrays, expressed with MoA can be transformed to DNF, which is equivalent to an ONF and in-turn a fragment-C loop-nest program. That program's computation dag is highly symmetric and need not be enumerated: its loop sizes (array shapes) are sufficient to describe it.

Many equivalent variants of the computation dag can be produced by dimension-lifting. This search space can be restricted to block sizes that match the architecture: size of the vector registers, number of cores, number of supercomputer nodes etc. But it can cover the combinatorial space of dimensions: each subset of the array dimensions can be "dimension lifted" i.e. split into blocks. For example a two-dimensional array can have its rows divided, or its columns, or both. The first two options produce a 3-dimensional array and the last one produces a 4-dimensional array.

For each computation dag derived in this manner, a cost estimate can be produced by applying some experimental estimate and architecture knowledge: parallel or sequential execution, load-store protocols and delays, communication or synchronization delays etc. Once the costs for each step and transfer are produced, they can be combined with the semiring operations: max for parallel parts of the graph, and + for successive parts.

Here are arguments for using the general theory of semirings in this context, the max+ algebra being the most typical example, but not general enough.

* The $N$ values are obviously too simplistic. If costs are measured in seconds, then fractions will appear so the base type can be that of rationals.

* If the architecture is hierarchical then so should the computation dags be. In that case, the semiring elements could be vectors of costs, combined with pointwise addition and vector-to-rational maximum.

* For abstract reasoning it may be possible to define the graphs from regular expressions. In which case the semiring elements become languages (sets) of abstract costs, + is replaced by concatenation and max is replaced by union.

* In modern processors it often happens that the clock rate changes over time, so the constant arithmetic rate hypothesis becomes unrealistic. Cost then become time series, + can be replaced by concatenation and max by the maximum of some integral or statistics on the time series.

For all the above reasons we introduce the following definition for parallel execution costs. It is fully axiomatic and we have provided arguments that it leads to a fruitful interaction between algebra and performance evaluation/optimization.

**Definition 24.** A parallel execution cost structure $(\mathcal{G}, (C, \oplus, \odot, \emptyset, \mathcal{Z}), \ominus)$ consists of
1. An infinite family $\mathcal{G}$ of finite directed acyclic graphs $g \in \mathcal{G}$ whose nodes $x$ are called instructions, whose edges $(x, y)$ are called (data) transfers and such that every $g \in \mathcal{G}$ has a single source and single sink.

2. A semiring $(C, \oplus, \odot, \emptyset, Z)$ of costs with a well-founded total order $<$ on $C$ such that $x \le (x \oplus y), x \le x \odot y$ and $x \le x \odot y$.

3. A cost function that associates an element of $C$ to every node and edge of any graph $g \in \mathcal{G}$.

4. An $\{\oplus, \odot\}$-coloring $\alpha$ of the nodes of every graph $g$ (except the sink) such that if $\alpha(x) = \odot$ then $x$ is called a sequential node and has exactly one outgoing edge, if $\alpha(y) = \oplus$ then $y$ is called a parallel node and has more than one outgoing edges.

The example execution graphs above illustrate trivial cases of the family of graphs, the two-coloring and the costs. General execution cost is defined on sequential-parallel graphs of any nesting depth and size, and the following general definition of graph cost.

**Definition 25.** *Given a parallel execution cost structure $(\mathcal{G}, (C, \oplus, \odot, \emptyset, \mathcal{Z}), \odot)$ the cost $C$ of a dag $g \in \mathcal{G}$ is defined recursively by its source node:*

* if $z$ is the graph sink then $C(z) = \text{cost}(x) \odot Z$.

* if $(x, x')$ is the unique outgoing edge from a sequential node $x$ then $C(x) = \text{cost}(x) \odot C(x')$.

* if $(y, y_1), (y, y_2), \dots, (y, y_k)$ are the outgoing edges from a parallel node $y$ then $C(y) = \text{cost}(y) \odot (C(y_1) \oplus \dots C(y_k))$.

The definition of graph cost we have used in our examples uses the max-+ semiring $(\mathbb{N}, \max, +, 0, 0)$, once the execution and transfer costs are evaluated to natural numbers. It is also possible to use a semiring of time-dependent natural numbers (to model variable clock rate in the architecture), a semiring of rational numbers (to model actual delays in seconds and not in instructions counts) or any complex combinations that the axioms allow.

Finally we give an abstract form of the effect of MoA dimension-lifting on computation costs, that may inspire more general applications and theorems.

**Definition 26.** *A dimension lifting on a cost structure is an internal map $L: \mathcal{G} \to \mathcal{G}$, such that (a) preserves edges and nodes ($L(g)$ is a permutation of $g$) (b) preserves the cost of nodes and edges. (c) reduces the length of every sequential chain*

**Theorem 2.** *For a dimension lifting $L$ on $(\mathcal{G}, (C, \oplus, \odot, \emptyset, \mathcal{Z}), \odot)$ and any graph $g \in \mathcal{G}$, then $C(L(g)) \le C(g)$.

**Proof** *Proof: by the recursive definition of (sub)graph cost the complete graph's cost $C(x)$ is an expression on operators $\oplus, \odot$. By definition of $L$, $L(x)$ contains the same semiring values but every sequential chain of nodes is shorter. By the properties $x \le x \odot y$ and $x \le y$ it follows that the cost of every sequential chain starting at a node $x_0$ is less in (the graph rooted at) $L(x_0)$ than in (the graph rooted at) $x_0$. Given a parallel node $y_0$, its local cost $\text{cost}(y_0)$ is preserved by definition of $L$. Moreover, by induction on the depth of the subgraph, the cost of its subgraphs $y_1, \dots, y_k$ are also decreased through $L$. As a result $C(L(y_0)) =$*

$$ \text{cost}(L(y_0)) \odot (C(L(y_1)) \oplus \dots C(L(y_k))) = $$

$$ \text{cost}(y_0) \odot (C(L(y_1)) \oplus \dots C(L(y_k))) \le $$

$$ \text{cost}(y_0) \odot (C(y_1) \oplus \dots C(y_k)). $$

QED.
# 5. Discussion

Another example of the abstract dimension-lifting uses so-called *binary transpose*. Binary transpose permutes chosen dimensions and their indices. It is defined by a permutation of the dimensions, for example [0, 2, 1, 3] if the array has four dimensions. In the example 1 we explained how it can combine with the reduce operator to allow the reduction of any dimension in the array, while pure reduction only reduces the first dimension.

It can also combine with the reshape operator to drive the algorithm mapping to processors. In other words, reshaping and permutations define the MoA expression's ONF so indirectly, the loop bounds of the equivalent fragment-C and finally the shape of an execution graph for this pseudo code.

In section 4 we gave examples of mapping a single loop to processors, and the block structures in those execution graphs can be created by simple transpose-reshape operations. Those examples were 1-dimensional but suppose we would want to map blocks that are two-dimensional. Assume an array $\zeta$ of shape $\langle 4, 4 \rangle$ mapped to two processors, we could reshape it to $\langle 2, 2, 2, 2 \rangle$ i.e. the shape of a square 2x2 array of 2x2 blocks $^6$. By permuting dimensions, i.e $\langle 0, 2, 1, 3 \rangle$ transpose of a $\langle 2, 2, 2, 2 \rangle$-reshaped $\zeta'$, we also permute the indices, thus creating blocks. Thus for any index $0 \le \tilde{i} < (\rho\zeta')$ in the reshaped array each component would be accessed by $\langle 0, 2, 1, 3 \rangle \psi\zeta'$.

Why is this blocking useful? Blocking identifies a size relevant to the speed of the closest processor memory, e.g. cache. The blocking size would be related to the cache (or whatever the closest memory unit is) and of course to the elements' data type and size. As levels of memory emerge in the formulation of costs for arbitrary machines, the blocking would continue recursively as the properties of linear algebra allow and as seen for example in Strassen matrix multiplication.

# 6. Conclusions and future work

We have shown how the MoA array algebra is a general and suitable basis for a mathematical theory of parallel computing. From algebra of algorithms to execution costs, a methodology has been defined and its mathematical basics defined.

Parallel compilation should become an experimental science and we want to define the basis for that evolution. There are enormous implications for that branch of applied mathematics, a blend of algebra, operations research and physics of computation including group theory.

Future research directions include (a) non-arithmetic preserving transform in MoA (Winograd- or Strassen-like optimizations that reduce the arithmetic operation count) to generalize the dimension-lifting transforms; (b) the mechanization of MoA and of MoA-based compilation; (c) the application of Vitanyi's [43] results that can relate the orbits of symmetry groups on MoA expressions to lower-bounds on lengths in space-time.

**Author Contributions:** Conceptualization G.H. and L.M.; methodology G.H and L.M.; formal analysis (algebra) L.M.; formal analysis (semantics) G.H.; high-performance programming and tuning L.M.; data analysis G.H.; original draft preparation, G.H and L.M.; writing—review and editing, G.H and L.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** The experimental work [27] that has led to this mathematical research have been funded by the Stony Brook University Ookami computing center and by CRIANN the *Centre Régional Informatique et d'Applications Numériques de Normandie* under project 2024004. Ookami is supported by the US NSF grant #1927880. CRIANN is supported by *Région Normandie*, the French government and the European union.

<footnote n="6">⁶ This type of blocking is ideal for algorithms like the Kronecker Product.</footnote>
**Data Availability Statement:** There is no dataset specific to this paper's content. Our existing datasets contain performance measurements and not the mathematical foundations presented here.

**Acknowledgments:** The authors thank their research partners for much inspiration and related work: Ian Grout, Allan Porterfield, John Mullins and others.

**Conflicts of Interest:** The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

**Abbreviations**

The following abbreviations are used in this article:

<table><tr><td>N</td><td>The totally-ordered set of non-negative integers</td></tr><tr><td>Z</td><td>The totally-ordered set of all integers</td></tr><tr><td>S*</td><td>The set of possibly-empty lists of elements of S</td></tr><tr><td>_++_</td><td>The concatenation operation on lists.</td></tr><tr><td>A<sup>□</sup></td><td>The set of finite, possibly-empty arrays of elements of A</td></tr><tr><td>HPC</td><td>High-Performance Computing</td></tr><tr><td>MoA</td><td>A Mathematics of Arrays</td></tr><tr><td>BSP</td><td>Bulk-Synchronous Parallel (computing)</td></tr><tr><td>GEMM</td><td>General Matrix Multiply</td></tr><tr><td>MM</td><td>Matrix Multiply</td></tr><tr><td>HP</td><td>Hadamard Product</td></tr><tr><td>KP</td><td>Kronecker Product</td></tr><tr><td>DNF</td><td>Denotational Normal Form</td></tr><tr><td>ONF</td><td>Operational Normal Form</td></tr><tr><td>ζ.ρ</td><td>Shape list of array ζ</td></tr><tr><td>ζ.κ</td><td>Content list of array ζ</td></tr><tr><td>ρ_</td><td>The shape operator.</td></tr><tr><td>δ_</td><td>The dimensionality operator.</td></tr><tr><td>τ_</td><td>The total size operator.</td></tr><tr><td>_↑_</td><td>Overloaded: list-take or array-take operator.</td></tr><tr><td>_↓_</td><td>Overloaded: list-drop or array-drop operator.</td></tr><tr><td>ι_</td><td>Overloaded: the index generator on non-negative integers or on arrays.</td></tr><tr><td>_[_]</td><td>Overloaded: flat indexing or one-integer indexing.</td></tr><tr><td>_Ψ_</td><td>The (generalized) indexing operator.</td></tr><tr><td>rav_</td><td>Ravel, the flattening operator.</td></tr><tr><td>γ(_,_)</td><td>The index conversion operator.</td></tr><tr><td>_ρ̂_</td><td>The reshape operator.</td></tr><tr><td>red⊕_</td><td>The reduce operator for binary operation ⊕.</td></tr></table>

**References**

1. Jenkins, M.A.; Mullin, L.R., A Comparison of Array Theory and a Mathematics of Arrays. In *Arrays, Functional Languages, and Parallel Systems*; Mullin, L.M.R.; Jenkins, M.; Hains, G.; Bernecky, R.; Gao, G., Eds.; Springer US: Boston, MA, 1991; pp. 237–267. https://doi.org/10.1007/978-1-4615-4002-1_14.

2. Markus, A.; Mullin, L. Implementing Mathematics of Arrays in Modern Fortran: Efficiency and Efficacy. Technical report, TechXriv, 2023. DOI: 10.36227/techrxiv.24058647.v1.

3. Kennedy, K.; Koelbel, C.; Zima, H. The rise and fall of High Performance Fortran: an historical object lesson. In Proceedings of the Proceedings of the third ACM SIGPLAN conference on History of programming languages, 2007, pp. 7–1.

4. Ritchie, D.M.; Johnson, S.C.; Lesk, M.; Kernighan, B.; et al. The C programming language. *Bell Sys. Tech. J* **1978**, *57*, 1991–2019.

5. Clarke, L.; Glendinning, I.; Hempel, R. The MPI message passing interface standard. In Proceedings of the Programming Environments for Massively Parallel Distributed Systems: Working Conference of the IFIP WG 10.3, April 25–29, 1994. Springer, 1994, pp. 213–218.
6. Dagum, L.; Menon, R. OpenMP: an industry standard API for shared-memory programming. *IEEE computational science and engineering* **1998**, 5, 46–55.

7. Farber, R. *CUDA application design and development*; Elsevier, 2011.

8. Farber, R. *Parallel programming with OpenACC*; Newnes, 2016.

9. foundation, P. Pytorch documentation. https://docs.pytorch.org/docs/stable/index.html.

10. Valiant, L.G. Why BSP computers?(bulk-synchronous parallel computers). In Proceedings of the [1993] Proceedings Seventh International Parallel Processing Symposium. IEEE, 1993, pp. 2–5.

11. McColl, W.F.; Tiskin, A. Memory-efficient matrix multiplication in the BSP model. *Algorithmica* **1999**, 24, 287–297.

12. Li, C.; Hains, G. SGL: Towards a bridging model for heterogeneous hierarchical platforms. *International Journal of High Performance Computing and Networking* **2012**, 7, 139–151. https://doi.org/10.1504/IJHPCN.2012.046371.

13. Gerbessiotis, A.V. Extending the BSP model for multi-core and out-of-core computing: MBSP. *Parallel Computing* **2015**, 41, 90–102.

14. Allombert, V.; Gava, F.; Tesson, J. Multi-ML: programming multi-BSP algorithms in ML. *International Journal of Parallel Programming* **2017**, 45, 340–361.

15. Legaux, J.; Hu, Z.; Loulergue, F.; Matsuzaki, K.; Tesson, J. Programming with BSP homomorphisms. In Proceedings of the European Conference on Parallel Processing. Springer, 2013, pp. 446–457.

16. Williams, S.; Waterman, A.; Patterson, D. Roofline: an insightful visual performance model for multicore architectures. *Communications of the ACM* **2009**, 52, 65–76.

17. Ofenbeck, G.; Steinmann, R.; Caparros, V.; Spampinato, D.G.; Püschel, M. Applying the roofline model. In Proceedings of the 2014 IEEE International Symposium on Performance Analysis of Systems and Software (ISPASS), 2014, pp. 76–85. https://doi.org/10.1109/ISPASS.2014.6844463.

18. Treibig, J.; Hager, G. Introducing a performance model for bandwidth-limited loop kernels. In Proceedings of the International Conference on Parallel Processing and Applied Mathematics. Springer, 2009, pp. 615–624.

19. Rajopadhye, S.V., Dependence analysis and parallelizing transformations. In *The compiler design handbook*; Srikant, Y.N.; Shankar, P., Eds.; CRC Press: Publication Location, 2003.

20. Mullin, L.M. A Mathematics of Arrays. PhD thesis, Syracuse University, Syracuse, NY, 1988.

21. Berkling, K. Arrays and the Lambda Calculus. Technical report, Syracuse University, 1990.

22. Hains, G.; Mullin, L.M.R. Parallel functional programming with arrays. *The Computer Journal* **1993**, 36, 238–245.

23. Hunt, H.; Mullin, L.; Rosenkrantz, D.; Raynolds, J. A Transformation-Based Approach for the Design of Parallel/Distributed Scientific Software: the FFT. *CoRR* **2008**, abs/0811.2535, [0811.2535].

24. Mullin, L.R.; Raynolds, J.E. Conformal Computing: Algebraically connecting the hardware/software boundary using a uniform approach to high-performance computation for software and hardware applications, 2008, [arXiv:cs.MS/0803.2386].

25. Grout, I.; Mullin, L. Realization of the Kronecker Product in VHDL using Multi-Dimensional Arrays. In Proceedings of the 2019 7th International Electrical Engineering Congress (iEECON), 2019, pp. 1–4.

26. Grout, I.; Mullin, L. Realizing Mathematics of Arrays Operations as Custom Architecture Hardware-Software Co-Design Solutions. *Information* **2022**, 13.

27. Mullin, L.; Hains, G. Towards Automatic, Predictable and High-Performance Parallel Code Generation. In Proceedings of the Parallel and Distributed Processing Techniques (CSCE/PDPTA2024); Arabnia, H.R.; Takata, M.; Deligiannidis, L.; Rivas, P.; Ohue, M.; Yasuo, N., Eds. Springer Nature Switzerland, 2025, pp. 109–119.

28. Hains, G.; Fenek, O. Machine Learning Pseudo-Natural Language for Temporal Logic Requirements of Embedded Systems. In Proceedings of the 2023 15th International Conference on Knowledge and Systems Engineering (KSE). IEEE, 2023, pp. 1–4.

29. Grout, I.; Mullin, L. Hardware Considerations for Tensor Implementation and Analysis Using the Field Programmable Gate Array. *Electronics* **2018**, 7. https://doi.org/10.3390/electronics7110320.
30. Thomas, S.; Mullin, L.; Swirydowicz, K.; Khan, R. Threaded Multi-Core GEMM with MoA and Cache-Blocking. In Proceedings of the 2021 World Congress in Computer Science, Computer Engineering, & Applied Computing (CSCE'21), 7 2021. 1148 1149 1150

31. Mullin, L.R. A uniform way of reasoning about array-based computation in radar: Algebraically connecting the hardware/software boundary. Digit. Signal Process. 2005, 15, 466-520. https://doi.org/10.1016/j.dsp.2005.02.003. 1151 1152 1153

32. Berkling, K. Arrays and the Lambda Calculus. Technical Report 93, Syracuse University, Syracuse, NY, USA, 1990. SU-CIS-90-22. 1154 1155

33. Ostrouchov, C.; Mullin, L. PythonMoA. https://labs.quansight.org/blog/2019/04/python-moa-tensor-compiler/. 1156 1157

34. Thomas, S.; Mullin, L.; Świrydowicz, K. Improving the Performance of DGEMM with MoA and Cache Blocking. In Proceedings of the Arrays2021. ACM, 2021. 1158 1159

35. Grout, I.; Mullin, L. Realization of the Kronecker Product in VHDL using Multi-Dimensional Arrays. In Proceedings of the 2019 7th International Electrical Engineering Congress (iEECON). IEEE, mar 2019. https://doi.org/10.1109/ieecon45304.2019.8938846. 1160 1161 1162

36. Hains, G.; Mullins, J.; Mullin, L. Language Foundations for HPC Array Structures. In Proceedings of the CSCE/FCS'25 - 21st Int. Conf. Foundations of Computer Science. Springer Nature Switzerland, 2025. To appear. 1163 1164 1165

37. Mullin, L.M. A Mathematics of Arrays 1988. 1166

38. Gaétan, H. Algorithmes et programmation parallèles-Théorie avec BSP et pratique avec OCaml; Editions Ellipses: Paris, 2018. 1167 1168

39. Hager, G.; Treibig, J.; Habich, J.; Wellein, G. Exploring performance and power properties of modern multi-core chips via simple machine models. Concurrency and computation: practice and experience 2016, 28, 189-210. 1169 1170 1171

40. Valiant, L.G. Why BSP computers?(bulk-synchronous parallel computers). In Proceedings of the [1993] Proceedings Seventh International Parallel Processing Symposium, 1993, pp. 2-5. 1172 1173

41. Gondran, M.; Minoux, M. Graphs, dioids and semirings: new models and algorithms; Springer, 2008. 1174

42. Mullin, L.; Phan, W. A Transformational Approach to Scientific Software: The Mathematics of Arrays (MoA) FFT with OpenACC, 2021. 1175 1176

43. Vitányi, P.M. Locality, communication, and interconnect length in multicomputers. SIAM Journal on Computing 1988, 17, 659-672. 1177 1178

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content. 1179 1180 1181 1182