Avoiding 3 nested for loops?

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

[Jan]

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]

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?

Answer 1

[Bob Thompson]

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.

[Riccardo Sorvillo]

Hi Bob Nbob,

Firstly, thank you for your reply.

Actually, I couldn't understand your suggestions, even though I must thank you 'cause I discovered a new interesting MATLAB built-in function as repmat, that would be useful someday in the future.

Problem is that, with reference to my code above, only some positions of M2 must change iteratively according to mtmd(j), while with your code it returns error due to line

   M2(idb(35,2),idb(35,2)) = M2(idb(35,2),idb(35,2))+mtmd2;

Moreover, when it comes to C2 changes due to ctmd, if I would use no loop, or even a single loop, I would lose out-of-diagonal values, i.e. ctmd(i,j) with . Am I missing something from your approach?

Moving from your line, I thought that a good idea could have been to use 3D vectors and splitting for cycles in multiple points, as follows.

M2 = zeros(105,105,length(safe)); K2 = zeros(105,105,length(safe));
C2 = zeros(105,105,2*length(safe)); KFF2 = zeros(ndof,ndof,length(safe));
MFF2 = zeros(ndof,ndof,length(safe)); CFF2 = zeros(ndof,ndof,2*length(safe));
K_tmd = zeros(2,2,length(safe)); C_tmd = zeros(2,2,2*length(safe));
for j = 1:length(safe)
    [M2(:,:,j),K2(:,:,j)] = assem(incidenze2,l2,m2,EA2,EJ2,gamma2,idb2);
    K_tmd(:,:,j) = [ktmd(j) -ktmd(j);-ktmd(j) ktmd(j)];
    M2(idb(35,2),idb(35,2),(j)) = M2(idb(35,2),idb(35,2),(j))+mtmd(j); % contribution of mtmd
end
    K2([idb(6,2) idb(35,2)],[idb(6,2) idb(35,2)],:) = K2([idb(6,2) idb(35,2)],[idb(6,2) id2(35,2)],:)+K_tmd; % contribution of ktmd
    MFF2 = M2(1:ndof,1:ndof,:); KFF2 = K2(1:ndof,1:ndof,:);
for j = 1:2*length(safe)
    if j<=length(safe)
C2(:,:,j) = alfa*M2(:,:,j)+beta*K2(:,:,j);
    else
C2(:,:,j) = alfa*M2(:,:,j-length(safe))+beta*K2(:,:,j-length(safe));
    end
end

I am actually losing some information on the way and yet I am missing the construnction of damping matrix. Till this point, however, I can now run the code and waiting less than 1 minute.

Anyway, is there any smarter way to do it? Otherwise, how can I continue from this point?

[Bob Thompson]

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]

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.