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

1 view (last 30 days)
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)';

Answers (2)

Jan
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.
  3 Comments
Jan
Jan on 8 Aug 2013
Edited: Jan on 8 Aug 2013
When you speak of optimizing, the main problem here is most likely not the multiplication and summation, but the expensive EXP() function. It does not get clear, if the generation of the test data is a complete dummy or essential part of the problem.
I C you would not replace the vector operations by summing up the single elements.
Kate
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?

Sign in to comment.


Jan
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.

Tags

Community Treasure Hunt

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

Start Hunting!