Performance comparison of parallel computing methods: element-wise operations on vectors

Introduction

The rapid growth of data volumes and the increasing complexity of computational models in physics, engineering, and data sciences make the problem of improving the performance of software code one of the central ones in modern numerical research. Reducing computing time is inevitably associated with the proper use of parallelism at all levels — from vector instructions of a single processor core to involving multiple GPU cores.

This example provides a comparison of the performance of parallel and vectorized computing methods designed for efficient array processing. An element-by-element operation on two vectors of large dimension, including a combination of algebraic and transcendental functions, is chosen as a reference problem.

Importing libraries

We will attach the necessary libraries.

Benchmarks — benchmark measurement of execution time and allocations (@btime);
Random — generation of initial random data.

In this example, we will measure the calculation time using the following methods:

CUDA — performing calculations on a GPU;
Flops — flexible parallel loops with transparent multithreading;
Folds — parallel piecemeal operations;
LoopVectorization — automatic vectorization of loops;
Polyester — low-level multithreading with minimal overhead;
ThreadsX — parallel processing of collections with automatic load balancing;
Tullio is Einstein tensor notation with automatic parallelization and vectorization.

import Pkg
Pkg.add(["BenchmarkTools", "CUDA", "FLoops", "Folds", "LoopVectorization", "Polyester", "Random", "ThreadsX", "Tullio"])
using BenchmarkTools, CUDA, FLoops, Folds, LoopVectorization, Polyester, Random, ThreadsX, Tullio

Initial data

To test the performance of piecemeal operations on vectors, we use this function.

This function was chosen as a representative sample of real scientific and engineering computing.: it combines easy arithmetic (multiplication) and transcendental functions of varying complexity ( and ) and the addition of three heterogeneous components. This allows you to evaluate memory bandwidth, vectorization quality, the ability to merge cycles without creating intermediate arrays, multithreaded load balancing on uneven tasks, PROCESSOR instruction scheduling efficiency, and GPU scalability in a single test. This combination of factors will allow you to objectively evaluate the performance of each of the presented methods.

Create random values of vectors and .

N = 10^7
a = randn(N)
b = randn(N)

Sequential and vectorized calculations

Sequential calculations

# Consistent
c = similar(a);
@btime begin
    @inbounds for i in eachindex(a)
        c[i] = a[i] * b[i] + sqrt(abs(a[i])) + sin(b[i])
    end
end

  5.218 s (119998978 allocations: 1.94 GiB)

The execution time of the sequential calculation cycle was 5.218 seconds. We use this indicator as a base for comparison with other methods.

Vectorized calculations

# Vectorized
c = similar(a);
@btime @. c = a * b + sqrt(abs(a)) + sin(b);

  235.731 ms (5 allocations: 128 bytes)

Using point translation (@.) allowed to delegate calculations to an optimized function broadcast, avoiding global scope issues. The execution time is 235.7 ms, which is ~22.1 times faster than sequential calculations.

Sequential calculations within a function

# Consistent within the function
c = similar(a);
function serial(c, a, b)
    @inbounds for i in eachindex(a)
        c[i] = a[i] * b[i] + sqrt(abs(a[i])) + sin(b[i])
    end
    return c
end

@btime serial(c, a, b);

  186.282 ms (0 allocations: 0 bytes)

Wrapping the loop in a function allowed the JIT compiler to derive stable types and generate native code. Time is 186.3 ms, acceleration is ~28.0 times.

Parallel computing

Polyester

# Polyester
c = similar(a);
function polyester!(c, a, b)
    @batch for i in eachindex(a)
        c[i] = a[i] * b[i] + sqrt(abs(a[i])) + sin(b[i])
    end
    return c
end

@btime polyester!($c, $a, $b);

  53.342 ms (0 allocations: 0 bytes)

The macro @batch efficiently distributed loop iterations between CPU threads with extremely low overhead. The execution time is 53.34 ms, which is in ~97.8 times faster than sequential calculations.

Floops

# Floops
c = similar(a);
function floops!(c, a, b)
    @floop for i in eachindex(a)
        c[i] = a[i] * b[i] + sqrt(abs(a[i])) + sin(b[i])
    end
    return c
end

@btime floops!($c, $a, $b);

  47.111 ms (35 allocations: 2.62 KiB)

Thanks to transparent multithreading and the ability to broadcast actions, the calculation execution time was 47.11 ms, which is ~110.7 times faster than sequential calculations.

Folds

# Folds
c = similar(a)
@btime Folds.map!(i -> $a[i]*$b[i] + sqrt(abs($a[i])) + sin($b[i]), $c, eachindex($a, $b));

  46.395 ms (42 allocations: 3.28 KiB)

Parallel convolution in a functional style without creating intermediate arrays: Folds.map! effectively distributes the load. Time is 46.40 ms, acceleration is ~112.5 times.

ThreadsX

# ThreadsX
c = similar(a)
@btime ThreadsX.map!(i -> $a[i]*$b[i] + sqrt(abs($a[i])) + sin($b[i]), $c, eachindex($a, $b));

  44.999 ms (253 allocations: 18.84 KiB)

Specialized parallel processing of collections provided the best balance among pure multithreaded solutions. Time is 45.00 ms, acceleration is ~116.0 times.

LoopVectorization

# LoopVectorization
c = similar(a);
function turbo!(c, a, b)
    @turbo for i in eachindex(a)
        c[i] = a[i] * b[i] + sqrt(abs(a[i])) + sin(b[i])
    end
    return c
end

@btime turbo!($c, $a, $b);

  32.258 ms (0 allocations: 0 bytes)

Warning: detected a stack overflow; program state may be corrupted, so further execution might be unreliable.

The execution time was 32.26 ms, which is ~161.8 times faster than a sequential calculation cycle.

Tullio

# Tullio
function tullio!(c, a, b)
    @tullio c[i] = a[i] * b[i] + sqrt(abs(a[i])) + sin(b[i])
end
@btime tullio!($c, $a, $b);

  8.288 ms (41 allocations: 1.67 KiB)

Einstein's notation with automatic parallelization and vectorization, being wrapped in a function, showed impressive results. The execution time is 8.288 ms with minimal allocations (1.67 KiB), which is ~629.6 times faster than sequential calculations.

CUDA

To use the CUDA library.jl, it is necessary to synchronize the version of the runtime and the CUDA driver. You may need to run this command and restart the Julia kernel.

CUDA.set_runtime_version!(CUDA.driver_version())

Let's do the calculations using CUDA.jl.

# CUDA
da = CuArray(a)
db = CuArray(b)
dc = similar(da)

function cuda_broadcast!(dc, da, db)
    @. dc = da * db + sqrt(abs(da)) + sin(db)
    return dc
end

cuda_broadcast!(dc, da, db)
CUDA.synchronize()

@btime begin
    cuda_broadcast!($dc, $da, $db)
    CUDA.synchronize()
end

  198.631 μs (101 allocations: 2.61 KiB)

The arrays were transferred to the graphics card's memory, and calculations were performed on GPU cores. The time is ~0.199 ms (198.6 microseconds), which is ~26,200 times faster than sequential calculations. An absolute measurement record.

Final table

Location	Method	Time	Acceleration
1	CUDA	0.199 ms	~26 200×
2	Tullio	10.04 ms	~520×
3	LoopVectorization	32.26 ms	~162×
4	ThreadsX	45.00 ms	~116×
5	Folds	46.40 ms	~113×
6	FLoops	47.11 ms	~111×
7	Polyester	53.34 ms	~98×
8	Consistent in function	186.3 ms	~28×
9	Vectorized	235.7 ms	~22×
10	Consistent	5218 ms	1×

Choosing the right method gives a 26,000-fold performance spread: from 5.2 seconds of sequential cycle to 199 microseconds on the GPU.

Conclusion

The measurements showed that different ways of performing the same task can give completely different operating times — from a few seconds to fractions of a millisecond. Even the simple design of calculations in the form of a function gives a noticeable acceleration. Special libraries for multithreaded calculations allow you to get even more benefits, and using a video card makes calculations almost instantaneous.

A simple practical conclusion follows from this: you should not stop at the first working version of the code. By trying several available tools, you can significantly reduce the calculation time, which means you can get results faster and solve larger tasks.