Memory efficient alternative for meshgrid?

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?

[Stephen Cobeldick]

https://www.mathworks.com/matlabcentral/fileexchange/74956-ndgridvecs

Answer 1

[Fabio Freschi]

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