Parfor Loop Error with Classification

When trying to convert the outer loop of the code below to a parfor loop, the error is:

Error: The variable z_new in a parfor cannot be classified.
See Parallel for Loops in MATLAB, "Overview" 

Following the instruction on another post, I used a temporary variable to represent z_new, but that did not work either. Can anyone see what needs to be changed to permit this to run with parfor?

x = linspace(0,6,100); y = x; [xx,yy] = meshgrid(x,y); 
z = sin(xx + 0.25*yy);

dx = 1; dy = 1; Dep = 100*ones(size(z));
z_new = z + Dep;

[nr,nc] = size(z);

nnx = round(Dep/dx,0); nny = round(Dep/dy,0);


parfor i = 1:nc
    for j = 1:nr
        
        for ii = i-nnx(j,i):i+nnx(j,i)
            for jj = j-nny(j,i):j+nny(j,i)
                iinew = ii;
                jjnew = jj;
                
                Thk = Dep.*Dep - dx*dx*(i-ii)*(i-ii) - dy*dy*(j-jj)*(j-jj);
                
                if (ii > 0 && ii <= nc && jj > 0 && jj <= nr && Thk(j,i) >= 0)
                    Th_tmp = z(j,i) + sqrt(Thk(j,i));
                    
                    if (iinew > 0 && iinew <= nc && jjnew > 0 && jjnew <= nr)
                    if z_new(jjnew,iinew) < Th_tmp
                        z_new(jjnew,iinew) = Th_tmp;
                    end
                    end
                end
            end
        end
    end
end

Answer 1

[Edric Ellis]

I think you can adapt this to work with parfor, but it's going to be a little bit of a stretch. Your variable z_new can be considered to be a parfor reduction variable. Normally, simple "reduction variables" are fine for computing a maximum value. Consider the simplest case I can think of:

largest = -Inf;
parfor i = 1:10
    newValue = rand();
    largest = max(newValue, largest);
end

This looks slightly different to your example in a couple of ways: firstly, there's no indexing into largest - we'll deal with that later; secondly - I've used max rather than if to convince parfor that largest is indeed a reduction variable.

So, next we need to deal with the indexing. This is a bit more complicated. We can do this like so:

largest = -Inf(1, 3); % Starting point
parfor i = 1:10
    % Allocate an update for this iteration of the parfor loop
    update = -Inf(size(largest));
    for j = 1:3
        % Fill out an element of 'update'
        update(j) = rand();
    end
    % Finally, apply the update all at once
    largest = max(largest, update);
end

The trick here is to create a single update array that we can use with a single "reduction" operation max. In your case, you'll need to allocate your equivalent of update just ahead of the loops over ii and jj and then apply it just after those loops.

[Paul Safier]

@Edric Ellis :Thanks for the suggestions! I've been wrestling with the code to reformulate it in a way that seems consistent with your example. However, I have yet to be successful in satisfying parfor's requirements--this time the error indicates a temporary variable "largest" is intended as a reduction variable...

Error: The variable 'largest' is perhaps intended as a reduction variable, but is actually an uninitialized temporary. For more information

Here is the updated code. I also simplified it a tad. It runs properly when the parfor is just for.

Can you see still a way to make largest acceptable to parfor?

x = linspace(0,6,10); y = x; [xx,yy] = meshgrid(x,y); 
z = sin(xx + yy);

dx = 1; dy = 1; Dep = 10;
z_new = z + Dep;

[nr,nc] = size(z);

nnx = Dep/dx; nny = Dep/dy;

 largest = z + Dep;
% Reminder: rows=j=y. cols=i=x.

parfor i = 1:nc
for j = 1:nr
              
    update = -Inf(size(largest));
        
        for ii = i-nnx:i+nnx
        for jj = j-nny:j+nny
                iinew = ii;
                jjnew = jj;
                
                
                Thk = Dep*Dep - dx*dx*(i-ii)*(i-ii) - dy*dy*(j-jj)*(j-jj); 
                
                
                 if (ii > 0 && ii <= nc && jj > 0 && jj <= nr && Thk >= 0)
                    Th_tmp = z(j,i) + sqrt(Thk);
                    update(jj,ii) = Th_tmp;
                 end

                
        end
        end
            largest = max(largest,update);
end

end

[Edric Ellis]

That's really close to working, but unfortunately, the rules for reduction variables are very stringent - you cannot access the reduction variable in any way at all other than the reduction expression. It took me a while to spot it, but it's the innocuous size(largest) that's tripping you over. Fortunately, in this case, it's simple to "hoist" that expression outside the loop like this:

szLargest = size(largest)
parfor ...
    update = -Inf(szLargest);
end

[Paul Safier]

@Edric Ellis : Thanks so much!! That worked. What a benign hang-up! Thank you for your help, I appreciate it!

[Paul Safier]

@Edric Ellis : on a side note, I converted the double loop to a single and that seems to work, and obviously MUCH faster:

ntot = nr*nc;
parfor k = 1:ntot
    
    % Determine i and j
        i = fix(k/nr);
    if k/nr - fix(k/nr) ~= 0
        i = i + 1;
    end
        j = nr - (i*nr - k);
        
        stuff...
end