Issue with parallel loop

Can somebody point out to me what is wrong with the following parallel structure, and/or tell me how to fix it. I keep getting a restricted indices error, but I don't see what is wrong with these indices for MX2. ind will always be an integer and it will always be unique across any unique set of f,j,i.

q = 45;
vec = randn(3,q,q,6);
MX2 = zeros(6*(q)^2,2);

parfor f = 1:6
    for j = 1:q
        for i = 1:q 
            ind = (f-1)*q^2+(j-1)*q+i;
            V = vec(:,i,j,f)-vec(:,i+1,j+1,f);
            MX2(ind,1) = norm(V);
            V2 = vec(:,i+1,j,f)-vec(:,i,j+1,f);
            MX2(ind,2) = norm(V2);

          end
      end

end
MX2 = MX2.^2;
MX2 = max(MX2(:));

Thanks in advance.

Answer 1

[Andrew]

All,

With the help of Walter's answer I have figured out the correct way to create this parallel loop structure and will report it here in case anyone is having a similar issue. However, Matt's answer here is actually a much better solution to this problem (and is actually how I am now implementing this section of code now) so if you are having a similar issue, consider if you can adopt Matt's advice instead.

Here is the full answer. The issue had to deal with both what Walter pointed out with the indexing, and with the way MX was being used. If we instead write

MX2 = zeros(6,1);

parfor f = 1:6
    MX3 = zeros(q^2,1);
    for j = 1:q
        for i = 1:q
            ind = (j-1)*q+i;
            V = vec(:,i,j,f)-vec(:,i+1,j+1,f);
            V2 = vec(:,i+1,j,f)-vec(:,i,j+1,f);
            MX3(ind) = max([norm(V),norm(V2)]);

          end
      end
      MX2(f) = max([MX3;0]);
  end

MX2 = max(MX2)^2;

then the parallel loop works. Note that we had to define a new temporary variable here, MX3, in order to get things to work. As Matt mentions, we can create a situation without the temporary variable MX3 where everything would theoretically work; however, because of the way things are defined and indexed MATLAB's parfor checker will still throw an error.

Anyway, I hope this helps anyone who is having a similar issue with parallel processing.

Answer 2

[Walter Roberson]

            V = vec(:,i,j,f)-vec(:,i+1,j+1,f);
            V2 = vec(:,i+1,j,f)-vec(:,i,j+1,f);
            MX2(ind,:) = [norm(V), norm(V2)];

You were not using "consistent indexing" for MX2: in one place you had ind,1 and in the other place you had ind,2 . So instead you create a vector and assign the vector into the slice.

[Andrew]

Thanks for the response Walter. While this fixes one issue, it does not resolve all of the issues. Because of this I will give you a vote, but I will not accept your answer. Instead I will provide an answer that fully fixes the parallel issues here.

Answer 3

[Matt J]

You'll need to read Classification of Variables in Parfor Loops. Basically, because MX2 is both defined prior to the parfor loop and assigned values within the loop, parfor has no choice but to classify it as a "sliced variable". But there are all kinds of restrictions on how you can index sliced variables that are discussed at the link. Even though your indexing is parallelizable in theory, there are simply limits on parfor's ability to recognize this. Hence, the restrictions.

In any case, the use of parfor is not appropriate/optimal here. The operations you are performing are abundantly vectorizable,

            V = vec(:,1:end-1,1:end-1,:)-vec(:,2:end,2:end,:);
            Vnorms = sqrt(sum(reshape(V,3,[]).^2));

            V2 = vec(:,2:end,1:end-1,:)-vec(:,1:end-1,2:end,f);
            V2norms = sqrt(sum(reshape(V2,3,[]).^2));

            MX2=[Vnorms(:),V2norms(:)];

[Andrew]

Thanks for your help. I was struggling to figure out how best to vectorize this but I'm unfamiliar with the operations that are valid on n-dimensional arrays so I hadn't taken time to really try anything out. I'll give you a vote, but I won't accept your answer because it doesn't entirely answer the question asked and instead proposes another (albeit much better) way to do what I asked.

Answer 4

[Matt J]

Here's a parfor-based solution that avoids multiple nested loops, though I still advocate the vectorized approach in my earlier Answer.

q = 45;
vec = randn(3,q+1,q+1,6);
MX2 = zeros(6*(q)^2,2);

    parfor ind = 1:6*q^2

           [i,j,f]=ind2sub([q,q,6],ind);

            V = vec(:,i,j,f)-vec(:,i+1,j+1,f);
            V2 = vec(:,i+1,j,f)-vec(:,i,j+1,f);

            MX2(ind,:) = [norm(V), norm(V2)];

     end

    MX2 = MX2.^2;
    MX2 = max(MX2(:));