Engee documentation
Notebook

Performance comparison of parallel computing methods: matrix product

Introduction

Matrix product is one of the central operations in numerical methods, machine learning, signal processing, and physical modeling. Unlike element-by-element operations on vectors, here the computational complexity increases as and the bottleneck is not only the memory bandwidth, but also the efficiency of using the cache, vector registers, and optimized linear algebra libraries.

Cache locking is a method of organizing calculations in which large matrices are divided into small blocks of such size that they fit entirely into the processor's cache memory. Instead of repeatedly accessing slow RAM, the program reuses blocks that have already been loaded into the fast cache, and thus reduces data access time.

Data reuse is a property of an algorithm in which the same elements of input arrays are used repeatedly during calculations while they are still in the cache or registers. The higher the reuse, the lower the proportion of RAM accesses and the closer the performance is to the peak computing capacity of the processor.

In this paper, the performance of various methods for calculating the product of two square matrices of size is compared. from the sequential three-cycle option to automatic vectorization, tensor notation, and optimized computation.

Connecting libraries

Adding the necessary libraries:

  • Benchmarks — benchmark measurement of execution time and allocations (@btime);
  • LinearAlgebra — standard mul!, calling optimized Basic Linear Algebra Routines (BLAS);
  • Octavian — multithreaded BLAS-like matrix multiplication;
  • LoopVectorization — automatic loop vectorization (@turbo);
  • Tullio — Einstein tensor notation with automatic parallelization and vectorization.
In [ ]: 
for pkg in (["BenchmarkTools", "LinearAlgebra", "Octavian", "LoopVectorization", "Tullio"])
    try
        eval(Meta.parse("using $pkg"))
    catch
        Pkg.add(pkg)
        eval(Meta.parse("using $pkg"))
    end
end

Setting the task

Let's choose the product of two square matrices as a reference problem. and size .

Each element of the result ​ is a scalar product of -th row on -th column :

which is written in matrix form as .

Initial data

Let's create two random square matrices and size and the matrix for the result. Size let's choose so that the time of sequential implementation remains measurable, but not too long for interactive work.

In [ ]: 
N = 1000
A = randn(N, N)
B = randn(N, N)
C = similar(A);

Comparison of methods

Fully consistent calculation

Sequential implementation through three nested loops, wrapped in a function for stable typing by the JIT compiler. The order of the cycles j → i → k provides sequential access to columns B, and the battery s allows you not to create in C at each iteration of the inner loop.

In [ ]: 
function последовательное(C, A, B)
    @inbounds for j in axes(C, 2)
        for i in axes(C, 1)
            s = zero(eltype(C))
            for k in axes(A, 2)
                s += A[i, k] * B[k, j]
            end
            C[i, j] = s
        end
    end
    return C
end

@btime последовательное(C, A, B);

  3.698 s (0 allocations: 0 bytes)

The execution time was 3.698 s with zero allocations. This is the base reference point, relative to which we will consider the acceleration of the other methods.

LoopVectorization

The same three-cycle algorithm, but with a macro @turbo, which automatically applies SIMD vectorization (a method of performing calculations in which one processor instruction is applied to several data elements simultaneously), loop unwrapping, and instruction permutation.

In [ ]: 
function turbo_matmul(C, A, B)
    @turbo for j in axes(C, 2)
        for i in axes(C, 1)
            C[i, j] = zero(eltype(C))
            for k in axes(A, 2)
                C[i, j] += A[i, k] * B[k, j]
            end
        end
    end
    return C
end

@btime turbo_matmul(C, A, B);

  64.135 ms (0 allocations: 0 bytes)

The time was 64.14 ms with zero allocations, which is ≈57.7 times faster than the sequential version. However, the accumulation in C[i, j] inside the inner loop, it creates extra memory accesses compared to the accumulator in the register.

Octavian

Octavian is an implementation of matrix multiplication that uses multithreading and cache locking. It does not depend on external BLAS and in some cases shows results that are competitive with OpenBLAS, especially on small and medium-sized matrices.

In [ ]: 
@btime Octavian.matmul!(C, A, B);

  34.552 ms (0 allocations: 0 bytes)

The time is 34.55 ms with zero allocations, acceleration is ≈107 times relative to the sequential version.

Tullio

Einstein's Notation: @tullio C[i, j] = A[i, k] * B[k, j]. By default, Tullio uses multithreading and vectorization, and also tries to use cache locking and rearrange the order of loops. For a matrix product, this gives a result that is noticeably superior to pure @turbo.

In [ ]: 
@btime @tullio C[i, j] = A[i, k] * B[k, j];

  12.993 ms (59 allocations: 2.22 KiB)

The time is 12.99 ms with minimal allocations (2.22 KiB for service structures). Acceleration of ≈284.6 times relative to the sequential version.

LinearAlgebra.mul!

Function mul! from LinearAlgebra for Matrix{Float64} calls an optimized BLAS implementation. This is a multithreaded, vectorized, and cache-locked operation.

In [ ]: 
@btime LinearAlgebra.mul!(C, A, B);

  8.841 ms (0 allocations: 0 bytes)

The time is 8.841 ms with zero allocations, acceleration is ≈418.3 times relative to the sequential version. This is the best result among all measured methods on the CPU.

Final table

#[PIPE]
[PIPE][PIPE]

Location

[PIPE]

Method

[PIPE]

Time

[PIPE]

Acceleration

[PIPE][PIPE]
[PIPE][PIPE]

1

[PIPE]

LinearAlgebra.mul!

[PIPE]

8.841 ms

[PIPE]

≈418×

[PIPE][PIPE]
[PIPE][PIPE]

2

[PIPE]

Tullio

[PIPE]

12.99 ms

[PIPE]

≈285×

[PIPE][PIPE]
[PIPE][PIPE]

3

[PIPE]

Octavian

[PIPE]

34.55 ms

[PIPE]

≈107×

[PIPE][PIPE]
[PIPE][PIPE]

4

[PIPE]

LoopVectorization

[PIPE]

64.14 ms

[PIPE]

≈58×

[PIPE][PIPE]
[PIPE][PIPE]

5

[PIPE]

Consistent

[PIPE]

3698 ms

[PIPE]

1×

[PIPE][PIPE]
[PIPE]#

The performance spread between the sequential implementation and the optimized BLAS is about 418 times.

Conclusion

The measurements show that for the matrix product, the choice of the method gives a performance spread of almost three orders of magnitude: from 3.7 seconds of the sequential three-cycle version to 8.8 milliseconds of the standard optimized function. LinearAlgebra.mul!. At the same time, the key conclusion is that simple parallelization or autovectorization is not enough here. Consistent @turbo speeds up execution by about 58 times, but still loses. LinearAlgebra.mul! by almost an order of magnitude, because it does not perform cache locking and data reuse.

Octavian and Tullio occupy an intermediate position. Octavian shows a result close to LinearAlgebra.mul!, but on this configuration it is slightly inferior to it; Tullio, thanks to the rearrangement of cycles and multithreading, turns out to be faster than Octavian and only 1.5 times slower. mul!. This makes Tullio a reasonable choice if for some reason you cannot use the standard linear algebra library, but you need performance close to optimal.

Practical conclusion: for a matrix product, in most cases, the most effective tools will be LinearAlgebra.mul!, which automatically invokes a highly optimized linear algebra library.