# How can I code this small piece of code in MEX ?

3 views (last 30 days)
Kate on 8 Aug 2013
My code is as follows:
clc;clear all
%----generate data----%
[X,Y] = meshgrid(-50:2:50, -50:2:50);
Z = X .* exp(-X.^2 - Y.^2);
X= X(:);
Y = Y(:);
Z =Z(:);
x = [X Y Z];
%---------------------%
%----- run computation -----%
[m, M] = size (x);
dim = M - 1;
K= zeros (m,m);
for it_dim=1:dim
tmp = x(:,it_dim+1) * ones(1,m) - ones(m,1) * x(:,it_dim+1)';
tmp = tmp .* tmp;
K = K + tmp;
end;
%---------------------------%
can anyone help me out in converting the 'FOR' loop portion of code into MEX code? I am having issues in particular with the following line in conversion to MEX C++:
tmp = x(:,it_dim+1) * ones(1,m) - ones(m,1) * x(:,it_dim+1)';
James Tursa on 8 Aug 2013
Please use the "{ } Code" button to format your code. Thanks.

Jan on 8 Aug 2013
Do you want to improve the speed only? Then try this at first:
tmp = bsxfun(@minus, x(:,it_dim+1), x(:,it_dim+1)');
From a mathematical point of view, (a - b)^2 can be written as a^2 - 2*a*b + b^2. This can reduce the number of multiplications massively, because the squaring can be performed for the vectors before inflating them by the multiplication by ONES.
Kate on 8 Aug 2013
Jan, the test data component is just a dummy for testing purposes. The 'run computation' component (i.e. the FOR loop) is what I want in MEX form. I am new to MEX programming and having difficulty especially with the indexing,etc in C++. Can you help out?

Jan on 9 Aug 2013
Edited: Jan on 9 Aug 2013
Perhaps this helps for the conversion to C:
K = zeros(m, m);
for it_dim = 1:dim
for i1 = 1:m
K(i1, i1) = K(i1, i1) + (x(i1, it_dim + 1) - x(i1, it_dim + 1)) ^ 2;
d = x(i1, it_dim + 1);
for i2 = i1 + 1:m
c = (d - x(i2, it_dim + 1)) ^ 2;
K(i1, i2) = K(i1, i2) + c;
K(i2, i1) = K(i2, i1) + c;
end
end
end
Btw., the original version is nicer, but this is faster:
Original/vector operations: Elapsed time is 0.737606 seconds.
Elemantary loops, mirroring: Elapsed time is 0.611425 seconds.