MATLAB running slow on MacBook pro

I have been running a MATLAB program on MacBook Pro for almost six hours now, and it is still not complete. It is cycling through three while loops (the outer two loops are n=855, the inner loop is n=500). Is this a surprise that it is taking this long? Is there anything I can do to increase the speed? I am including the code below, as well as the variable data types underneath that.

while i < (numAtoms + 1)
    pointAccessible = ones(numPoints,1);
    j = 1;
    while j <(numAtoms + 1)
        if (i ~= j)
            k=1;
            while k < (numPoints + 1)
                if (pointAccessible(k) == 1)
                    sphereCoord = [cell2mat(atomX(i)) + p + sphereX(k), cell2mat(atomY(i)) + p + sphereY(k), cell2mat(atomZ(i)) + p + sphereZ(k)];
                    neighborCoord = [cell2mat(atomX(j)), cell2mat(atomY(j)), cell2mat(atomZ(j))];
                    coords(1,:) = [sphereCoord];
                    coords(2,:) = [neighborCoord];
                    if (pdist(coords) < (atomRadius(j) + p))
                        pointAccessible(k)=0;
                    end
                end
                k = k + 1;
            end
        end
        j = j+1;
    end
    remainingPoints(i) = sum(pointAccessible);
    i = i +1;
end

Variable Data Types:
numAtoms = 855
numPoints = 500
p = 1.4
atomRadius = <855 * 1 double>
pointAccessible = <500 * 1 double>
atomX, atomY, atomZ = <1 * 855 cell>
sphereX, sphereY, sphereZ = <500 * 1 double>
remainingPoints = <855 * 1 double>

Answer 1

[Ken Atwell]

Your code's performance is suffering from repeated error-checking code in pdist. This is rather hacky and something of an advanced maneuver, but starting with Andreas's code, if you do the following, execution wall time will drop by 75% or so. I did this with R2013b, other versions of MATLAB may differ in the specifics.

1. Copy the file pdistmex.mexmaci64 to wherever your working folder is (where your script is stored). This file can be found in the MATLAB distribution under "matlabroot", in the subfolder toolbox/stats/stats/private. You need to do this so you can call this function directly.
2. Change this code:

 pdist(coords)

to:

 pdistmex(coords', 'euc', [])

If this seems like magic, single-step through the code and into the implementation of pdist so understand where this is coming from.

If you need further performance gains, this code is fundamentally parallel-friendly. Change your outmost 'while' loop to a 'for' loops, and construct 'coords' a little differently and you can make the outer loop a 'parfor' loop using Parallel Computing Toolbox.

 coords = [1 + p + sphereX(k), 1 + p + sphereY(k), 1 + p + sphereZ(k); 1, 1, 1];

[Andreas Goser]

"fundamentally parallel-friendly". Nice phrase :-)

[Ken Atwell]

This variant on your code, using Parallel Computing Toolbox, ran in 20 minutes on my 8-core Mac Pro.

 clear all
numAtoms = 855;
numPoints = 500;
p = 1.4;
atomRadius = rand(855,1);
pointAccessible = rand(500,1);
sphereX = rand(500,1);
sphereY = rand(500,1);
sphereZ = rand(500,1);
remainingPoints = rand(855,1);
tic
parfor i=1:numAtoms
    pointAccessible = ones(numPoints,1);
    j = 1;
    while j <(numAtoms + 1)
        if (i ~= j)
            k=1;
            while k < (numPoints + 1)
                if (pointAccessible(k) == 1)
                    coords = [ 1 + p + sphereX(k), 1 + p + sphereY(k), 1 + p + sphereZ(k); 1, 1, 1]
                    if (pdistmex(coords', 'euc', []) < (atomRadius(j) + p))
                        pointAccessible(k)=0;
                    end
                end
                k = k + 1;
            end
        end
        j = j+1;
    end
    remainingPoints(i) = sum(pointAccessible);
end
toc

Answer 2

[Andreas Goser]

The following code on my machine needs 18.5 s for each outer loop. So it will need 4.5 hours (Win7 64 bit, R2013a). Please compare a test and lets see if it is an issue with your environment or the code. We can go from there.

numAtoms = 855
numPoints = 500
p = 1.4
atomRadius = rand(855,1);
pointAccessible = rand(500,1);
%atomX, atomY, atomZ = <1 * 855 cell>
sphereX = rand(500,1);
sphereY = rand(500,1);
sphereZ = rand(500,1);
remainingPoints = rand(855,1);
i=1
tic
while i < (numAtoms + 1)
    i
    pointAccessible = ones(numPoints,1);
    j = 1;
    while j <(numAtoms + 1)
        if (i ~= j)
            k=1;
            while k < (numPoints + 1)
                if (pointAccessible(k) == 1)
                    sphereCoord = [1 + p + sphereX(k), 1 + p + sphereY(k), 1 + p + sphereZ(k)];
                    neighborCoord = [1, 1, 1];
                    coords(1,:) = [sphereCoord];
                    coords(2,:) = [neighborCoord];
                    if (pdist(coords) < (atomRadius(j) + p))
                        pointAccessible(k)=0;
                    end
                end
                k = k + 1;
            end
        end
        j = j+1;
    end
    sum(pointAccessible)
    remainingPoints(i) = sum(pointAccessible);
    i = i +1;
    toc
end
toc

[Lindz]

Hello:

I ran it yesterday starting at 3:30pm. I went to bed at 1am, it was not complete. When I woke up at 8am, it was complete. So, it takes anywhere from 9.5-16.5 hours.

Thank you...