Avoiding 3 nested for loops?

12 views (last 30 days)
Hello world!
I know this question has been already widely treated, but I am still struggling to grasp how to solve these types of problems.
Below the lines of my code that are super slow (mean time to finish of the order of 10 minutes). I guess that the inversion of matrix at each loop is slowing things down tremendously, but that the inner loop is nearly impossible to avoid (at least in my opinion); I see that the outer two can be vectorilized somehow, or at least worked out so to be simplified... but I don't know how to do it. The if loops should not be so influent, but I may be wrong.
ndof = 100; % number of degrees of freedom
vect_f = 0:.02:7; % frequency vector
omeg_n = 5.3; % natural freq.
safe = 1/5:.1:1; % safety factor
limamp = 1; % here's actually max(FRF(vect_f)), but for time being...
Mtot = 5000;
limM = 1.2*Mtot;
alpha = .15;
beta = 2*10-4;
amp2 = zeros(length(safe),length(safe)); % max FRF amplitude with TMD
redct = zeros(length(safe),length(safe)); % percentage of reduction with TMD
xx2 = zeros(ndof,length(vett_f));
F_A2 = zeros(ndof,1);
M2 = zeros(105,105);
K2 = zeros(105,105);
K_ktmd = zeros(2,2);
idb = randi(35*3,[35,3]); % here's actually a matrix of indeces that is calculated in previous lines of my code
mtmd = 100.*safe;
ktmd = mtmd.*(omeg_n)^2;
htmd = .3.*safe;
ctmd = htmd'*mtmd*2*(omeg_n);
Mtot2 = Mtot+mtmd;
F_A2(30,1) = 1; % input force
for j = 1:length(safe)
for i = 1:length(safe)
if Mtot2(j)>limM % constraint on the total mass
disp('Mass chosen is too high.');
return
else
[M2,K2] = assem(incidenze,l,m,EA,EJ,gamma,idb); % here's a function needed to create M2 and K2
M2(idb(35,2),idb(35,2)) = M2(idb(35,2),idb(35,2))+mtmd(j); % contribution of mtmd
K_ktmd = [ktmd(j) -ktmd(j);-ktmd(j) ktmd(j)];
K2([idb(6,2) idb(35,2)],[idb(6,2) idb(35,2)]) = ... % contribution of ktmd
K2([idb(6,2) idb(35,2)],[idb(6,2) idb(35,2)])+K_ktmd;
MFF2 = M2(1:ndof,1:ndof);
KFF2 = K2(1:ndof,1:ndof);
C2 = alfa*M2+beta*K2;
C_ctmd = [ctmd(j,i) -ctmd(j,i); -ctmd(j,i) ctmd(j,i)];
C2([idb(6,2) idb(35,2)],[idb(6,2) idb(35,2)]) = ... % contribution of ctmd
C2([idb(6,2) idb(35,2)],[idb(6,2) idb(35,2)])+C_ctmd;
CFF2 = C2(1:ndof,1:ndof);
for k = 1:length(vett_f)
A = -vett_f(k)^2*MFF2+sqrt(-1)*vett_f(k)*CFF2+KFF2;
xx2(:,k) = A\F_A2;
end
clear k
FRF_yA2 = xx2(idb(6,2),:);
amp2(j,i) = max(abs(FRF_yA2(1:ceil(length(FRF_yA2)/2))));
if amp2(j,i)<=limamp % this if has been added to trace the algorithm: returns multiple lines fastly, but there's a huge gap from 1 cycle to the other
redct(j,i) = (1-amp2(j,i)/amp)*100;
disp(['Amplitude vibration is reduced of ',num2str(redct(j,i)),...
'% with mmtd = ',num2str(.02.*safe(j)),'M_tot and hmtd = ',num2str(htmd(i))]);
else
redct(j,i) = (1-amp2(j,i)./amp)*100;
disp('Reduction of amplitude vibration is too low.')
end
end
end
end
clear j i
The function "assem.m" is a function built by my professor, which reads (don't know if you need it, so I upload just in case):
function [M,K] = assem(incidenze,l,m,EA,EJ,gamma,idb)
% Checking consistency input data
[n_el,nc]=size(incidenze);
if nc ~= 6
disp('Error: number of columns of incidence matrix different from 6')
return
end
if length(l) ~= n_el
sprintf('Error: number of elements in l differenc from n')
return
end
if length(m) ~= n_el
sprintf('Error: number of elements in m differenc from number of elements')
return
end
if length(EA) ~= n_el
sprintf('Error: number of elements in EA differenc number of elements')
return
end
if length(EJ) ~= n_el
sprintf('Error: number of elements in EJ differenc number of elements')
return
end
if length(gamma) ~= n_el
sprintf('Error: number of elements in alpha differenc number of elements')
return
end
if min(min(incidenze)) ~= 1
sprintf('Error: dof sequence does not start from 1')
return
end
n_doftot=max(max(idb));
% Assembling matrices M and K
M=zeros(n_doftot,n_doftot);
K=zeros(n_doftot,n_doftot);
for k=1:n_el
[mG,kG] = el_tra(l(k),m(k),EA(k),EJ(k),gamma(k));
for iri=1:6
for ico=1:6
i1=incidenze(k,iri);
i2=incidenze(k,ico);
M(i1,i2)=M(i1,i2)+mG(iri,ico);
K(i1,i2)=K(i1,i2)+kG(iri,ico);
end
end
end
Thank you in advance for your help, I appriacete it.
Riccardo
  2 Comments
Jan
Jan on 19 Jul 2019
Edited: Jan on 19 Jul 2019
Using sprintf to display an error message is fragile, because the code is not stopped. This is a very dangerous behavior of code, because proceeding to run after an error produces an undefined state. In addition creating the string does display it implicitly, but using a disp would be the direct way. Prefer a clean and strict error() for productive code.
Replace e.g.
disp('Mass chosen is too high.');
return
by
error('Mass chosen is too high.');
Riccardo Sorvillo
Riccardo Sorvillo on 21 Jul 2019
Thank you very much for the hint: avoiding dangerous behaviours is always more than welcome.
However, my main issue is the time needed to end the loop: any suggestion to reduce it?

Sign in to comment.

Accepted Answer

Bob Thompson
Bob Thompson on 18 Jul 2019
At first glance, I would agree that you can vectorize the two outer loops. After you create ctmd you can modify mtmd to be a repeat of the single row for the same number of columns. This for sure allows you to get rid of the first loop, because then you can call the entire row at once, and I think you can get rid of the second loop because you should be able to do all columns at once as well. I would suggest trying it.
mtmd = 100.*safe;
ktmd = mtmd.*(omeg_n)^2;
htmd = .3.*safe;
ctmd = htmd'*mtmd*2*(omeg_n);
mtmd2 = repmat(mtmd,[length(mtmd),1]); % Double check this creation, but should just duplicate the rows several times
Mtot2 = Mtot+mtmd;
% Remove first and second loops, and remove j and i index calls. Replace mtmd with mtmd2.
  3 Comments
Bob Thompson
Bob Thompson on 19 Jul 2019
I tend to use repmat most often when I am working with strings, reading text files, or other repetative data processing.
I'm not really sure what you mean by a 'smarter' way of doing things. Turning it into a 3D array was pretty good. I fully admit I am not an expert on all available options for manipulating matrices, but being able to eliminate one of the loops is already a big step.
I'm also not really sure what you mean by continuing. Have you not produced roughly the same results as before? How are things different that what you were expecting? I understand that you have 3D matrices now, but adjusting what you had before should simply be a matter of indexing.
Riccardo Sorvillo
Riccardo Sorvillo on 21 Jul 2019
I was able to improve only half of the code with those lines, still the FRF with the inversion of the matrix needed to be carry out, as well as the contribution of ctmd has to be inserted in C2 matrix.
While I was trying to figure it out, I launched my main code, waited for 20 mins or so and then, when it finished, stored the results in some ready-to-use .mat files. For time being, my advisor is fine with it, and so am I. Whenever I'll have the need to improve the code I will return on to it and post the final lines optimized.
In my opinion, starting from either 3D matrices or repmat and with the help of functions like ind2sub the final results can be achieved.
Thank you all for the help, apprieciated it.

Sign in to comment.

More Answers (0)

Products


Release

R2015b

Community Treasure Hunt

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

Start Hunting!