How to solve "parfor cannot be classified" issue?

Greetings,

I tried to convert my program with nested for-loops to parfor loop and got error "The variable z in a parfor cannot be classified". I read documents about "Troubleshoot Variables in parfor-Loops". But I don't totally understand it and even I tried different ways, the program still does not work.

Which rule I violate ( the slicing or nested for-loops )? And how to solve it, pack it as a function?

function xhat = eGram(Func,ParaVector,X0,scaling,flag) 

Tol = 1E-7;

% gramian parameters
n = ParaVector(1);
k = ParaVector(2);

% initialization for T for observability
T(:,:,1) = eye(n,n);
Cm = logspace(log10(0.01),log10(1),5) .';

e = eye(n,n);
s = length(Cm);
r = size(T,3);
xhat = zeros(n,n,r*s);

Ts = 0.1;
nu = 600;
[~,~,Y0] = Func(X0);

if scaling
    Y0_ = scale(Y0,Y0,Tol);
else
    Y0_ = Y0;
end

for l=1:r
    for m=1:s
        z = zeros(n,nu*k);
        parfor i=1:n    % !!! Here said variable z does not use
        
            % apply perturbed initial condition
            initvalue = X0 + Cm(m)*T(:,:,l)*e(:,i);
   
            [t,X,Y] = Func(initvalue);
            if scaling
                Y_ = scale(Y,Y0,Tol);
            else
                Y_ = Y;
            end
            
            if flag==0
                for iii = 1:k
                    z(i,nu*(iii-1)+1:nu*iii) = (Y_(:,iii) - Y0_(1,iii))';
                end
            else
                 for iii = 1:k
                    z(i,nu*(iii-1)+1:nu*iii) = (Y_(:,iii) - Y0_(:,iii))';
                end
            end
                
        end                
                
        chsi = z*z';
        xhat(:,:,(l-1)*s+m) = 1/(r*s*Cm(m)^2)*Ts*T(:,:,l)*chsi*T(:,:,l)';
        
        if any(isnan(xhat))
            disp(['l is ' num2str(l)]);
            disp(['m is ' num2str(m)]);
            disp(['i is ' num2str(i)]);
        end
        
    end
    
end


end

Answer 1

[Walter Roberson]

                 zt = zeros(1, nu*k);
                 for iii = 1:k
                    zt(1,nu*(iii-1)+1:nu*iii) = (Y_(:,iii) - Y0_(1,iii))';
                 end
                 z(ii, :) = zt;

You can probably vectorize. It looks like it would be something like

                 zt = Y_ - Y0_(1,:);
                 z(ii,:) = zt(:);

[Xin Shen]

Good suggestions!

But do you know which rule of parallel computation I violate?

[Walter Roberson]

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

Form of Indexing. Within the first-level of indexing for a sliced variable, exactly one indexing expression is of the form i, i+k, i-k, or k+i. The index i is the loop variable and k is a scalar integer constant or a simple (non-indexed) broadcast variable. Every other indexing expression is a positive integer constant, a simple (non-indexed) broadcast variable, a nested for-loop index variable, colon, or end.

You are using a range calculated from a nested for-loop index variable, rather than just a simple index variable.

You are safest if you use code of the form

temp = sliced_variable(parfor_index, :);
update temp
sliced_variable(parfor_index, :) = temp;

In the case (such as you have) where you are initializing the sliced variable to all 0 or a constant, you can often reduce that to

temp = zeros(1, size(sliced_variable,2)) + constant_initializer;
update temp
sliced_variable(parfor_index, :) = temp;

and if every location in temp is written to by the code without needing to read from it, then you can just:

create temp with appropriate code
sliced_variable(parfor_index, :) = temp;

which is what the vectorized code I showed does.

[Xin Shen]

Thank you! I followed Walter Roberson's suggestion and solved it.

I have read these links you sent to me, but frankly speaking, these pages are the most confusing help files of Matlab. Even I read them several times, when I write my programs based on parallel computation and I still frequently get errors. Is there any good books or websites explainning these variable classification issues clearly?

Thank you!

[Walter Roberson]

I recomend keeping it simple:

1. Try making the parfor index the last index into any variable. In particular, avoid making it the first index, as first index will give you the least efficient memory transfer.
2. For each variable that you are reading from, or which you are doing a read/modify/write for, extract all of the variable indexed by the parfor value into a temporary variable, so what you get in the temporary variable is just the parts that have to do with the parfor index
3. For each variable that you are writing into but not as a read/modify/write cycle, initialize a temporary variable that is the size of the data related to the slice indexed by the parfor index.
4. Manipulate the temporary variables as needed, including writing to them
5. For each read/modify/write variable, assign the temporary variable back over the slice indexed by the parfor index; for each variable writing into but not read/modify/write, assign the temporary variable into the output variable indexed by the parfor index.

If you have some very large variables you might need to modify this sequence. If you have particular need for distributed variables instead of sliced variables, then you might need to modify this sequence.