Help regarding parallelizing code

for i=1:(length(water)/3)-1
for j=i+1:length(water)/3
donor_mol_index=[i j]; %between 2 water molecules either can be donor/acceptor
for p=donor_mol_index
donor_h_index=[1 3];% water atom data 1st and last atoms are hydrogen atoms
for k=donor_h_index
hb_ww_count=hb_ww_count+1;
%di-donor index, ai-acceptor index, donor is the molecule which gives the h-atom
if p==i
di=i; ai=j; %if donor molecule index equals i
else
di=j;ai=i;%if donor molecule index equals j
end
%dist_ww(cnt_steps,hb_ww_count)=sqrt(sum((water_mol(di).atom_data(k,3:5)-water_mol(ai).atom_data(2,3:5)).^2));
v1=water_mol(di).atom_data(k,3:5)-water_mol(di).atom_data(2,3:5);%vector connecting donor hydrogen-donor oxygen
v2=water_mol(di).atom_data(k,3:5)-water_mol(ai).atom_data(2,3:5);%vector connecting donor hydrogen-acceptor oxygen
angle=acos(dot(v1, v2) / (norm(v1) * norm(v2)));
angle_ww(cnt_steps,hb_ww_count)=radtodeg(angle);
if sum((water_mol(di).atom_data(k,3:5)-water_mol(ai).atom_data(2,3:5)).^2)<=dist_hb^2 %distance between donor hydrogen and acceptor oxygen
v1=water_mol(di).atom_data(k,3:5)-water_mol(di).atom_data(2,3:5);%vector connecting donor hydrogen-donor oxygen
v2=water_mol(di).atom_data(k,3:5)-water_mol(ai).atom_data(2,3:5);%vector connecting donor hydrogen-acceptor oxygen
angle=acos(dot(v1, v2) / (norm(v1) * norm(v2)));
angle=radtodeg(angle);
if angle<=angle_hb
hb_ww(cnt_steps,hb_ww_count)=1;
end
else
hb_ww(cnt_steps,hb_ww_count)=0;
end
end
end
end
end
I would like to parallelize the outer loop. Can I put all the code inside the first for loop in a function and supply the loop variable i and the structures to the function ? hb_ww_count is a reduction variable. I use this variable for indexing another array within my code.
Please provide some suggestions. I am sort of lost.

2 个评论

Start with simplifying your code: why do you calculate v1 and v2 and the angle twice? Change:
v1=water_mol(di).atom_data(k,3:5)-water_mol(di).atom_data(2,3:5);%vector connecting donor hydrogen-donor oxygen
v2=water_mol(di).atom_data(k,3:5)-water_mol(ai).atom_data(2,3:5);%vector connecting donor hydrogen-acceptor oxygen
angle=acos(dot(v1, v2) / (norm(v1) * norm(v2)));
angle_ww(cnt_steps,hb_ww_count)=radtodeg(angle);
if sum((water_mol(di).atom_data(k,3:5)-water_mol(ai).atom_data(2,3:5)).^2)<=dist_hb^2 %distance between donor hydrogen and acceptor oxygen
v1=water_mol(di).atom_data(k,3:5)-water_mol(di).atom_data(2,3:5);%vector connecting donor hydrogen-donor oxygen
v2=water_mol(di).atom_data(k,3:5)-water_mol(ai).atom_data(2,3:5);%vector connecting donor hydrogen-acceptor oxygen
angle=acos(dot(v1, v2) / (norm(v1) * norm(v2)));
to:
% Before the loops:
dist_hb2 = dist_hb ^ 2; % Once only
len = floor(length(water)/3);
n = 4 * (len - 1) * len / 2;
hb_ww = zeros(cnt_steps, n); % Pre-allocate!!!
angle_ww = zeros(cnt_steps, n);
...
v1 = water_mol(di).atom_data(k,3:5) - water_mol(di).atom_data(2,3:5);%vector connecting donor hydrogen-donor oxygen
v2 = water_mol(di).atom_data(k,3:5) - water_mol(ai).atom_data(2,3:5);%vector connecting donor hydrogen-acceptor oxygen
angle = acos(dot(v1, v2) / (norm(v1) * norm(v2))); % Not accurate, prefer atan2
angle_ww(cnt_steps, hb_ww_count) = angle * (180 / pi);
if sum(v2 .^ 2) <= dist_hb2 %distance between donor hydrogen and acceptor oxygen
if angle<=angle_hb
hb_ww(cnt_steps, hb_ww_count)=1;
end
end
end
With avoiding repeated calculation and a proper pre-allocation the code is at least efficient without parallelizing.
preallocation of all the arrays and structures have been done but not shown here.

请先登录,再进行评论。

回答(2 个)

Hello Jan,
I have preallocated those variables as well as the structures . I didn't show the complete code.

2 个评论

Please post comments inthe section for comments. Not showing the relevant part of the code, especially the pre-allocation, is not helpful for the readers.
Even without the double angle calculation, each inner loop takes around 0.5 secs on a 2.6 ghz core. i would have to run approx 2000 such inner loops and 2000 such outer loops

请先登录,再进行评论。

类别

帮助中心File Exchange 中查找有关 Historical Contests 的更多信息

产品

版本

R2015b

重新打开:

2018-12-22

Community Treasure Hunt

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

Start Hunting!

Translated by