Vectorizing nested for loops

I have looked at other examples of vectorizing for loops, but I don't know how to apply those solutions to what I have.

Essentially, in this part of the code, I am randomly generating xyz coordinates for particles, and storing them in a coordinate array. I then want to calculate the distance, d, between particles. If this distance is less than a certain number based on the radius, R, and a constant, T, then I want to consider those particles in contact. I have a spanMatrix set up where a 1 indicates two particles are connected.

I am generating hundreds to thousands of particles, and this part of my code becomes very slow. Could anyone please help me speed this up? Thank you!

% Place an entry into a coordinate array with the random xyz coordinates
    coordArray(numParticles,1) = x(1);
    coordArray(numParticles,2) = y(1);
    coordArray(numParticles,3) = z(1);
    
    % For every element in the coordinate array
    for i = 1:size(coordArray)
        for j = 1:i
            % for every (i,j) pair do distance calculation
            d=(sqrt((coordArray(i,1)-coordArray(j,1)).^2 + (coordArray(i,2)-coordArray(j,2)).^2 + (coordArray(i,3)-coordArray(j,3)).^2));
            
            % then check if the distance is close enough
            if d <= (2*R + T)
                % then these two span
                spanMatrix(i,j) = 1;
                spanMatrix(j,i) = 1;
            end % end if d
        end % end for j
    end % end for i

Answer 1

[Guillaume]

I am generating hundreds to thousands of particles

That's not many, it's only when you get in the order of a hundred thousand that the following starts to become problematic. The following generates a full symmetric matrix:

%coordArray: a Nx3 matrix
%generates a NxN full matrix of distance (using euclidean distance)
pointdist = sqrt(sum((permute(coordArray, [1, 3, 2]) - permute(coordArray, [3, 1, 2])) .^2, 3));
spanMatrix = pointdist <= 2*R + T;

A thousand points would only require around 8 MB of memory to store each full matrix pointdist and spanMatrix. Even 10,000 points is only around 800 MB. 100,000 points on the other hand requires around 75 GB. If you get to that number of points, then you can't vectorise the whole but can certainly vectorise the inner loop:

%coordArray: a Nx3 matrix
coordpairs = [];  %to store coordinate pairs that will be used to construct the final sparse matrix
for ptidx = 1:size(coordpairs, 1)
    pointdist = sqrt(sum((coordArray - coordArray(ptidx, :)) .^2, 2);
    isclose = find(pointdist <= 2*R + T);
    coordpairs = [coordpairs; [repmat(ptidx, numel(isclose), 1), isclose]];  %#ok<AGROW>
end
spanMatrix = sparse(coordpairs(:, 1), coordpairs(:, 2), 1)

Finally, note that in your original code:

for i = 1:size(coordArray)

works but most likely only by accident. Note that the size function returns a vector and that the colon operator only use the first element of the matrix. It's better to always be explicit:

for i = 1:size(coordArray, 1)

Answer 2

[darova]

Use pdist2() to calculate combinations of distances

D = pdist2( ... )
cond = d <= (2*R + T);        % indices of matrix you want
spanMatrix(cond) = 1;

[darova]

Like my answer please