Can this code be generic and fast?

I am having issues making a function fast and generic. I have attached test data, a test script, and a number of code snippets explaining my attempts to find a satisfactory solution. I will explain what the code does, and I will then explain the problem. To run and time code snippet 1, open the test script and comment out code snippet 2 and 3. The test data can be downloaded from the following link;

https://drive.google.com/file/d/1yGzfuY8tMyviXGFL3dDhrQSSsAD2Gk5d/view?usp=sharing

What does the code do?

The code calculates the Euclidean distance between particles (sometimes referred to as nodes) in 1/2/3-dimensional space. Particles are connected by 'bonds' and connectivity information is stored in a bond list. The bond list consists of two columns, the first column contains the ID for node 'i' and the second column contains the ID for node 'j'. The code also stores the X, Y, Z component of every bond. This function must be called at every time step in a particle simulation and could be called up to a million times.

The Problem

I would like my code to be fast and generic, capable of taking 1D/2D/3D data without any user changes. I have attached some code snippets below and I will discuss the advantages and disadvantages of each one.

Code snippet 1: approximate execution time 0.13s

This is the fastest solution I have found by a considerable margin. This code iterates over the bond list and the X, Y, Z components of every bond are calculated and stored in separate 1D vectors. This code is fast, but it is not generic. Perhaps the best solution would be to create separate functions for 1D/2D/3D problems?

% Calculate the deformed bond components using a for loop
for kBond = 1:nBonds                           
    
   nodei = BONDLIST(kBond,1);
   nodej = BONDLIST(kBond,2);

   deformedBondComponentX(kBond) = deformedCoordinates(nodej,1) - deformedCoordinates(nodei,1); % X-component of deformed bond
   deformedBondComponentY(kBond) = deformedCoordinates(nodej,2) - deformedCoordinates(nodei,2); % Y-component of deformed bond
   deformedBondComponentZ(kBond) = deformedCoordinates(nodej,3) - deformedCoordinates(nodei,3); % Z-component of deformed bond

end


% Length of deformed bond
deformedLength = deformedBondComponentX.^2 + deformedBondComponentY.^2 + deformedBondComponentZ.^2;  % Move outside of for loop - optimises speed
deformedLength = sqrt(deformedLength);                                                               % sqrt outside of for loop - optimises speed

stretch = (deformedLength - undeformedLength) ./ deformedLength;

Code snippet 2: approximate execution time 0.50s

This solution is very slow, but it is generic. The code iterates over the bond list and a colon-operator is used to handle N-dimensional data.

% Calculate the deformed bond components using a for loop
for kBond = 1:nBonds                           
    
   nodei = BONDLIST(kBond,1);
   nodej = BONDLIST(kBond,2);

   deformedBondComponents(kBond,:) = deformedCoordinates(nodej,:) - deformedCoordinates(nodei,:); 

end

% Calculate the deformed length of every bond
deformedLength = deformedBondComponents.^2;
deformedLength = sum(deformedLength,2);   
deformedLength = sqrt(deformedLength);   


stretch = (deformedLength - undeformedLength) ./ deformedLength;

Code snippet 3: approximate execution time 0.30s

This code is completely vectorised. It is faster than code snippet 2 but it is still significantly slower than code snippet 1.

% Calculate the deformed X,Y,Z component of every bond 
deformedBondComponents = deformedCoordinates(BONDLIST(:,1),:) - deformedCoordinates(BONDLIST(:,2),:); 

% Calculate the deformed length of every bond
deformedLength = deformedBondComponents.^2;
deformedLength = sum(deformedLength,2);   
deformedLength = sqrt(deformedLength);   

stretch = (deformedLength - undeformedLength) ./ deformedLength;

Any help to make this code generic and fast would be much appreciated.

Thanks,

MH

Answer 1

[darova]

Can't understand: why completely vectorised function is slower than the first with loops? Maybe someone can explain?

Used switch statement (always remember about preallocating. Took me 0.28s and 3.62s with and without)

See attached files

[MCH]

Hi Darova,

Thanks for your reply. Yes that is one way to do it and is effectively the same as creating separate functions for 1D/2D/3D data. This function is contained within a for loop that is stepping forward in time and may be called a million times, so I would move the switch or if statement outside of the time stepping loop to prevent unecessary repetition.

Yes I am unsure why the vectorised function is slow. My initial thought is that memory access to deformedCoordinates is slow because of the scattered nature of the bond list. For example, node '1' and node '52,873' might share a bond.

I also do not completely understand why the colon-operator is so slow. There is some discussion here...

https://uk.mathworks.com/matlabcentral/answers/396757-inefficient-or-slow-use-of-colon-operator-within-a-matrix

Thanks again,

MH