How can I use vectorization to speed up the following code

2 次查看（过去 30 天）
Oluwaseun Lijoka
Oluwaseun Lijoka 2015-11-1
评论： Oluwaseun Lijoka 2015-11-1
I have the following code. x, y and t are variables. M is a big number and I discovered that the double loop will take a longer time. I need help in making it faster, all other example have seen are partially helpful. Thanks a lot.
function f=exact(x,y,t,M,j)
f=0;
for i=1:M
for j=1:M
f=f+(16*(sin(j*pi)*sin(i*pi)*l*i^2*pi^3+4*sin(j*pi)*cos(i*pi)*j*i*pi^2+2*sin(i*pi)*cos(j*pi)*i^2*pi^2+2*sin(j*pi)*j*i*pi^2-2*i^2*pi^2*sin(i*pi)...
-6*sin(j*pi)*sin(i*pi)*j*pi+8*cos(i*pi)*cos(j*pi)*i*pi-8*i*pi*cos(i*pi)+4*cos(j*pi)*i*pi-12*sin(i*pi)*cos(j*pi)-4*i*pi+12*sin(j*pi)))/(i^4*pi^7*j^3)...
*cos(sqrt(i^2+j^2)*pi*t).*sin(i*pi*x).*sin(j*pi*y);
end
end
  2 个评论
Jan
Jan 2015-11-1
编辑：Jan 2015-11-1
The expected size of M matters: Completely different approaches are required for M=4 and M=4*1e7 . "Faster" ist not clear enough: Do you have a fixed time limit, such that it is worth to spend an hour to create a C-Mex function? Or doy you want to save some minutes only, such that spending an hour would mean a loss of total time?
The input "j" is not used, because it is overwritten by the FOR loop counter. So perhaps this is a typo: The variable "l" is undefinde. Then accelerating buggy code is fruitless.
Jan
Jan 2015-11-1
编辑：Jan 2015-11-1
The formula looks very strange: You evaluate a lot of "sin(i*pi)" and "cos(j*pi)", but "i" and "j" are intergers. Then rounding errors dominate the result: Mathematically "sin(i*pi)" is exactly zeros, but the limited precision leads to a value of 1.2246e-016. Nevertheless, you can omit all terms, which contain a "sin(i*pi)" or "sin(j*pi)", while the cos-expressions can be replaced by a multiplication by 1 - which can be omitted.
I have the strong impression, that you did not post the real code. But optimizing pseudo-code is meaningless. Please provide a valid input also.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Jan
Jan 2015-11-1
编辑：Jan 2015-11-1
Although I'm in doubt that the shown code is correct, I've cleaned it up a little bit moving repeated calculations outside the inner loop:
function f=exact(x,y,t,M,l) % "l" instead of "j"?!
pi3 = pi^3;
pi2 = pi^2;
pi7 = pi^7;
f = 0;
for i=1:M
i2 = i^2;
sin_i_pi = sin(i * pi);
cos_i_pi = cos(i * pi);
sin_i_pi_x = sin(i * pi * x);
C1 = sin_i_pi * i2;
C2 = 2 * i2 * pi2 * sin_i_pi;
C3 = 8 * i * pi * cos_i_pi;
C4 = 4 * i * pi;
C5 = i^4 * pi7;
C6 = 4 * cos_i_pi * i * pi2;
C7 = 8 * cos_i_pi * i * pi;
C8 = 6 * sin_i_pi * pi;
for j=1:M
sin_j_pi = sin(j * pi);
cos_j_pi = cos(j * pi);
f = f + (sin_j_pi * C1 * l * pi3 + ...
sin_j_pi * j * C6 + ...
2 * C1 * cos_j_pi * pi2 + ...
2 * sin_j_pi * j * i * pi2 - ...
C2 - ...
sin_j_pi * j * C8 + ...
cos_j_pi * C7 - ...
C3 + ...
4 * cos_j_pi * i * pi - ...
12* sin_i_pi * cos_j_pi - ...
C4 + 12 * sin_j_pi) / ...
(C5 * j^3) * ...
cos(sqrt(i2 + j^2) * pi * t) .* sin_i_pi_x .* sin(j*pi*y);
end
end
f = f * 16;
This takes 1.10 sec instead of 3.15 sec on my old Matlab R2011b/64/Win7/Core2Duo.
But as written in my comment already: All sin terms can be omitted, because they are 0, while the cos expression is 1, such that the formula can be simplified substantially.
  3 个评论
显示 1更早的评论
Jan
Jan 2015-11-1
The "l" is undefined in the first term:
sin(j*pi) * sin(i*pi) * l * i^2 * pi^3
What about the vanishing value of "sin(i*pi)"?
Oluwaseun Lijoka
Oluwaseun Lijoka 2015-11-1
Dear Jan, The code is an exact solution to a solved wave equation in which when M is tending to infinity, that is why I stated M to be a big number. I want to compare this exact solution with a numerical approximation in terms of error and convergence. Looping i and j over M each time could be waste of time. Here is the whole code:
function f=exact(x,y,t,M,j)
if (nargin<5) j=0; end f=0; % tic if (j==0) % exact soln for m=1:M for l=1:M f=f+((16*(sin(l*pi)*sin(m*pi)*l*m^2*pi^3+4*sin(l*pi)*cos(m*pi)*l*m*pi^2+2*sin(m*pi)*cos(l*pi)*m^2*pi^2+2*sin(l*pi)*l*m*pi^2-2*m^2*pi^2*sin(m*pi)... -6*sin(l*pi)*sin(m*pi)*l*pi+8*cos(m*pi)*cos(l*pi)*m*pi-8*m*pi*cos(m*pi)+4*cos(l*pi)*m*pi-12*sin(m*pi)*cos(l*pi)-4*m*pi+12*sin(m*pi)))/(m^4*pi^7*l^3))... .*cos(sqrt(m^2+l^2)*pi*t).*sin(m*pi*x).*sin(l*pi*y);
end
end
elseif (j==1) % derivative w.r.t time
for m=1:M
for l=1:M
f=f+((16*(sin(l*pi)*sin(m*pi)*l*m^2*pi^3+4*sin(l*pi)*cos(m*pi)*l*m*pi^2+2*sin(m*pi)*cos(l*pi)*m^2*pi^2+2*sin(l*pi)*l*m*pi^2-2*m^2*pi^2*sin(m*pi)...
-6*sin(l*pi)*sin(m*pi)*l*pi+8*cos(m*pi)*cos(l*pi)*m*pi-8*m*pi*cos(m*pi)+4*cos(l*pi)*m*pi-12*sin(m*pi)*cos(l*pi)-4*m*pi+12*sin(m*pi)))/(m^4*pi^7*l^3))...
*(-sqrt(m^2+l^2)*pi).*sin(sqrt(m^2+l^2)*pi*t).*sin(m*pi*x).*sin(l*pi*y);
end
end
elseif (j==2) % derivative w.r.t variable x for m=1:M for l=1:M f=f+((16*(sin(l*pi)*sin(m*pi)*l*m^2*pi^3+4*sin(l*pi)*cos(m*pi)*l*m*pi^2+2*sin(m*pi)*cos(l*pi)*m^2*pi^2+2*sin(l*pi)*l*m*pi^2-2*m^2*pi^2*sin(m*pi)... -6*sin(l*pi)*sin(m*pi)*l*pi+8*cos(m*pi)*cos(l*pi)*m*pi-8*m*pi*cos(m*pi)+4*cos(l*pi)*m*pi-12*sin(m*pi)*cos(l*pi)-4*m*pi+12*sin(m*pi)))/(m^4*pi^7*l^3))... *(m*pi).*cos(sqrt(m^2+l^2)*pi*t).*cos(m*pi*x).*sin(l*pi*y); end
end
elseif (j==3) % derivative w r t variable y
for m=1:M
for l=1:M
f=f+((16*(sin(l*pi)*sin(m*pi)*l*m^2*pi^3+4*sin(l*pi)*cos(m*pi)*l*m*pi^2+2*sin(m*pi)*cos(l*pi)*m^2*pi^2+2*sin(l*pi)*l*m*pi^2-2*m^2*pi^2*sin(m*pi)...
-6*sin(l*pi)*sin(m*pi)*l*pi+8*cos(m*pi)*cos(l*pi)*m*pi-8*m*pi*cos(m*pi)+4*cos(l*pi)*m*pi-12*sin(m*pi)*cos(l*pi)-4*m*pi+12*sin(m*pi)))/(m^4*pi^7*l^3))...
*(l*pi).*cos(sqrt(m^2+l^2)*pi*t).*sin(m*pi*x).*cos(l*pi*y);
end
end
end Thanks a lot

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by