How to convert two nested for-loops to one parfor loop.

I have the following code. I want to get it in one parfor loop.

clear; clc;

% number of points:
Nx = 80;
Ny = 90;



xs = linspace(-2*pi/(3),4*pi/(3),Nx);
ys = linspace(-2*pi/(sqrt(3)),2*pi/(sqrt(3)),Ny);


% Allocate memory
ZZ = zeros(Nx,Ny,8);
XX = zeros(Nx,Ny); 
YY = zeros(Nx,Ny);


for ix = 1:Nx
    x = xs(ix);

    for iy = 1:Ny
        y = ys(iy);

        FUN = fun(x,y);

        % sort eigenvalues:
        [~,D] = eig(FUN); 
        [D,I] = sort(diag(real(D)),'descend');
        evals = diag(D);

        % store data:
        ZZ(ix,iy,:) = diag(evals);
        XX(ix,iy) = x;
        YY(ix,iy) = y;

    end
end

%% Plot figure
figure;
tiled = tiledlayout(2,2,"TileSpacing","tight","Padding","compact");
for q = 1:2:8
    nexttile
    surf(XX,YY,ZZ(:,:,q),'LineStyle','none','FaceColor','interp');
    view(2)
    axis([-2*pi/(3) 4*pi/(3) -2*pi/(sqrt(3)) +2*pi/(sqrt(3))])
    box on
    grid off
    axis square
    colorbar
end




%%
function out = fun(x,y)

out = [                                                                    3,                                                                     0, - 2*cos(x/4 + (3^(1/2)*y)/4) - (3*cos((3*x)/4 - (3^(1/2)*y)/4))/5,                                                                     0, - 2*cos(x/4 - (3^(1/2)*y)/4) - (3*cos((3*x)/4 + (3^(1/2)*y)/4))/5,                                                                     0,                             - 2*cos(x/2) - (3*cos((3^(1/2)*y)/2))/5,                                                                     0
    0,                                                                    -3,                                                                     0, - 2*cos(x/4 + (3^(1/2)*y)/4) - (3*cos((3*x)/4 - (3^(1/2)*y)/4))/5,                                                                     0, - 2*cos(x/4 - (3^(1/2)*y)/4) - (3*cos((3*x)/4 + (3^(1/2)*y)/4))/5,                                                                     0,                             - 2*cos(x/2) - (3*cos((3^(1/2)*y)/2))/5
    - 2*cos(x/4 + (3^(1/2)*y)/4) - (3*cos((3*x)/4 - (3^(1/2)*y)/4))/5,                                                                     0,                                                                    -1,                                                            -2*2^(1/2),                             - 2*cos(x/2) - (3*cos((3^(1/2)*y)/2))/5,                                                                     0, - 2*cos(x/4 - (3^(1/2)*y)/4) - (3*cos((3*x)/4 + (3^(1/2)*y)/4))/5,                                                                     0
    0, - 2*cos(x/4 + (3^(1/2)*y)/4) - (3*cos((3*x)/4 - (3^(1/2)*y)/4))/5,                                                            -2*2^(1/2),                                                                     1,                                                                     0,                             - 2*cos(x/2) - (3*cos((3^(1/2)*y)/2))/5,                                                                     0, - 2*cos(x/4 - (3^(1/2)*y)/4) - (3*cos((3*x)/4 + (3^(1/2)*y)/4))/5
    - 2*cos(x/4 - (3^(1/2)*y)/4) - (3*cos((3*x)/4 + (3^(1/2)*y)/4))/5,                                                                     0,                             - 2*cos(x/2) - (3*cos((3^(1/2)*y)/2))/5,                                                                     0,                                                                    -1,                                                  2^(1/2) - 6^(1/2)*1i, - 2*cos(x/4 + (3^(1/2)*y)/4) - (3*cos((3*x)/4 - (3^(1/2)*y)/4))/5,                                                                     0
    0, - 2*cos(x/4 - (3^(1/2)*y)/4) - (3*cos((3*x)/4 + (3^(1/2)*y)/4))/5,                                                                     0,                             - 2*cos(x/2) - (3*cos((3^(1/2)*y)/2))/5,                                                  2^(1/2) + 6^(1/2)*1i,                                                                     1,                                                                     0, - 2*cos(x/4 + (3^(1/2)*y)/4) - (3*cos((3*x)/4 - (3^(1/2)*y)/4))/5
    - 2*cos(x/2) - (3*cos((3^(1/2)*y)/2))/5,                                                                     0, - 2*cos(x/4 - (3^(1/2)*y)/4) - (3*cos((3*x)/4 + (3^(1/2)*y)/4))/5,                                                                     0, - 2*cos(x/4 + (3^(1/2)*y)/4) - (3*cos((3*x)/4 - (3^(1/2)*y)/4))/5,                                                                     0,                                                                    -1,                                                  2^(1/2) + 6^(1/2)*1i
    0,                             - 2*cos(x/2) - (3*cos((3^(1/2)*y)/2))/5,                                                                     0, - 2*cos(x/4 - (3^(1/2)*y)/4) - (3*cos((3*x)/4 + (3^(1/2)*y)/4))/5,                                                                     0, - 2*cos(x/4 + (3^(1/2)*y)/4) - (3*cos((3*x)/4 - (3^(1/2)*y)/4))/5,                                                  2^(1/2) - 6^(1/2)*1i,                                                                     1
    ];

end

[Matt J]

It is not clear why you expect fun(x,y) to be an 8x8 matrix. As posted, fun() returns a 1x12 vector.

Answer 1

[Matt J]

It does not seem advisable to use parfor. Everything in your code is vectorizable. However, here is what a parfor approach could look like:

[XX,YY,ZZ]=meshgrid(xs,ys);  
[M,N]=size(XX);

ZZ=zeros(M*N,8);

parfor i=1:M*N
 
        ZZ(i,:)=sort( eig(  fun(XX(i),YY(i))  ),'descend');    
    
end

ZZ=reshape(ZZ,[M,N,8]);