Fastest way to initialize a double array of a non-standard value

Question

MATLAB provides functions for preallocation/initializing arrays with common values like 0 or 1. However, if we want the array to have some arbitrary double value, there are various ways to do it, and it's not obvious which one is preferable.

This problem is not new - it was previously discussed in places such as this blog post and in this answer. However, experience shows that software (specifically, MATLAB and its execution engine) and hardware change with time, and so it is possible that the best approach would be different on different systems. Unfortunately, the aforementioned sources do not provide benchmarking code, which might be the ultimate (and timeless) way to answer this question.

I'm looking for a benchmark I could run that would tell me the fastest approach to use on my system, considering that I might be using both "regular" double arrays and gpuArray double arrays, of various sizes.

Answer 1

function allocationBenchmark(arrSz)
if nargin < 1
  arrSz = 1000;
end

%% RAM
t = [];
disp('--------------- Allocations in RAM ---------------')
t(end+1) = timeit(@()v1(arrSz), 1);
t(end+1) = timeit(@()v2(arrSz), 1);
t(end+1) = timeit(@()v3(arrSz), 1);
t(end+1) = timeit(@()v4(arrSz), 1);
t(end+1) = timeit(@()v5(arrSz), 1);
t(end+1) = timeit(@()v6(arrSz), 1);
t(end+1) = timeit(@()v7(arrSz), 1);
t = 1E3 * t; % conversion to msec
disp(t); disp(" ");
[~,I] = min(t);
disp("Conclusion: method #" + I + " is the fastest on the CPU!"); disp(" ");

%% VRAM
if gpuDeviceCount == 0, return; end
t = [];
disp('--------------- Allocations in VRAM --------------')
t(end+1) = NaN; % impossible (?) to run v1 on the gpu
t(end+1) = gputimeit(@()v2gpu(arrSz), 1);
t(end+1) = gputimeit(@()v3gpu(arrSz), 1);
t(end+1) = gputimeit(@()v4gpu(arrSz), 1);
t(end+1) = gputimeit(@()v5gpu(arrSz), 1);
t(end+1) = gputimeit(@()v6gpu(arrSz), 1);
t(end+1) = gputimeit(@()v7gpu(arrSz), 1);
t = 1E3 * t; % conversion to msec
disp(t); disp(" ");
[~,I] = min(t);
disp("Conclusion: method #" + I + " is the fastest on the GPU!");

end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% RAM
function out = v1(M)
% Indexing on the undefined matrix with assignment:
out(1:M, 1:M) = pi;
end

function out = v2(M)
% Indexing on the target value using the `ones` function:
scalar = pi;
out = scalar(ones(M));
end

function out = v3(M)
% Using the `zeros` function with addition:
out = zeros(M, M) + pi;
end

function out = v4(M)
% Using the `repmat` function:
out = repmat(pi, [M, M]);
end

function out = v5(M)
% Using the ones function with multiplication:
out = ones(M) .* pi;
end

function out = v6(M)
% Default initialization with full assignment:
out = zeros(M);
out(:) = pi;
end

function out = v7(M)
% Using the `repelem` function:
out = repelem(pi,M,M);
end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% VRAM
function out = v2gpu(M)
scalar = gpuArray(pi);
out = scalar(gpuArray.ones(M));
end

function out = v3gpu(M)
out = gpuArray.zeros(M, M) + gpuArray(pi);
end

function out = v4gpu(M)
out = repmat(gpuArray(pi), [M, M]);
end

function out = v5gpu(M)
out = gpuArray.ones(M) .* gpuArray(pi);
end

function out = v6gpu(M)
% Default initialization with full assignment:
out = gpuArray.zeros(M);
out(:) = gpuArray(pi);
end

function out = v7gpu(M)
% Using the `repelem` function:
out = repelem(gpuArray(pi),M,M);
end

Running the above (e.g. with an input of 5000) results in the following:

--------------- Allocations in RAM ---------------
  110.4832  328.1685   48.7895   47.9652  108.8930   93.0481   47.9037

Conclusion: method #7 is the fastest on the CPU!

--------------- Allocations in VRAM --------------
       NaN   37.0322   17.9096   14.2873   17.7377   16.1386   16.6330

Conclusion: method #4 is the fastest on the GPU!

... which tells us the best (or equivalent) method to use in every case.