MATLAB

Loop ordering

Subcategory: Core

In Fortran, MATLAB, and R, arrays are stored in column-major order, meaning that the entries from the left-most index are stored contiguously in memory. This is often referred to as the left-most index varying quickest. This is different from many other popular scientific programming languages such as C, C++, and Python, which use a row-major ordering. The implication for writing loops is that it is always best to order loops by which indices vary quickest. In this way, the loop nest will traverse contiguous memory, which can be done efficiently.

In Fortran, MATLAB, and R, arrays are stored in column-major order, meaning that the entries from the left-most index are stored contiguously in memory. This is often referred to as the left-most index varying quickest. This is different from many other popular scientific programming languages such as C, C++, and Python, which use a row-major ordering. The implication for writing loops is that it is always best to order loops by which indices vary quickest. In this way, the loop nest will traverse contiguous memory, which can be done efficiently.

Example Code

Consider the following Fortran example involving a loop over a 3D array, which uses the recommended column-major ordering:

integer, parameter :: n = 500
real, dimension(n, n, n) :: a, b, c
integer :: i, j, k

! Initialisation of a and b omitted

! Hand-coded c = a + b using column-major
do k = 1, n
  do j = 1, n
    do i = 1, n
      c(i,j,k) = a(i,j,k) + b(i,j,k)
    end do
  end do
end do

Without compiler optimisations, this will likely run several times faster than the (not recommended) equivalent row-major version:

integer, parameter :: n = 500
real, dimension(n, n, n) :: a, b, c
integer :: i, j, k

! Initialisation of a and b omitted

! Hand-coded c = a + b using row-major
do i = 1, n
  do j = 1, n
    do k = 1, n
      c(i,j,k) = a(i,j,k) + b(i,j,k)
    end do
  end do
end do

With compiler optimisations, simple loops such as this where the result of the sum are unused would however likely be optimised out by the compiler. In practice you wouldn’t normally have a redundant loop in your code, so this shouldn’t be a concern.

A further comment on the simple loops above is that the same could be achieved with an array operation, which has a more compact notation:

integer, parameter :: n = 500
real, dimension(n, n, n) :: a, b, c

! Initialisation of a and b omitted

! Compute c = a + b with an array operation
c(:,:,:) = a + b

Benchmark

Below these examples have been extended into a full benchmark, to demonstrate the impact. An additional checksum is used to ensure the code we’re benchmarking isn’t optimised away by the compiler.

Benchmark Code

```fortran program benchmark_loops implicit none integer, parameter :: n = 1000 real, allocatable :: a(:,:,:), b(:,:,:), c(:,:,:) real :: t1, t2 integer :: i real :: checksum allocate(a(n,n,n), b(n,n,n), c(n,n,n)) ! Initialize data a = 1.0 b = 2.0 c = 0.0 ! Warm-up runs (optional but recommended) call col_major_add(a, b, c, n) call row_major_add(a, b, c, n) ! Benchmark column-major loop order call cpu_time(t1) call col_major_add(a, b, c, n) call cpu_time(t2) print *, "Column-major time (s): ", t2 - t1 checksum = sum(c) print *," checksum: ", checksum ! Reset result array to avoid reuse artifacts c = 0.0 ! Benchmark row-major loop order call cpu_time(t1) call row_major_add(a, b, c, n) call cpu_time(t2) print *, "Row-major time (s): ", t2 - t1 checksum = sum(c) print *," checksum: ", checksum contains subroutine col_major_add(a, b, c, n) integer, intent(in) :: n real, intent(in) :: a(n,n,n), b(n,n,n) real, intent(out) :: c(n,n,n) integer :: i, j, k do k = 1, n do j = 1, n do i = 1, n c(i,j,k) = a(i,j,k) + b(i,j,k) end do end do end do end subroutine col_major_add subroutine row_major_add(a, b, c, n) integer, intent(in) :: n real, intent(in) :: a(n,n,n), b(n,n,n) real, intent(out) :: c(n,n,n) integer :: i, j, k do i = 1, n do j = 1, n do k = 1, n c(i,j,k) = a(i,j,k) + b(i,j,k) end do end do end do end subroutine row_major_add end program benchmark_loops ```

Compiled with GFortran and compiler optimisations enabled, column-major accesses were 2.4x faster under Ubuntu!

$ gfortran -O3 benchmark.f90
$ ./a.out
 Column-major time (s):    3.98375320    
  checksum:    67108864.0    
 Row-major time (s):       9.72134018    
  checksum:    67108864.0

Under WSL on different hardware a greater 7.1x speedup was seen, so results are likely to vary, but should remain positive.

$ gfortran -O3 benchmark.f90
$ ./a.out
 Column-major time (s):   0.440809250
  checksum:    67108864.0
 Row-major time (s):       3.11280060
  checksum:    67108864.0

Technical Details

The underlying reason that order has such a big impact here is due to how computer memory operates. When a processor loads a variable from RAM into it’s caches, it doesn’t load an individual variable (likely 4 bytes), it loads a full cache line (likely 64 bytes).

Therefore, variables stored consecutively in memory will be loaded at the same time if they fall within the same cache line. This means the processor does not need to go to RAM to access the variable, which is orders of magnitude higher latency. Due to how Fortran lays out memory, column-major ordering achieves this.

In contrast, if you operate as though memory is row-major, consecutively accessed variables will be stored a long way apart from each other in memory. Hence each individual access will need to load a full cache line from RAM. Due to this high turnover, it’s likely by the time a second variable would be accessed from any cache line, that the cache line has been evicted to be replaced with a fresher load.

If you’re working with 4-byte types you could be doing 16x the RAM accesses, with 8-byte types it’s still 8x!

Furthermore, in order for the compiler to perform vectorisation, it needs to be able to apply vector instructions to a full cache line of variables. If you are not operating consecutively on variables within a cache line, it can be much harder for the compiler to detect that vectorisation is appropriate.

Preallocation

Subcategory: Core

When you grow an array dynamically inside a for loop, MATLAB must repeatedly allocate new memory and copy the existing contents. This can significantly degrade performance. Preallocating reserves sufficient contiguous memory upfront, avoiding costly resizing operations later.

When you grow an array dynamically inside a for loop, MATLAB must repeatedly allocate new memory and copy the existing contents. This can significantly degrade performance. Preallocating reserves sufficient contiguous memory upfront, avoiding costly resizing operations later.

Example Code

To take advantage of Preallocation, you should allocate your array outside of the for-loop before using it.

The example below creates an array that grows with every iteration of the loop

tic
x = 0;
N = 10000000;
for k = 2:N
   x(k) = x(k-1) + 5;
end
toc

Elapsed time is 0.380283 seconds.

If we allocate the array before entering the loop, the code runs almost 10x faster

tic
x = zeros(1,N);
N = 10000000;
for k = 2:N
   x(k) = x(k-1) + 5;
end
toc

Elapsed time is 0.040541 seconds.

The Technical Detail

Prior to MATLAB R2010b, memory allocation worked like this:

The first time through the loop, the variable x hasn’t been assigned to yet, so MATLAB creates it through the assignment to x(1). The second time through the loop assigns to x(2), which doesn’t exist yet. So MATLAB allocates enough new memory for two elements and then copies x(1) and assigns x(2) into the new space. The third time through the loop assigns to x(3), which doesn’t exist yet. So MATLAB allocates enough new memory for three elements and then copies x(1) and x(2) and assigns x(3) into the new space.

See the pattern? On the k-th time through the loop, MATLAB has to make a copy of k elements into newly allocated memory. The time required to copy k elements is proportional to k. The time required to execute the entire loop is therefore proportional to the sum of the integers from 1 to n, which is n*(n+1)/2.

Bottom line: memory copying during the assignment statement causes the time required to execute the loop to be proportional to n^2. For large N, this is very slow!

From R2010b onwards, MATLAB uses a much more efficient memory allocation strategy when it encounters loops like this so you are not ‘punished’ as much as you used to be if your code follows this pattern. However, it is still significantly less efficient than simply creating the array once and writing to it in the loop.

Vectorisation

Subcategory: Core

MATLAB is optimised to take advantage of vectorisation for mathematical operations between arrays, whereby the processor executes one instruction across multiple variables simultaneously. Vectorisation can perform mathematical operations many times faster than a for loop and even take advantage of conditional logic.

MATLAB is optimised to take advantage of vectorisation for mathematical operations between arrays, whereby the processor executes one instruction across multiple variables simultaneously. Vectorisation can perform mathematical operations many times faster than a for loop and even take advantage of conditional logic.

Example Code

To take advantage of vectorisation within MATLAB, operations should be performed at an array-level, rather than by looping over the individual elements.

The below example multiplies every element of an array by 2 using a for loop.

% Allocate an array of 1 million elements
vals = 1:1000000;
% Multiply all elements by 2
for ii=1:length(vals)
    vals(ii) = vals(ii) * 2;
end

To use vectorisation, it can be rewritten as a single line of code

% Allocate an array of 1 million elements
vals = 1:1000000;
% Multiply all elements by 2
vals = vals * 2;

Not all of your MATLAB code will be this simple though, maybe you have conditional logic which restricts which elements the operation should be performed on.

% Allocate an array of 1 million elements
vals = 1:1000000;
% Multiple by 2 all elements divisible by 3, divide other elements by 2
for ii=1:length(vals)
    if mod(ii, 3) == 0
        vals(ii) = vals(ii) * 2;
    else
        vals(ii) = vals(ii) / 2;
    end
end

The above example can be written as

% Allocate an array of 1 million elements
vals = 1:1000000;
% Produce a boolean array k of whether elements are divisible by 3
k = mod(vals, 3) == 0
% Conditionally divide or multiply elements by 2
vals(k) = vals(k) * 2;
vals(~k) = vals(~k) / 2;

Using vectorisation is both more concise and more performant.

The Technical Detail

Vector instructions, which MATLAB can take advantage of, enable a CPU to apply the same operation to multiple data elements in parallel with a single thread.

Modern CPUs use SIMD (Single Instruction, Multiple Data) instructions to operate on several values packed into a register. Since a typical CPU cache line is 64 bytes, and standard data types like 32-bit or 64-bit floats and integers are 4 or 8 bytes each, 8–16 values can fit within a single cache line and be processed together by a single vector instruction.

However, to take advantage of this, the data must be aligned in memory. Laid out such that it fits neatly into the expected cache line boundaries. MATLAB arrays are explicitly designed for numerical performance, they allocate memory with alignment guarantees, ensuring that vector instructions can be used. Therefore, you just need to make sure where possible you’re using operations which support vectorisation.

GPU Arrays

Subcategory: Parallel Computing Toolbox

MATLAB makes it easy to run matrix operations on an NVIDIA GPU. There is a cost for moving data to and from the GPU, but if your entire algorithm can be done in matrix or vector operations, and your data is big enough, it can be faster.