How to rewrite parfor loop using arrayfun to compute on GPU?

Is there a way to rewrite this parfor/for loops to compute on GPU using arrayfun? And will there be any difference?

parfor i = 1:size(A,1)

    for j = i+1:size(A,1)

        for k = j+1:size(A,1)

            sumMatrix = [sumMatrix; [(A(i) + A(j) + A(k)) A(i) A(j) A(k)]];

        end
    end
end

A is a vector of integers.

Answer 1

[Joss Knight]

Seems like you're just trying to add together all the elements of A in every permutation, with the proviso that k > j > i. So really you shouldn't be using loops at all. Just form the list of numbers you want to add and prune out all the ones you don't want.

[Ak, Aj, Ai] = ndgrid(A);
[maskK, maskJ, maskI] = ndgrid(1:numel(A));
mask = maskJ > maskI & maskK > maskJ;
Acat = [Ai(mask) Aj(mask) Ak(mask)];
sumMatrix2 = [sum(Acat,2) Acat];

In answer to your last question, yes, it makes an astronomical difference. With even 100 numbers the looping method takes several minutes while the vectorized method takes a fraction of a second. Although you could speed up the looping method a massive amount by preallocating the output.

[Stepan Ulyanin]

Hi! Your version works perfectly on arrays of moderate sizes, but I tried to run the code on an array with 1000000 elements and it resulted in a 100% utilization of memory and crush afterwards, do you have any ideas?

I also changed the code a little since there are some bounding conditions, so I check for a condition now and if it works then I add to the array, if it doesn't I disregard the entry. My CPU runs it fairly fast but when I run for 1 000 000 000 000 numbers it takes a long time, and I wondered if there is more efficient algorithm or if use of GPU would make the computation faster?

Here is the code I am using right now:

parfor i = 1:(N-2)

    %draw the progress bar
    parfor_progress;

    %sum int at i with int at i+1 and with int at i+2
    for j = i+1:(N-1)

        for k = j+1:N

            %sum of 3 perfect squares
            sum = arrayOfPerfectSquares(i) + arrayOfPerfectSquares(j) + arrayOfPerfectSquares(k);

            %sumMatrix = [sumMatrix; [(A(i) + A(j) + A(k)) A(i) A(j) A(k)]];

            % ---------
            %if at least one of the perfect squares * 3 are equal to the sum ->
            %append the sum and the squares to the matrix
            if (3 * arrayOfPerfectSquares(i) == sum || 3 * arrayOfPerfectSquares(j) == sum || 3 * arrayOfPerfectSquares(k) == sum) && (arrayOfPerfectSquares(i) + arrayOfPerfectSquares(j) == 2 * arrayOfPerfectSquares(k) || arrayOfPerfectSquares(i) + arrayOfPerfectSquares(k) == arrayOfPerfectSquares(j) || arrayOfPerfectSquares(j) + arrayOfPerfectSquares(k) == 2 * arrayOfPerfectSquares(i))

                sumMatrix = [sumMatrix; [sum arrayOfPerfectSquares(i) arrayOfPerfectSquares(j) arrayOfPerfectSquares(k)]];

            end

            %----------

        end
    end
end

I am sorry but I am fairly new to MATLAB. Thanks for your response!