Keep a temporary variable after parfor

I am just wondering how do I keep my value z after the parfor has run?

% x (real(c))
tic;
ticBytes(gcp);
dx=0.01; % step
x1=-2;
x2=2;

% y (imag(c))
dy=0.01; % step
y1=-1.5;
y2=1.5;
x=x1:dx:x2;
y=y1:dy:y2;
[X,Y] = meshgrid(x,y); % for vectorized calculation
c=X+1i*Y;
R=1000; % if abs(z)>R then z is considered as infinity
n=1000; % maximal number of iterations, 100 is close to infinity


parfor nc=1:n
    z=zeros(size(c)); % starts from zeros
    z=z.^2+c; % vectorized
end

I=abs(z)<R; % logical image
imagesc(x,y,I);
colormap((jet));
set(gca,'XColor', 'none','YColor','none');
toc;
tocBytes(gcp);

[Rik]

The content of the loop doesn't depend on the loop variable and is overwritten every iteration. Since parfor cannot guarantee an order, z is not clearly defined (even though it actually is).

You should first make your code work with for before you swap it with parfor.

[Oisin Kelly]

The code does work as a regular For loop, I had z=zeros(size(c)); above the for loop but to get it to work in the parfor i moved it inside the loop. Parfor loop works when i comment out the plotting parts of the code. But i appreciate your reply.

[Rik]

The problem with your code is that this:

for nc=1:n
    z=zeros(size(c)); % starts from zeros
    z=z.^2+c; % vectorized
end

has the same result as this:

z=c;

Why do you want to use a parfor loop? One main requirement of a parallel loop is that every iteration is independent.

[Raymond Norris]

If I understand you correctly, your code use to look like this

z=zeros(size(c)); % starts from zeros
for nc=1:n
    z=z.^2+c; % vectorized
end

In order to speed it up, you rewrote the for-loop to be parfor. However, z is a temporary value (in the loop) and must be defined on the left hand size before referencing it on the right hand side. You then pushed the zeros call into the parfor, but z (a temporary variable) doesn't get passed back and therefore can't be plotted. Do I have that right?

We could probably go into the intracacies of sliced and reducation variables, but I want to make sure this is the actual code we're working with. If so, you ought not need to the loop to beging with, no? Won't this do:

% x (real(c))
tic;
dx=0.01; % step
x1=-2;
x2=2;

% y (imag(c))
dy=0.01; % step
y1=-1.5;
y2=1.5;
x=x1:dx:x2;
y=y1:dy:y2;
[X,Y] = meshgrid(x,y); % for vectorized calculation
c=X+1i*Y;
R=1000; % if abs(z)>R then z is considered as infinity
n=1000; % maximal number of iterations, 100 is close to infinity

z=z.^2+c; % vectorized
    
I=abs(z)<R; % logical image
imagesc(x,y,I);
colormap((jet));
set(gca,'XColor', 'none','YColor','none');
toc

[Oisin Kelly]

Thank you for your reply but I tried the code you suggested and it just plots a red square. I need the for loop to complete the iterations of the Mandelbrot set as it needs to filter out the values that are within and those that are not.

z=zeros(size(c)); % starts from zeros

for nc=1:n
    z=z.^2+c; % vectorized
end 

You're correct when you say this is the code that I originally used before the parfor , I am trying to use all the cores in my system to see the difference it makes but what im starting to understand is the code I have based on Pauls answer below it would be next to impossible without doing a complete overhaul on my code?

[Paul Hoffrichter]

@Raymond Norris

>> go into the intracacies of sliced and reducation variables

I have an idea where you are going with this, and would appreciate whether you think it would work in the parfor loop. Given the square of s within the parfor loop, I did not think that s could be broken into chunks for each of the parfor workers since the workers would not know how to initialize s properly in each of the workers' processes. I try to explain why I cam to this conclusion in my posts below.

[Rik]

Why do you want to remove your question?

(Just in case you decide to be rude and edit away your question: I have made a capture so I can trivially restore that.)

Answer 1

[Paul Hoffrichter]

>> Parfor loop works when i comment out the plotting parts of the code.

I tried what you said and only experienced what you wrote when I forgot to clear the z-matrix after the previous run. Otherwise, it did not work. Were you leaving around your prevoius z?

Your loop appears to be trying to accumulate a value across an iteration of the loop, regardless of the iteration order. Using Matlab terminology, you wanted your z-matrix to be a reduction variable, rather than a temporary variable.

"Reduction Variables: Variables that accumulates a value across iterations of the loop, regardless of iteration order"

Ref: https://www.mathworks.com/help/parallel-computing/troubleshoot-variables-in-parfor-loops.html

"Another common cause of uninitialized temporaries can arise when you have a variable that you intended to be a reduction variable. However, if you use it elsewhere in the loop, then it is classified as a temporary variable."

Ref: https://www.mathworks.com/help/parallel-computing/temporary-variable.html

Consider these two loops:

s = 0;
constant = 1;
for ii = 1:100
   s = s + constant; % i.e., s(n) = s(n-1) + constant;
end
>> s
s =
   100

s = 0;
constant = 1;
parfor ii = 1:100
   s = s + constant; % i.e., s(n) = s(n-1) + constant;
end
>> s
s =
   100

No Matlab errors. Despite the current value of s in the loop being dependent upon the previous value, the Matlab parfor analyzer could see that each parfor worker could initialize its value of s to 0, and then add together the results after all the workers have accumulated their sums. The variable s in this case is called a "reduction" variable.

But in your case, you have:

   s = s.^2 + constant;

Because of the squaring, only one worker would know how to initialize the value of s. To initialize the other workers' s-values, a worker would have to compute it knowing its ii-chunk. But that defeats the purpose of parfor.

To see the issue:

s1 = s0.^2 + c;
s2 = s1.^2 + c = (s0.^2 + c).^2 + c = s0.^4 + 2*c*s0.^2 + c.^2 + c; % looking hard to compute for another worker
s3 = s2.^2 + c = (s0.^4 + 2*c*s0.^2 + c.^2 + c).^2 + c; % a parfor worker would give up on initial conditions for s.

[Oisin Kelly]

First I would like to thank you for you're reply and for such detail you went into to help me.

I am not too experienced when it comes to coding let alone doing all the intracate details. I wanted to use parfor to use all my 4 cores on my system but starting to see that isn't possible with the code I have from reading your answer. I did not realize it was keeping a previous z value that is a mistake on my behalf. I did notice that my code ran twice as fast when it did work but that could be luck im guessing based on your reply.

From what I read there you're saying its being used as a reduction value until it is being called in the next line of code below it which causes it to be a temporary value? Interesting.

I am learning a lot today some of which is very overwhelming , I understand having different tasks for each cores but just unsure how to divide them up. Here was me thinking it would just be a simple as adding in a small line or clicking a button in matlab that enabled all the cores to be used. I have a lot to learn but again thank you for your reply.

[Paul Hoffrichter]

The conclusion is that you cannot have

s = s.^2 + constant;

in the parfor loop, since the parfor analyzer is unable for each of its workers to efficiently compute initial conditions for s due to the squaring of s inside the parfor loop. So, s, instead of being a reduction variable, is just a temporary variable.

[Paul Hoffrichter]

@Oisin Kelly,

Sorry that I did not read you reply before I wrote down my second post with my conclusion.

>> you're saying its being used as a reduction value until it is being called in the next line of code

What I was trying to say is that it appeared to me that you were hoping that z would be a reduction variable, in which case, the parfor would work. But, in fact, it never achieved that status because of the z-squaring. This squaring makes the separate independent workers impossible to initialize s in their own worker's process.

Here is one way, in the future, to try to manually figure out whether a parfor can work. See if you can break down a for loop in just, say, 2 chunks (e.g., 1:N/2 and N/2+1:N), and the see how you should initialize the variables whos values are needed after the two workers have completed their portion of the parfor-loop. Then, at the parfor end-line, can you figure out how to put the two results together to form the desired answer. Even in the simple case where

z = z + constant; % i.e., z(n) = z(n-1) + constant;

it is not immediately obvious (to me) that it can be done, but with a pencil and paper, I am able to prove it possible. Or, you can keep this as a referece for the classifications of variables inside a parfor-loop:

Loop Variables       Loop indices
Sliced Variables     Arrays whose segments are operated on by different iterations of the loop
Broadcast Variables  Variables defined before the loop whose value is required inside the loop, but never assigned inside the loop
Reduction Variables  Variables that accumulates a value across iterations of the loop, regardless of iteration order
Temporary Variables  Variables created inside the loop, and not accessed outside the loop

https://www.mathworks.com/help/parallel-computing/troubleshoot-variables-in-parfor-loops.html#bvjfik_

One last note. Your line of code in the parfor loop is that of a difference equation. If you can figure out how to make z a closed function of (x,y), then you will not need any loop at all.

[Raymond Norris]

Reduction variables are broken into two parts. Let's take the following example

z = 10;
constant = 5;
parpool(4);
parfor idx = 1:400
    z = z + constant;
end

parfor will send n chunks to each of the 4 workers. Let's assume that first chunk is n=50. Each worker is going to add constant to the addition identity matrix (i.e. 0) n times.

worker 1

local_z = 0;    % Because 0 is the addition identity matrix
for ii=1:50
    local_z = local_z+constant
end

worker 2

local_z = 0;
for ii=1:50
    local_z = local_z+constant
end

ditto for worker 3 and worker 4.

Once a worker is done, it'll get its next chunk, let's say n=25

local_z = 0;
for ii=1:25
    local_z = local_z+constant
end

Eventually, all 400 iterations are complete and each worker has a local summation. Let's assume each worker has done the same, then each local sum is 500. Then parfor does a final summation of the initial value (10) and each of the workers: 10 + 500 + 500 + 500 + 500 = 2010.

This is similar for

z = z * constant;

where we use the multiplication identity matrix (1).