Easy Multi-Precision Acceleration via Vectorization

Easy Multi-Precision Acceleration via Vectorization Chansup Byun, Piotr Luszczek, Jeremy Kepner, LaToya Anderson, William Arcand, David Bestor, William Bergeron, Alex Bonn, Daniel Burrill, Vijay Gadepally, Michael Houle, Matthew Hubbell, Hayden Jananthan, Michael Jones, Lauren Milechin, Julie Mullen, Andrew Prout, Albert Reuther, Antonio Rosa, Charles Yee, Peter Michaleas

I. I NTRODUCTION Easy-to-use high level programming languages (e.g., Matlab [1], Python [2], and Julia [3],) enable millions of programmers and play an important role in computational science and high performance computing. Accelerators, most notably GPUs (Graphics Processing Units), have become key enablers for a wide range of computationally intensive workloads, particularly in AI (Artificial Intelligence) [4], [5]. The capability of CPUs (Central Processing Units) continues to increase Research was sponsored by the Department of the Air Force Artificial Intelligence Accelerator and was accomplished under Cooperative Agreement Number FA8750-19-2-1000. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Department of the Air Force or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein.

1.00E+17 1017

Blackwell

Floating Point Operation per Second

+

E

E M

1.00E+16 1016

Hopper +

E

E

E

180x

E M M M

1015 1.00E+15

Volta

1014 1.00E+14

+

E

E

E

E

E M M M M M M M M M M

2015

2020

2025

Year

2030

1.00E+14 1014

Bytes per Second

Blackwell

1.00E+13 1013

Hopper

bandwidth

1.00E+12 1012

Volta

Abstract—High level programming languages and accelerators are powerful enablers for a wide range of applications. Accelerators offer significant performance via diverse numeric precision and increased memory bandwidth and capacity. These enable new algorithms to be explored. Construction of these algorithms can be enabled by the ready accessibility of these performance features in high-level languages. Achieving scalable vertical (within a compute node), precision (over various numeric precisions), and temporal (over different generations of hardware) performance while retaining productivity requires effective abstractions. Vectorized coding that expresses operations as highlevel matrix mathematics has historically been an important method for achieving portable performance. This paper examines the performance availability across a wide range of precisions and hardware using simple vectorized coding techniques applied to several example programs that highlight bandwidth, real and complex matrix multiplication, and 2D convolution. The systems evaluated span nearly a decade of computing hardware, enabling direct comparison of performance gains over this period. This work demonstrates that near hardware peak performance is readily accessible from high-level programming environments. In particular for memory bandwidth, real and complex matrix multiplication in double/single/half precision, and 2D convolutions using the Fast Fourier Transform (FFT). Furthermore, algorithms that trade decreased compute intensity for increased memory intensity have significant potential. Index Terms—GPUs, distributed arrays, parallel processing, vertical scaling, horizontal scaling

Bytes

arXiv:2611.05558v1 [cs.PF] 3 Nov 2026

MIT

1.00E+11 1011

8.5x

70x

8.5x

capacity

1010 1.00E+10 2015

2020

Year

2025

2030

Fig. 1. Precision, Performance, Bandwidth, and Capacity Trends. Increasing performance of Nvidia GPUs as they offer optimized hardware for lower precisions (top) and larger memory (bottom). (top) + denotes a sign bit, E denotes an exponent bit, and M denotes a mantissa bit.

as the number of processing cores on a CPU grows. Accelerators offer significant performance via diverse numeric precision and increase memory bandwidth. A key challenge is maintaining the ease-of-use of high level programming languages while achieving scalable performance across different dimensions: vertical (within a compute node), precision (FP64, FP32, FP16, ...), and temporal (over different generations of hardware). Vertical scaling within nodes often relies on threading approaches (e.g., OpenMP [6], pthreads [7], and Cilk [8]) that leverage shared memory model. Threading allows libraries (e.g., MKL [9], MAGMA [10], SPIRAL [11], and GraphBLAS [12], [13]) written in low-level languages (e.g., C [14], C++ [15], and Fortran [16]) to deliver parallel performance at

the function level that can be readily called by a high level programming language. Vertical scaling on accelerators often follows a similar approach leverage technologies that are designed for more heterogenous accelerator architectures (e.g., CUDA [17] and OpenCL [18]). Accelerators have created a veritable zoo of precisions as they pursue the relatively low precision requirements of AI (see Figure 1). Added to the traditional double (FP64), single (FP64), and integer (int64, int32, int16, int8) types are half (FP16), quarter (FP8), eighth (FP4), brain float (BF16), and tensor float (TF32), which takes BF16 input and produces FP32 output). The potential performance gains of these new types are significant and it is important for algorithm developers to have performant access to them in high level languages. Accelerators have also increased both memory bandwidth and memory capacity of their fastest programmable memories (see Figure 1). The combined memory bandwidth and capacity improvements increasingly enables algorithms that trade decreased compute for increased memory usage. Specifically, algorithms that do fewer operations per byte but access those bytes more (e.g., sparse and spectral methods). Vectorization [19], [20] is an important optimization technique that is particularly relevant to high-level languages that often use large arrays as a core construct. Operating on large arrays can often be the key to leveraging the underlying high performance libraries that are built with the aforementioned vertical scaling methodologies. Using vectorization to access accelerator performance begins with good memory bandwidth performance as high level languages can often use more memory than lower level languages. The Stream benchmark is a common benchmark for measuring bandwidth that vectorizes well [21]. This work builds on our prior work with this benchmark [22]. Diverse precision performance first requires the availability of the corresponding types in the high level language. Next, the extent that high level languages makes it easy to access the underlying hardware performance can be assessed. Matrix multiply is a good tool for exploring both as matrix multiply is a defined mathematical operation for diverse types (both real and complex) and the underlying hardware is often design to accelerate this operation. Together, Stream vector operations and matrix multiplication are a good foundation for building many real applications, particularly in the area of sensor processing. Two important sensor processing applications are beamforming and image convolution. Beamforming is the basis of many sensing and communication applications, including: radar, mineral exploration, medical imaging (MRI, ultrasound, . . . ), wireless communications, land/sea/undersea/air/space manned/unmanned navigation and traffic management, satellite imaging, weather monitoring, oceanography, and radio astronomy. Correspondingly, we have a developed and published a “Beamformer” mini-application that exercises the most computationally intensive aspects of beamforming [23]. Beamforming relies on complex matrix multiplication over diverse precisions and should track the corresponding underlying matrix multiplication performance.

Image (2D) convolution is a foundation of most image processing that is widely use in many of the same applications as beamforming. 2D convolution is also critical to convolutional neural networks (CNNs) widely used in AI applications. Likewise, we have a developed and published a “Blurimage” mini-application that exercises the most computationally intensive aspects of 2D image convolution [23]. Blurimage relies on 2D convolution over diverse precisions (both real and complex) and can be implemented with compute intensive direct methods and memory intensive spectral methods using the Fast Fourier Transform (FFT) [24]–[26]. This paper presents a brief case study on the performance of accelerated vectorized codes over diverse precisions in high level languages, using Stream, matrix multiplication, Beamformer, and Blurimage. These benchmarks are implemented in Matlab and Python and leverage the pMatlab [23] and pPython [27] distributed array libraries for multi-process parallelism when required. Performance is measured on a invidual nodes of varying CPU and GPU hardware spanning almost a decade of technology. The remainder of the paper is organized as follows. We first describe basic vectorization and parallelization. We then present an in-depth look at the matrix multiplication benchmark. Subsequent sections describe the application benchmarks, the benchmarking hardware, and the performance results, followed by conclusions and a discussion of potential future work. II. V ECTORIZATION & PARALLELIZATION The essence of vectorization in high level languages is removing “for loops” and replacing them with higher level operations on whole vectors or matrices that can take advantage of underlying hardware optimized libraries. In the case of Stream a good example is replacing for i = 1 : N do a(i) = b(i) + q*c(i) with a = b + q*c where a, b, c are N element vectors and q is a scalar. In the case of the Stream benchmark, there is an also outer loop that is run Nt times. Similarly, matrix multiplication of N ×N matrices A, B, C can be vectorized by replacing for i = 1 : N , j = 1 : N , k = 1 : N do C(i,j) += A(i,k) * B(k,j) with C=A*B Code vectorization can dramatically increase performance, reduces code size, simplifies debugging, and is often highly portable across generations of hardware. Furthermore, these operations work well across a wide range of types and

block columns

block rows

map grid: 1x4

map grid: 4x1

Block Columns & Rows

block rows with overlap

map grid: 2x2

map grid: 1x4 overlap: NK

N0

N0

N0

A:

N0xN0

N0

0

1

…

NP-1

B:

N0xP(N)

N0

0

1

…

NP-1

C:

N0xP(N)

N0

0

1

…

N/NP

N/NP

C = A*B

NP-1

N/NP

Fig. 3. Parallel Matrix Multiplication Design. A is a normal N0 ×N0 matrix replicated across all processes PID ∈ {0, ..., NP − 1}. B and C are N0 ×N distributed matrices with N0 ×N/NP local matrix on each process. If B and C have the same parallel map, then the resulting program will require no communication.

III. M ATRIX M ULTIPLICATION Fig. 2. Distributed Array Mappings (adapted from [23]). Different parallel mappings of a two-dimensional array. Arrays can be broken up in any dimension. A block mapping means that each PID holds a contiguous piece of the array. Overlap allows the boundaries of an array to be stored on two neighboring PID s.

can provide a simple way to access underlying accelerated hardware of performance. For example, simply changing the type of the vectors or matrices is typically all that is required to leverage different underlying computing hardware. When comparing CPU and GPU programming it is important to use optimal approaches for each. A common difference is that CPUs often achieve maximum performance when many instances of the same program are run on a multi-core CPU, while GPUs often to achieve maximum performance with running one instance of a program per GPU. To easily explore these differences requires a method of to easily run multiple instances of the programs. The distributed array model in many ways is the most natural parallel programming model for vectorized high level languages because it is implemented using arrays, which are often the core data type in a high level language. Briefly, the distributed array model creates distributed arrays in which each process (PID ∈ {0, ..., NP − 1}) stores or owns a piece of the whole array. Additional information is stored in the array so that every PID knows which parts of the array the other PID s have. How the arrays are broken up among the PID s is specified by a map [23], [27], [28]. For example, Figure 2 shows a matrix broken up by rows, columns, rows and columns, and columns with some overlap. Computations on distributed matrices are usually performed using the “owner computes” rule, which means that each PID is responsible for doing a computation on the data it is storing locally. Maps can become quite complex and express virtually arbitrary distributions. One example of this complexity is the case in which a boundary of an array is required by more than one PID and will be implicitly communicated to complete the computation. Another example are pipelines which can be implemented by mapping different arrays to different sets of PID s.

The matrix multiplication benchmark constructs several large matrices and multiplies them. Because the amount of computation performed per matrix element is large, performance is primarily constrained by the underlying hardware arithmetic performance. The matrix multiplication benchmark is defined over three real (or complex) matrices chosen so that their total storage occupies a significant fraction of the available processor memory. The basic algorithm for matrix multiplication is shown in Algorithm 1 (see caption for details). Algorithm 1 (Serial Matrix Multiplicaiton). The algorithm creates three N 2 element matrices and performs multiplication on those matrices. M AT M UL(A, B, C : RN ×N ) C=A*B Although matrix multiplication is conceptually simple, it serves as a useful barometer for estimating the performance of many programs on a given hardware platform. The goal of this benchmark is to measure performance in floating point operations per second. It constructs matrices using various data types to exploit the low precision arithmetic units available in recent GPU hardware, thereby demonstrating the significant performance improvements enabled by modern GPU designs. The parallel matrix multiplication program can be implemented quite easily using distributed arrays (see Figure 3). As long as all the arrays use the same distributed mapping, the resulting program requires no interprocessor communication, since each of the Np processes is designed to perform an independent matrix multiplication on its own local matrices. The no communication requirement can be strictly enforced by making the program use the .loc construct which specifies that only the part of the array that is local to the PID is used. This approach is generally recommended as it guarantees that no communication will take place. If a mapping that required communication were accidentally used, then the resulting program will either produce an error or will fail to validate. This mapping results in the parallel algorithm shown using P () notation in Algorithm 2, which specifies which dimensions of the array are to be divided among different PID s. A

visual representation of the design for the parallel matrix multiplication algorithm is shown in Figure 3. Algorithm 2 (Parallel Matrix Multiplication). The parallel algorithm creates two N0 ×P (N ) distributed arrays and one N0 ×N0 local arrays and performs multiplication operations on just the local part of those arrays. P () denotes the dimension being divided amongst different PID . M AT M UL(A : RN0 ×N0 , B, C : RN0 ×P (N ) ) C.loc = A * B.loc Typically N = N0 NP so the amount of work per process remains constant as the number of processes NP is increased.

A. Code

Code excerpts for distributed array matrix multiplication implementations in Matlab and Python are shown in Code Listings 1 and 2. The parallel code can be implemented by simply converting the matrices from regular arrays to distributed arrays. This same code is run by every instance with the only difference between these instances is that each has a unique PID that is used to determine which parts of a distributed array belong to each instance (i.e., B.loc and C.loc). The descriptions of Code Listings 1 and 2 are as follows: Lines 1–3 of the Python code import the required Python packages. The CuPy package highlighted in blue enables running on GPUs. Line 4 of the Python code defines the data type of the distributed arrays and other variables used in the code. This is a convenient way to varying the precision of the matrices used in the matrix multiplication so that we can experiment the performance with different precision, especially low precision on GPUs. As written, this Matlab code requires additional functions (single(), half(), complex()) be applied to each matrix to convert to these types. Line 9 creates an alternating vector of 1 and -1 for initializing A so that resulting matrix multiply stays within the dynamic range of lower precisions. Line 11 creates the parallel map BCmap specifying that the N0 ×N element distributed arrays are to have their columns divided equally among all Np over processes with PID = 0, ..., Np-1.

Code Listing 1 (Matlab Parallel Matrix Multiplication). 1 2 3 4 5 lgN0 = 15; lgN = 15; 6 N0 = 2.∧ lgN0; N = Np .* 2.∧ lgN; 7 A0 = 1.0/N0; B0 =1.0; C0 = 1.0; 8 AsignVec = ones(1,N0); 9 AsignVec(2:2:end) = -1.0; 10 11 BCmap = map([1 Np],{},0:Np-1); % Map. 12 % Allocate and initialize matrices. 13 A = (zeros(N0, N0)+A0) .* AsignVec; 14 Bloc = local(zeros(N0, N, BCmap))+B0; 15 Cloc = local(zeros(N0, N, BCmap))+C0; 16 A = gpuArray(A); 17 Bloc = gpuArray(Bloc); 18 Cloc = gpuArray(Cloc); 19 tic; 20 Cloc = A * Bloc; % Multiplication. 21 t = toc; Code Listing 2 (Python Parallel Matrix Multiplication). 1 import numpy as np 2 import cupy as cp 3 import pPython as GPC 4 DT=np.float64 5 lgN0 = 15; lgN = 15 6 N0 = 2.**lgN0; N = Np .* 2.**lgN 7 A0 = 1.0/N0; B0 =1.0; C0 = 1.0 8 AsignVec = cp.ones((N0)) 9 AsignVec[1::2] = -1.0 10 AsignVec = AsignVec.astype(DT) 11 BCmap = Dmap([1,Np],{},range(Np)) # Map. 12 # Allocate and initialize distributed matrices. 13 A =(cp.zeros((N0,N0),dtype=DT)+A0)*AsignVec 14 Bloc=local(zeros(N0,N,map=BCmap,gpu=1,dtype=DT) 15 Cloc=local(zeros(N0,N,map=BCmap,gpu=1,dtype=DT) 16 17 Bloc = Bloc + B0 18 Cloc = Cloc + C0 19 tic = timer() 20 Cloc = A @ Bloc # Multiplication. 21 t = timer()-tic Lines 12–18 create a normal array A and distributed arrays B and C; extract their local parts; and then initialize them to A0, B0, and C0. Because the program requires no communication, the distributed arrays B and C are never actually allocated. Instead, only their local parts Bloc and Cloc are created. This simplifies the program and reduces the total memory required to run the program by eliminating unnecessary duplication. The Matlab Parallel Computing Toolbox (PCT) function gpuArray() and the Python CuPy functions cp.ones(), cp.zeros, and cp.array() (used by pPython zeros()),

X0

sensor response

X1

beamform

X2

accumulate and sum

Standard operation counts for N element vector addition, multiplication of N × matrices, and convolution of a N × matrix with a Nkern ×Nkern . Operation counts are also provided for alternative memory intensive algorithms use sparse or spectral (FFT) methods. The reduction in operation count is the ratio of the standard to the alternative.

Nt

0

1

Nf

TABLE I O PERATION C OUNTS

X3

2

3

Nb

Ns

Nb

Fig. 4. Parallel Beamformer Design. Each stage of the algorithm interacts primarily with two arrays: one for input and one for output. The distributed arrays X0,1,2 are all decomposed in the Nf dimension so that no communication is required in the sensor response, beamform, and write data stages. The final stage, accumulation and summation, requires that all results be sent back to PID = 0. Nt , Ns , and Nb set the sizes of the various arrays used in the calculation.

copy A, Bloc, and Cloc to the GPU memory; subsequent operations on these variables will be automatically performed on the GPU. Lines 19–21 execute the matrix multiplication operations on Aloc and Bloc, and store the results in Cloc. The time for the matrix multiplication operation is summed for each trial. For GPUs a PCT wait command or CuPy synchronize command may need to be called prior to each tic or toc command to ensure the operations are completed on the GPU. In terms of the differences between Matlab and Python, our preliminary observation is that Matlab makes more choices on behalf of the user to significantly lower the barrier to entry. Matlab GPU programming is significantly easier than Python GPU programming. The Matlab performance profiler dramatically simplifies both CPU and GPU performance tuning. There is no equivalent to the Matlab performance profiler in Python. That said, Python tends to expose more performance tuning knobs which can be used to increase performance. In general, we used Matlab to find a good implementation and then converted the design to Python. IV. A PPLICATION B ENCHMARKS We also examined two application benchmarks, Beamformer and Blurimage. These benchmarks are representative of real-world sensor applications used across multiple domains. The Beamformer benchmark involves complex matrix multiplications and serves as an essential kernel in many applications. A visual representation of the parallel Beamformer design is shown in Figure 4. The Blurimage benchmark is a simple image-filtering program that repeatedly convolves an image with a filter (or kernel) using 2D convolution. This is one of the most common operations in image processing. An example of the input and output images produced by the Blurimage application is shown in Figure 5. The input Nx ×Ny matrix consists of a set of random points that are then convolved with a Nkern ×Nkern kernel of all 1’s. A checkerboard patten of 1 and -1 values can be added to the input so that the resulting matrix stays within the dynamic range of lower precisions. Convolution can be implemented using both direct and spectral

Operation Vector Addition Matrix Multiplication 2D Convolution

Standard N

Alternative k

Standard/Alternative N/k

(sparse)

N3

k2 N

(N/k)2

(sparse)

N2 Nkern2

N2 log(N2) (FFT)

Nkern2 / log(N2)

(FFT) algorithms. Our preliminary experiments indicate that FFT convolution is significantly faster than direct convolution even using optimized direct convolvers from AI libraries. Thus, all subsequent results reported here are for FFT based convolution.

Fig. 5. Blurimage Input and Output. The input Nx ×Ny matrix consists of a set of random points that then convolve with a Nkern ×Nkern kernel of all 1’s.

Operation counts for vector addition, matrix multiplication, and 2D convolution operations are presented in Table I for standard and alternative memory intensive algorithms using sparse or spectral (FFT) methods. Where applicable, the constant factor for multiply accumulate on real values is 2 and for complex values it is 8. As shown in Table I the alternative algorithms provide significant benefits for matrix and 2D convolution operations. It is worth mentioning that it has become popular to report the performance of accelerator hardware with SpFlops (Sparse Flops), which uses dense operation counts when performing sparse operations. In some sense this gives credit for work that was not done. However, it has also become a popular way to highlight algorithmic improvements in a way that is more broadly appreciated. The same approach can also be applied to 2D convolutions. Fortunately, it is relatively easy to approximately convert between the two approaches by multiplying/dividing by the appropriate ratio (see column 4 of Table I). V. C OMPUTER H ARDWARE Our team has developed a high-productivity scalable platform—the MIT SuperCloud—for providing scientists and engineers the tools they need to analyze large-scale dynamic

TABLE II C OMPUTER H ARDWARE S PECIFICATIONS MIT SuperCloud maintains a diverse set of hardware running an identical modern software stack providing a unique platform for comparing performance over different eras. GPUs are listed below their host systems. NVIDIA B200 and RTX PRO 6000 GPU systems are hosted at the MGHPCC (Massachusetts Green High Performance Computing Center) AI Computing Resource (AICR). Node Label

Processor

Memory

Era

Part

Clock

Cores

Part

Size

amd-e9-64c

2026

Dual AMD EPYC 9575

3.3 GHz

128

DDR6

2304 GB

- b200

2026

Dual Nvidia B200

2.6 GHz

HBM3e

540 GB

- rtx6000

2026

Dual Nvidia RTX600

2.6 GHz

GDDR6

192 GB

amd-e9-24c

2024

Dual AMD EPYC 9254

2.9 GHz

DDR5

750 GB

- h100nvl

2024

Dual Nvidia H100 NVL

1.7 GHz

HBM3

188 GB

xeon-p8

2020

Dual Xeon Platinum 8260

2.4 GHz

48

DDR4

192 GB

xeon-g6

2018

Dual Xeon Gold 6248

2.5 GHz

40

DDR4

384 GB

- v100

2018

Dual Nvidia V100

1.2 GHz

HBM2

64 GB

48

TABLE III S INGLE N ODE B ENCHMARK PARAMETERS Benchmark parameters were chosen to balance consistency with the memory capacity of the different hardware. Parameters with /NP show the dimension that was scaled with increasing number of processes to keep the relative problem size constant for each process. Node Label

Stream (Nt, N/NP)

Matrix Multiplication (N0, N/NP)

Beamformer (Nt, Ns, Nb, Nf/NP)

Blurimage (Nx/NP, Ny, Nkern)

amd-e9-64c

40, 229

215, 215

9600, 23040, 10240, 4

214, 214, 28

- b200

1000, 230

215, 215

9600, 23040, 10240, 4

214, 214, 28

- rtx6000

1000, 230

215, 215

9600, 23040, 10240, 4

214, 214, 28

amd-e9-24c

40, 229

215, 215

9600, 23040, 10240, 4

214, 214, 28

- h100nvl

1000, 230

215, 215

9600, 23040, 10240, 4

214, 214, 28

xeon-p8

80, 227

215, 215

9600, 23040, 5760, 4

214, 214, 28

xeon-g6

40, 228

215, 215

9600, 23040, 5760, 4

214, 214, 28

- v100

1000, 229

214, 215

9600, 23040, 5760, 4

214, 214, 28

data [29]–[31]. The MIT SuperCloud provides interactive analysis capabilities accessible from high level programming environments (Python, Julia, Matlab) that scale to thousands of processing nodes. MIT SuperCloud maintains a diverse set of hardware running an identical software stack that allows direct comparison of hardware from different eras. A typical benchmarking run can be launched in a few seconds using the MIT SuperCloud triples-mode hierarchical launching system [31] enabling rapid interactive benchmarking. The launch parameters were [Nnode Nppn Ntpp ], corresponding to Nnode nodes, Nppn Matlab or Python processes per node, and Ntpp OpenMP threads per process. The total number of processes is given by NP = Nnode Nppn . On each node, each of the Nppn processes and their corresponding Ntpp threads were pinned to adjacent cores to minimize interprocess contention and maximize cache locality for the OpenMP threads [32]. Within each Matlab or Python process, the OpenMP parallelism is used as provided by their math libraries. At the end of the processing the results were aggregated using asynchronous file-based messaging [33]. Triples mode makes it easy to explore vertical scaling by examining

combinations of processes and threads on a node, and temporal scaling by running on diverse hardware from different eras. The computing hardware includes many different types of nodes that have been deployed over nearly a decade (see Table II). All nodes are multicore x86 compatible. The MIT SuperCloud maintains the same modern software across all nodes, which allows for direct comparison of hardware performance differences. NVIDIA B200 and RTX PRO 6000 GPU systems are hosted at the MGHPCC (Massachusetts Green High Performance Computing Center) AI Computing Resource (AICR). Matlab and Python benchmark distributed array implementations for the different hardware configurations (see Table II) were run with the parameters listed in Table III. The size parameters N , N0 , Ns , Nb , Nf , Nx , and Ny were chosen with a base value which is the maximum PCT and CuPy currently allow for GPU calculations. Within a node N was then scaled with NP (i.e., constant N/NP ) to the maximum memory allowed on the node for Stream and matrix multiplication. A similar approach was used for scaling Nf /NP and Nx /NP for Beamformer and Blurimage. Once the node memory maximum was hit, N was kept constant (i.e., decreasing N/NP ) and Nt was increased to keep the relative run time around a few hundred seconds. VI. P ERFORMANCE R ESULTS Table IV presents the best aggregated within-node benchmark results for all hardware configurations listed in Table II. Stream results are reported in GBytes/s. All other benchmarks are reported in GFlop/s. Blurimage (FFT) results use direct convolution operations counts (see Table I) to highlight the algorithmic performance gains; underlying FFT operation counts 2 can be estimated by dividing by Nkern /(log(Nx ) log(Ny )). The following labeling convention is used: bold much better than FP64 MatMul, blue best precision, * TF32 enabled, ** beyond typical user, ≪FP64 much slower than FP64 MatMul, u unsupported. In general, accessible performance near the underlying hardware limits is available for real/complex FP64/FP32 across all the benchmarks with FP32 offering significant performance boosts over FP64. When TF32 performance is accessible it offers significant additional performance. These performance benefits are available with simple vectorized code that runs well on both CPUs and GPUs. FP16 performance via TF32 is significantly faster, but otherwise generally performs at a FP32 or FP64 level or in some cases is much slower. INT16/FP8/INT8 is more sporadically supported, when supported it generally provides FP64 level performance. We did explore what is required to achieve fully FP8 performance on the h100nvl from Python and while significant performance was achieved, the methods used are generally beyond a typical Python user. More specifically, the Stream benchmark performs well across all systems in both Matlab and Python. Good Stream performance is often very important in high level languages that tend to do more data copies during the course of execution

TABLE IV B ENCHMARK R ESULTS .

Stream results are reported in GBytes/s. All other benchmarks are reported in GFlop/s. Blurimage (FFT) results use direct convolution operations counts (see 2 Table I) to highlight the algorithmic performance gains; underlying FFT operation counts can be estimated by dividing by Nkern /(log(Nx ) log(Ny )). The following labeling convention is used: bold much better than FP64 MatMul, blue best precision, * TF32 enabled, ** beyond typical user, ≪FP64 much slower than FP64 MatMul, u unsupported.

Hardware Year 2018 2020 2024 2018 2024 2018 2020 2024 2018 2024 2026 2026 Language Matlab Matlab Matlab Matlab Matlab Python Python Python Python Python Python Python (all units Giga) Type CPU CPU CPU GPU GPU CPU CPU CPU GPU GPU CPU GPU Benchmark Precision / HW xeon-g6 xeon-p8 e9-24c v100 h100nvl xeon-g6 xeon-p8 e9-24c v100 h100nvl e9-64c b200 Stream FP64 80 89 249 967 4,257 70 65 243 976 3,962 472 5,939 MatMul FP64 2,204 2,713 2,218 14,041 103,479 2,947 2,947 2,610 10,200 94,534 7,336 71,113 MatMul FP32 4,684 5,808 4,002 27,728 97,200 5,116 5,116 5,082 27,701 *611,809 15,958 *981,428 MatMul FP16 <<FP64 <<FP64 <<FP64 u *686,129 <<FP64 <<FP64 <<FP64 115,056 892,747 <<FP64 1,183,118 MatMul INT16 u u u u u 2,360 2,972 2,606 u 96,864 7,475 71,368 MatMul FP8 u u u u u u u u u **1,752,992 u **u MatMul INT8 u u u u u 2,383 2,877 2,543 u 104,256 6,566 71,907 MatMul FP64 (complex) 2,628 2,990 2,262 13,936 104,017 2,745 2,745 2,153 13,906 101,062 6,744 73,702 MatMul FP32 (complex) 5,216 6,182 4,642 27,239 97,699 5,448 5,448 4,866 27,819 *262,673 16,793 *1,053,564 MatMul FP16 (complex) <<FP64 <<FP64 <<FP64 u *295,656 5,687 5,687 4,842 27,786 *263,004 16,992 *1,057,609 MatMul INT16 (complex) u u u u u 2,292 2,424 2,153 u 102,332 6,540 73,696 MatMul FP8 (complex) u u u u u u u u u u u u MatMul INT8 (complex) u u u u u 2,299 2,761 2,158 u 100,887 6,517 73,364 Beamformer FP64 (complex) 2,091 2,516 2,200 13,817 77,149 3,377 2,466 2,611 13,648 71,548 5,742 71,077 Beamformer FP32 (complex) 4,546 5,432 4,643 27,125 86,048 6,893 5,420 5,681 27,284 *224,593 19,050 *943,963 Beamformer FP16 (complex) <<FP64 <<FP64 <<FP64 u *250,729 6,789 5,199 5,674 27,462 *259,062 19,402 *1,031,615 Beamformer INT16 (complex) u u u u u u u u u u u u Beamformer FP8 (complex) u u u u u u u u u u u u Beamformer INT8 (complex) u u u u u u u u u u u u Blurimage (FFT) FP64 17,992 21,996 79,093 216,172 538,203 56,298 65,735 254,841 249,993 1,715,884 269,666 2,809,502 Blurimage (FFT) FP32 30,173 36,541 129,510 345,396 705,916 97,031 113,070 239,985 592,569 2,921,901 329,655 4,407,889 Blurimage (FFT) FP16 27,502 32,894 113,965 u u 76,756 100,579 167,625 572,697 2,760,761 290,920 4,116,390 Blurimage (FFT) INT16 18,261 23,194 78,058 216,017 504,129 59,318 65,046 220,130 247,449 1,606,448 283,411 2,687,933 Blurimage (FFT) FP8 u u u u u u u u u u u u Blurimage (FFT) INT8 18,625 23,020 79,222 229,887 509,873 58,339 66,124 188,527 248,064 1,622,321 189,375 2,690,303 Blurimage (FFT) FP64 (complex) 60,424 74,441 265,580 881,918 2,130,639 53,740 63,263 300,733 861,863 3,165,512 321,758 6,050,128 Blurimage (FFT) FP32 (complex) 118,680 138,015 481,290 1,350,745 2,812,656 78,020 269,803 469,185 2,092,681 5,477,044 629,101 9,621,569 Blurimage (FFT) FP16 (complex) 109,527 131,764 461,739 u u 132,663 275,068 554,294 2,096,923 5,476,084 777,272 9,631,764 Blurimage (FFT) INT16 (complex) 60,891 73,908 263,376 864,854 2,020,847 69,402 133,967 297,222 862,451 3,158,478 427,753 6,055,115 Blurimage (FFT) FP8 (complex) u u u u u u u u u u u u Blurimage (FFT) INT8 (complex) 62,341 76,883 271,335 878,376 2,032,205 55,698 141,024 289,669 862,385 3,115,789 477,472 6,055,312

2026 Python GPU rtx6000 1,770 3,031 *400,070 631,911 3,041 **u 3,044 3,034 *417,892 *420,265 3,172 u 3,155 3,114 *386,255 *408,858 u u u 368,256 1,634,442 1,597,221 366,600 u 367,418 792,188 2,762,242 2,762,608 792,331 u 792,314

that low-level languages. Matrix multiplication (MatMul) generally performs as described above. The main observation is that FP32 is offers significant performance boost and FP16 offers even greater performance when TF32 can be used. Beamformer tracks complex matrix multiplication and shows that matrix multiplication is accessible within real applications contexts. Perhaps most interesting is the significant performance achievable with Blurimage (FFT). This highlights the potential compute intensive vs memory intensive algorithmic tradeoffs now possible with very large and very fast memory. VII. C ONCLUSION A wide range of applications are currently being enabled by high level programming languages and GPU accelerators. Retaining productivity requires effective abstractions for achieving performance within a compute node (vertical scaling) and over different generations of hardware (temporal scaling). One such abstraction that enables high level programming to achieve highly scalable performance are distributed arrays. By deriving parallelism from data locality, which naturally leads to high memory bandwidth efficiency, distributed arrays provide a straightforward approach for achieving parallel performance. Using the matrix multiplication and two applications, Beamformer and Blurimage, benchmarks, distributed array performance has been explored on a variety of hardware. Direct comparison of hardware improvements in combined GPU memory bandwidth and capacity over this period shows a roughly 70× increase over nearly a decade. Over the same interval, the matrix multiplication benchmark improved by about 35×, Beamformer by about 32×, and Blurimage by about 8×. With substantial effort, we achieved matrix multiplication performance exceeding 1 PetaFlop for FP8 on an h100nvl GPU. Future work will seek to extend the distributed array programming model to horizontal scaling across computing nodes and other high level languages and other hardware platforms. ACKNOWLEDGEMENT The authors wish to acknowledge the following individuals for their contributions and support: Guillermo Morales, Bob Bond, Alan Edelman, Peter Fisher, John Goodhue, Jeff Gottschalk, Chris Hill, Joss Knight, Charles Leiserson, Kristen Malvey, Chasen Milner, Sanjeev Mohindra, Sandy Pentland, Heidi Perry, Christian Prothmann, Steve Rejto, Scott Ruppel, Daniela Rus, Pete Saunders, Matt Schuyler, Mark Sherman, Scott Van Broekhoven. R EFERENCES [1] D. J. Higham and N. J. Higham, MATLAB Guide. SIAM, 2016. [2] G. Van Rossum and F. L. Drake Jr, Python Tutorial, vol. 620. Centrum voor Wiskunde en Informatica Amsterdam, The Netherlands, 1995. [3] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, “Julia: A fresh approach to numerical computing,” SIAM review, vol. 59, no. 1, pp. 65– 98, 2017. [4] A. Reuther, P. Michaleas, M. Jones, V. Gadepally, S. Samsi, and J. Kepner, “AI and ML accelerator survey and trends,” in 2022 IEEE High Performance Extreme Computing Conference (HPEC), pp. 1–10, IEEE, 2022.

[5] A. Reuther, P. Michaleas, M. Jones, V. Gadepally, S. Samsi, and J. Kepner, “Lincoln ai computing survey (LAICS) update,” in 2023 IEEE High Performance Extreme Computing Conference (HPEC), pp. 1–7, IEEE, 2023. [6] L. Dagum and R. Menon, “OpenMP: an industry standard API for shared-memory programming,” IEEE computational science and engineering, vol. 5, no. 1, pp. 46–55, 1998. [7] B. Nichols, D. Buttlar, and J. Farrell, Pthreads Programming: A POSIX Standard for Better Multiprocessing. ” O’Reilly Media, Inc.”, 1996. [8] R. D. Blumofe, C. F. Joerg, B. C. Kuszmaul, C. E. Leiserson, K. H. Randall, and Y. Zhou, “Cilk: An efficient multithreaded runtime system,” ACM SigPlan Notices, vol. 30, no. 8, pp. 207–216, 1995. [9] E. Wang, Q. Zhang, B. Shen, G. Zhang, X. Lu, Q. Wu, Y. Wang, E. Wang, Q. Zhang, B. Shen, et al., “Intel Math Kernel Library,” HighPerformance Computing on the Intel® Xeon Phi: How to Fully Exploit MIC Architectures, pp. 167–188, 2014. [10] E. Agullo, J. Demmel, J. Dongarra, B. Hadri, J. Kurzak, J. Langou, H. Ltaief, P. Luszczek, and S. Tomov, “Numerical linear algebra on emerging architectures: The PLASMA and MAGMA projects,” in Journal of Physics: Conference Series, vol. 180, p. 012037, IOP Publishing, 2009. [11] F. Franchetti, T. M. Low, D. T. Popovici, R. M. Veras, D. G. Spampinato, J. R. Johnson, M. Püschel, J. C. Hoe, and J. M. Moura, “SPIRAL: Extreme performance portability,” Proceedings of the IEEE, vol. 106, no. 11, pp. 1935–1968, 2018. [12] J. Kepner, P. Aaltonen, D. Bader, A. Buluç, F. Franchetti, J. Gilbert, D. Hutchison, M. Kumar, A. Lumsdaine, H. Meyerhenke, S. McMillan, J. D. Owens, C. Yang, M. Zalewski, T. Mattson, and J. Moreira, “Mathematical foundations of the GraphBLAS,” in High Performance Extreme Computing Conference (HPEC), IEEE, 2016. [13] T. A. Davis, “Algorithm 1000: SuiteSparse: GraphBLAS: Graph algorithms in the language of sparse linear algebra,” ACM Transactions on Mathematical Software (TOMS), vol. 45, no. 4, pp. 1–25, 2019. [14] B. W. Kernighan and D. M. Ritchie, The C programming language. prentice-Hall, 1988. [15] B. Stroustrup, The C++ programming language. Pearson Education, 2013. [16] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, “Numerical recipes in fortran 77,” The art of scientific computing, vol. 1, 1992. [17] J. Sanders and E. Kandrot, CUDA by example: an introduction to general-purpose GPU programming. Addison-Wesley Professional, 2010. [18] A. Munshi, B. Gaster, T. G. Mattson, and D. Ginsburg, OpenCL programming guide. Pearson Education, 2011. [19] X. Cai, H. P. Langtangen, and H. Moe, “On the performance of the Python programming language for serial and parallel scientific computations,” Scientific Programming, vol. 13, no. 1, pp. 31–56, 2005. [20] N. Birkbeck, J. Levesque, and J. N. Amaral, “A dimension abstraction approach to vectorization in matlab,” in International Symposium on Code Generation and Optimization (CGO’07), pp. 115–130, IEEE, 2007. [21] J. McCalpin, “STREAM: Sustainable memory bandwidth in high performance computers, http://www.cs.virginia.edu/stream/,” 2005. [22] J. Kepner, C. Byun, L. Anderson, W. Arcand, D. Bestor, W. Bergeron, A. Bonn, D. Burrill, V. Gadepally, R. Haney, M. Houle, M. Hubbell, H. Jananthan, M. Jones, P. Luszczek, L. Milechin, G. Morales, J. Mullen, A. Prout, A. Reuther, A. Rosa, C. Yee, and P. Michaleas, “Easy acceleration with distributed arrays,” in 2025 IEEE High Performance Extreme Computing Conference (HPEC), p. 1–8, IEEE, Sep. 2025. [23] J. Kepner, Parallel MATLAB for Multicore and Multinode Computers. SIAM, 2009. [24] W. Cochran, J. Cooley, D. Favin, H. Helms, R. Kaenel, W. Lang, G. Maling, D. Nelson, C. Rader, and P. Welch, “What is the fast fourier transform?,” Proceedings of the IEEE, vol. 55, no. 10, pp. 1664–1674, 1967. [25] H. J. Nussbaumer, “The fast fourier transform,” in Fast Fourier transform and convolution algorithms, pp. 80–111, Springer, 1981. [26] X. Liu, J. Pool, S. Han, and W. J. Dally, “Efficient sparse-winograd convolutional neural networks,” arXiv preprint arXiv:1802.06367, 2018. [27] C. Byun, W. Arcand, D. Bestor, B. Bergeron, V. Gadepally, M. Houle, M. Hubbell, H. Jananthan, M. Jones, A. Klein, et al., “pPython performance study,” in 2023 IEEE High Performance Extreme Computing Conference (HPEC), pp. 1–7, IEEE, 2023.

[28] J. Lebak, J. Kepner, H. Hoffmann, and E. Rutledge, “Parallel VSIPL++: An open standard software library for high-performance parallel signal processing,” Proceedings of the IEEE, vol. 93, no. 2, pp. 313–330, 2005. [29] J. Kepner, W. Arcand, W. Bergeron, N. Bliss, R. Bond, C. Byun, G. Condon, K. Gregson, M. Hubbell, J. Kurz, et al., “Dynamic distributed dimensional data model (d4m) database and computation system,” in 2012 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 5349–5352, IEEE, 2012. [30] V. Gadepally, J. Kepner, L. Milechin, W. Arcand, D. Bestor, B. Bergeron, C. Byun, M. Hubbell, M. Houle, M. Jones, et al., “Hyperscaling internet graph analysis with d4m on the mit supercloud,” in 2018 IEEE High Performance extreme Computing Conference (HPEC), pp. 1–6, IEEE, 2018. [31] A. Reuther, J. Kepner, C. Byun, S. Samsi, W. Arcand, D. Bestor, B. Bergeron, V. Gadepally, M. Houle, M. Hubbell, M. Jones, A. Klein, L. Milechin, J. Mullen, A. Prout, A. Rosa, C. Yee, and P. Michaleas, “Interactive supercomputing on 40,000 cores for machine learning and data analysis,” in 2018 IEEE High Performance extreme Computing Conference (HPEC), pp. 1–6, Sep. 2018. [32] C. Byun, J. Kepner, W. Arcand, D. Bestor, W. Bergeron, M. Hubbell, V. Gadepally, M. Houle, M. Jones, A. Klein, L. Milechin, P. Michaleas, J. Mullen, A. Prout, A. Rosa, S. Samsi, C. Yee, and A. Reuther, “Optimizing xeon phi for interactive data analysis,” in 2019 IEEE High Performance Extreme Computing Conference (HPEC), pp. 1–6, 2019. [33] C. Byun, J. Kepner, W. Arcand, D. Bestor, B. Bergeron, V. Gadepally, M. Houle, M. Hubbell, M. Jones, A. Klein, P. Michaleas, J. Mullen, A. Prout, A. Rosa, S. Samsi, C. Yee, and A. Reuther, “Large scale parallelization using file-based communications,” in 2019 IEEE High Performance Extreme Computing Conference (HPEC), pp. 1–7, 2019.

Related documents

• PrismStencil: Transforming Stencil Computations into Matrix Multiplications#1343893
• MANOJAVAM: A Scalable, Unified FPGA Accelerator for Matrix Multiplication and Singular Value Decomposition in Principal Component Analysis#155256
• Snowmass 2013 Computing Frontier: Accelerator Science#1015716
• EvoCG: Self-Evolving Precision Control for Faster Conjugate Gradient#1527898
• GaborFlow: Hardware-Software Co-Optimization of Gabor Filter Accelerator Design#1325899