25 views (last 30 days)

Show older comments

I am generating a meshgrid to be able to calculate my result fast:

% x, y, z are some large vectors

[a,b,c] = meshgrid(x,y,z);

% s, t are constants, M some matrix

result = (((c*s - b*t).^2)./(a.^2 + b.^2 + c.^2)).*M;

This is actually working quite nicely. Unfortunately, for very large x,y,z, the meshgrid function is running out of memory.

How do I rewrite the meshgrid function to be memory efficient?

I had thought of three loops like this:

result = zeros(length(x), length(y), length(z));

for i = 1:lenght(x)-1

for j = y = 1:lenght(y)-1

for k = z = 1:lenght(z)-1

b = ??

c = ??

result(i,j,k) = (((c*s - b*t).^2)./(x(i)^2 + y(j)^2 + z(k).^2));

end

end

end

result = result.*M;

What are the values for b and c?

How can I turn the outer for into a parfor?

Fabio Freschi
on 11 Jun 2020

This is how to make the three-loop version analogous to the meshgrid version

% some dummy values

N = 300;

x = linspace(1,10,N);

y = linspace(1,10,N);

z = linspace(1,10,N);

s = 1;

t = 1;

M = rand(N,N,N);

%% meshgrid

tic

% x, y, z are some large vectors

[a,b,c] = meshgrid(x,y,z);

% s, t are constants, M some matrix

result = (((c*s - b*t).^2)./(a.^2 + b.^2 + c.^2)).*M;

toc

%% three-loops

tic

% preallocation

result2 = zeros(length(x), length(y), length(z));

% note the order of the for-loop indices to mimic meshgrid

for iz = 1:length(x)

for jx = 1:length(y)

for ky = 1:length(z)

result2(iz,jx,ky) = M(iz,jx,ky)*(((z(iz)*s - y(ky)*t).^2)./(x(jx)^2 + y(ky)^2 + z(iz).^2));

end

end

end

toc

% check results

norm(result(:)-result2(:))./norm(result(:))

However I don't see how you can avoid running out of memory: meshgrid is creating a N*N*N (with my notation) matrix, if it runs out of memory, also the preallocation of the result matrix will

result2 = zeros(length(x), length(y), length(z));

It is however true that in the second version you only have 2 matrices with dimensions N*N*N (M and result2) whereas in the first case you have 5 (a, b, c, M, result).

Note that according to my tests, the meshgrid version with vectorization is always faster than the version with three loops

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!