Can a parallel for loop write to a common matrix?

Question

I'm not certain what it means for the iterations of a parallel for loop to be independent. Is the following an example of two valid parallel for loops? They write and read the same matrix, but the matrix indices are unique to each iteration.

X = zeros(64);
parfor i = 1:64^2
    X(i) = i;
end
parfor i = 1:64
    X(i,:) = X(i,:) .* randn(1,64);
end

Answer 1

As far as parfor is concerned, the following three statements can be regarded as equivalent:

1) The iterations of a parfor loop must be independent.

2) No iteration of a parfor loop may depend on the outcome of any other iteration.

3) The iterations of a parfor loop must be able to be performed in any order (from @Oli)

How do these statements compare to a regular loop? In a typical loop from 1 to 8, the 4th iteration, for example, may depend on iterations 1, 2, and 3, since the software can be certain these iterations have already occurred by the time we reach iteration number 4. It must NOT depend on iterations 5, 6, 7, and 8, since the software can be certain these iterations will not have occurred.

In a parfor loop, as @Oli states, the loops may occur in any order. They may occur in the following order, for example, 7 3 4 1 2 5 8 6. Or any permutation of these 8 numbers. This implies something very important: There is no way of knowing before the fact which iteration will occur first. To see this, just chuck an fprintf('Up to iteration %d of %d\n', t, T) inside your parfor loop, where t is the loop subscript and T is the loop upper bound.

The above statement immediately implies the following conclusion: Since any iteration might occur first, it is critical that no iteration depends on the outcome of any other iteration. I'll conclude the answer with some examples:

X = ones(8, 8)
parfor n = 1:8
    X(:,n) = X(:,n) .* (3 * ones(8,1));
end

In this example, (3 * ones(8,1)) clearly does not depend on any other iteration - being constant with respect to the loop counter. Similarly X(:, n) does not depend on any iteration other than the nth. EDIT: I previously was using randn in the above example - see the discussion in comments provided by @AndrewJanke for why this was a bad idea. What about this situation:

X = ones(8, 8);
parfor n = 1:8
    X(:,n) = X(:,n) + (n + 1);
end

This is also perfectly valid. Although there is an n + 1 in the expression, this is not the same as depending on iteration number n + 1. Rather it is simply assigning the integer value of the current iteration number, plus 1, to X.

Finally, consider:

X = ones(8, 1);
parfor n = 2:8
    X(n, 1) = X(n-1, 1) + 1;
end

This would be perfectly valid in a regular loop, since iteration number n-1 will always occur before iteration n (assuming we're looping forwards). But in a parfor loop, this will cause an error, since iteration number n might occur before iteration number n-1. The lingo Matlab uses to describe the problem here is called "slicing". Imagine X to be sliced up by the loop iterations. Then in the nth iteration, you may only ever refer to the nth slice of X.

A final point, if I ever have doubts about a parfor loop, I read the section in the documentation entitled: "Parallel for loops in Matlab - overview" (sorry, can't find the corresponding webpage - unusual for Matlab documentation) It describes all the possible variable classifications inside loops, and the restrictions a parfor loop places on each classification. What I've discussed in this answer is really only the tip of the iceberg. For example, statements such as n = n + 1 are also invalid in a parfor loop, since n is the loop variable, and assignments to the loop variable are not allowed.