Can this code be generic and fast?

2 次查看(过去 30 天)
MCH
MCH 2019-4-24
评论: MCH 2019-4-24
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;
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

回答(1 个)

darova
darova 2019-4-24
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
  1 个评论
MCH
MCH 2019-4-24
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...
Thanks again,
MH

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 MATLAB Compiler 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by