Help with this triple nested for loop

Hi there, I have the following five matrices:

A1 = ones(5,5)*1;
A2 = ones(5,5)*2;
A3 = ones(5,5)*3;
A4 = ones(5,5)*4;
A5 = ones(5,5)*5;

I want a total of 96 seperate matrices and want to repeat A1,A2,A3,A4 and A5 after every 32 loops, Hence, I will have at the end 96 of 5 x 5 matrices consiting of A1, A2, A3, A4 and A5. So, the first 6 matrices (1-6) are of A1, the second lot of (7-12) matrices are A2, the third lot (13-18) are A3, the fouth (19 - 24) are A4 and A5 (25 - 32). After the 32nd matrix I want to go back to A1 and repeat as before till we get 96 matrices.

So, I have attempted the code below:

for i = 1:96
    for j = 1:3
        for k = 1:32
            U(:,:,k,j) = A1;
            Ux(:,:,i) = U(:,:,k,j);
            if k > 6
                U(:,:,k,j) = A2;
                Ux(:,:,i) = U(:,:,k,j);          
            end
            if k > 12
                U(:,:,k,j) = A3;
                Ux(:,:,i) = U(:,:,k,j);    
            end
            if k > 18
                U(:,:,k,j) = A4;
                Ux(:,:,i) = U(:,:,k,j);
            end
            if k > 24
                U(:,:,k,j) = A5;
                Ux(:,:,i) = U(:,:,k,j);
            end
        end
    end
end

Could someone help me acheive this please?

Many thanks in advance, guys!

,,

Answer 1

[Matt J]

I don't know what you're trying to do with Ux, but U is just,

U = cat(3, ...
    repelem(cat(3,A1,A2,A3,A4),1,1,6), ...
    repelem(A5,1,1,8));

U = repmat(U,1,1,1,3);

The reason the intent of Ux is not clear in the current code is because Ux(:,:,i) is overwritten for every j and k. Its final value is always the value at k = 32, namely A5.

[Scott Banks]

Hi Matt, the Ux is "supposed" to be where I store each of the 96 matrices that consist of the A1,A2,A3,A4 and A5.

I copy and pasted your code and its only giving me 32 matrices. I need 96, but that repeat in the pattern and order of the 32. So, like 3 lots of 32.

[Scott Banks]

I'm getting my desired out come doing this:

A1 = ones(5,5)*1;
A2 = ones(5,5)*2;
A3 = ones(5,5)*3;
A4 = ones(5,5)*4;
A5 = ones(5,5)*5;


    for j = 1:3
        for k = 1:32
            U(:,:,k,j) = A1;
             
            if k > 6
                U(:,:,k,j) = A2;
                
                
            end
                if k > 12
                    U(:,:,k,j) = A3;
                   
                     
                end
                    if k > 18
                        U(:,:,k,j) = A4;
                        
                        
                    end
                        if k > 24
                            U(:,:,k,j) = A5;
          end
        end
    end

[Torsten]

A1 = ones(5)*1;
A2 = ones(5)*2;
A3 = ones(5)*3;
A4 = ones(5)*4;
A5 = ones(5)*5;
U = zeros(5,5,32,3);
for j = 1:3
    for k = 1:6
        U(:,:,k,j) = A1;
    end
    for k = 7:12
        U(:,:,k,j) = A2;
    end
    for k = 13:18
        U(:,:,k,j) = A3;
    end
    for k = 19:24
        U(:,:,k,j) = A4;
    end
    for k = 25:32
        U(:,:,k,j) = A4;
    end
end

[Matt J]

I'm getting my desired out come doing this:

That gives identical results to what I posted:

A1 = ones(5,5)*1;
A2 = ones(5,5)*2;
A3 = ones(5,5)*3;
A4 = ones(5,5)*4;
A5 = ones(5,5)*5;

isequal(version1(A1,A2,A3,A4,A5), ...
        version2(A1,A2,A3,A4,A5))

ans = logical
   1

function U=version1(A1,A2,A3,A4,A5)

    for j = 1:3
        for k = 1:32
            U(:,:,k,j) = A1;
             
            if k > 6
                U(:,:,k,j) = A2;
                
                
            end
                if k > 12
                    U(:,:,k,j) = A3;
                   
                     
                end
                    if k > 18
                        U(:,:,k,j) = A4;
                        
                        
                    end
                        if k > 24
                            U(:,:,k,j) = A5;
          end
        end
    end


end


function U=version2(A1,A2,A3,A4,A5)

    U = cat(3, ...
        repelem(cat(3,A1,A2,A3,A4),1,1,6), ...
        repelem(A5,1,1,8));
    
    U = repmat(U,1,1,1,3);
 
end

[Scott Banks]

Hi guys,

this was part of a much bigger task, and I am hoping if it's okay to ask a follow up?

So I have this section of code from my script:

%% Transformation from local to global matrix

T = cat(3, Tcx, Tbx, Tdx1, Tdx2);

for i = 1:352
  for j = 5
    for k = 32
       Kmc(:,:,k,j) = T(:,:,i)'*Kc(:,:,1)*T(:,:,i);
     if k > 6
        Kmc(:,:,k,j) = T(:,:,i)'*Kc(:,:,2)*T(:,:,i);
     end
     if k > 12
        Kmc(:,:,k,j) = T(:,:,i)'*Kc(:,:,3)*T(:,:,i);
     end
     if k > 18
        Kmc(:,:,k,j) = T(:,:,i)'*Kc(:,:,4)*T(:,:,i);
     end
     if k > 24
        Kmc(:,:,k,j) = T(:,:,i)'*Kc(:,:,5)*T(:,:,i);
     end
    end
  end
     if i > 160
       Km(:,:,i) = T(:,:,i)'*Kb*T(:,:,i);
     end
     if i > 288
       Km(:,:,i) = T(:,:,i)'*Kd*T(:,:,i);
     end
       if i > 320
       Km(:,:,i) = T(:,:,i)'*Kd*T(:,:,i);
       end
end

I don't get any errors, and I end up with a 495 x 495 x 32 x 5 for Kmc which is what I want. However, the Kmc matrices are just filled with zeros. There are no values in them at all.

Hence, I am wondering why this is?

Many thanks, guys!

[Torsten]

Kmc(:,:,32,5) is changed here:

     if k > 24
        Kmc(:,:,k,j) = T(:,:,i)'*Kc(:,:,5)*T(:,:,i);
     end
     

but we don't know the right-hand side

     T(:,:,i)'*Kc(:,:,5)*T(:,:,i);

[Matt J]

this was part of a much bigger task, and I am hoping if it's okay to ask a follow up?

Yes, but does that mean the original question has been answered. If so, please Accept-click the answer.

[Scott Banks]

@Torsten

Hi Torsten, it quite a large amount of code, but all this should give you what the right-hand side should equal.

E = 199947961.502;
Ic = [29997.921/100^4; 22528.963/100^4; 17510.628/100^4; 14268.001/100^4; 11407.468/100^4];
Ib = 4.82E-04;
Ac = [0.0212858; 0.0168; 0.0136; 0.0113; 0.0093];
Ab = 0.0115528;
Ad = 0.0254469;
L = 6;
h = 3.5;
alpha = atand(3.5/6);
beta = atand(6/3.5);
d = sqrt(h^2 + L^2);
H = 112;
ns = 112/h;

n = (ns+1)*5*3

%% Local stiffness matrices

for i = 1:5

Kc(:,:,i) = [E*Ac(i)/h 0 0 -E*Ac(i)/h 0 0;
    0 12*E*Ic(i)/h^3 6*E*Ic(i)/h.^2 0 -12*E*Ic(i)/h^3 6*E*Ic(i)/h^2;
    0 6*E*Ic(i)/h^2 4*E*Ic(i)/h 0 -6*E*Ic(i)/h^2 2*E*Ic(i)/h;
    -E*Ac(i)/h 0 0 E*Ac(i)/h 0 0;
    0 -12*E*Ic(i)/h^3 -6*E*Ic(i)/h^2 0 12*E*Ic(i)/h^3 -6*E*Ic(i)/h^2;
    0 6*E*Ic(i)/h^2 2*E*Ic(i)/h 0 -6*E*Ic(i)/h^2 4*E*Ic(i)/h];

end

Kb = [E*Ab/L 0 0 -E*Ab/L 0 0;
    0 12*E*Ib/L^3 6*E*Ib/L^2 0 -12*E*Ib/L^3 6*E*Ib/L^2;
    0 6*E*Ib/L^2 4*E*Ib/L 0 -6*E*Ib/L^2 2*E*Ib/L;
    -E*Ab/L 0 0 E*Ab/L 0 0;
    0 -12*E*Ib/L^3 -6*E*Ib/L^2 0 12*E*Ib/L^3 -6*E*Ib/L^2;
    0 6*E*Ib/L^2 2*E*Ib/L 0 -6*E*Ib/L^2 4*E*Ib/L]

Kd = zeros(6,6);
Kd(1,1) = E*Ad/d;
Kd(1,4) = -E.*Ad/d;
Kd(4,1) = -E.*Ad/d;
Kd(4,4) = E*Ad/d

%% Transformation matrices

Tc1 = zeros(6,18);
Tc1(2,1) = -1;
Tc1(1,2) = 1;
Tc1(3,3) = 1;
Tc1(5,16) = -1;
Tc1(4,17) = 1;
Tc1(6,18) = 1;

Tb1 = zeros(6,6);
Tb1(1,1) = 1;
Tb1(2,2) = 1;
Tb1(3,3) = 1;
Tb1(4,4) = 1;
Tb1(5,5) = 1;
Tb1(6,6) = 1;

Td1 = zeros(6,21);
Td1(1,1) = cosd(alpha);
Td1(1,2) = sind(alpha);
Td1(4,19) = cosd(alpha);
Td1(4,20) = sind(alpha);

Td2 = zeros(6,15);
Td2(1,1) = cosd(beta + 90);
Td2(1,2) = sind(beta + 90);
Td2(4,13) = cosd(beta + 90);
Td2(4,14) = sind(beta + 90);

%% Counters for member nodes for columns

vecC = 1:5;
countC = 0:ns;
vecN1 = vecC + 5.*countC'
nM = size(vecN1,2) * (size(vecN1,1)-1);
nC = 3 * max(vecN1(:));
Tcx = zeros(size(Tc1,1),nC,nM);
j = 0

% For loop to complete the transformation matrices for columns

for i = 1:5
    R = vecN1(:,i);
    for k = 1:ns
        j = j + 1;
        C = (3*R(k) - 2):(3*R(k+1));
        Tcx(:,C,j) = Tc1;
    end
end

%% Counters for member nodes for beams

vecB = 6:10;
countB = 0:ns-1;
vecN2 = [vecB + 5.*countB']';
nM = 128;
nC = 3 * max(vecN2(:));
Tbx = zeros(size(Tb1,1),nC,nM);
j = 0

% For loop to complete the transformation matrices for beams

for i = 1:size(vecN2,2)
    R = vecN2(:,i);
    for k = 1:4
        j = j + 1;
        C = (3*R(k) - 2):(3*R(k+1));
        Tbx(:,C,j) = Tb1;
    end
end

%% counters for member nodes for diagonal braces

countD = 0:ns;
vecD1 = [2 + 5.*countD]'
vecD2 = [3 + 5.*countD]'
vecD3 = [4 + 5.*countD]'

for i = 1:length(vecD1)-1
  vecN3(:,i) = [vecD1(i);vecD2(i+1)];
  vecN4(:,i) = [vecD3(i);vecD2(i+1)];
  end

% For loop to complete the transformation matrices for diagonal braces

nM = 32;
nC = n;
Tdx1 = zeros(size(Td1,1),nC,nM);
Tdx2 = zeros(size(Td2,1),nC,nM);
j = 0


for i = 1:size(vecN3,2)
    R = vecN3(:,i);
    D = vecN4(:,i);
    for k = 1
        j = j + 1;
        C = (3*R(k) - 2):(3*R(k+1));
        B = (3*D(k) - 2):(3*D(k+1));
        Tdx1(:,C,j) = Td1;
        Tdx2(:,B,j) = Td2;
    end
end

%% Transformation from local to global matrix

T = cat(3, Tcx, Tbx, Tdx1, Tdx2);

for i = 1:352
  for j = 5
    for k = 32
       Kmc(:,:,k,j) = T(:,:,i)'*Kc(:,:,1)*T(:,:,i);
     if k > 6
        Kmc(:,:,k,j) = T(:,:,i)'*Kc(:,:,2)*T(:,:,i);
     end
     if k > 12
        Kmc(:,:,k,j) = T(:,:,i)'*Kc(:,:,3)*T(:,:,i);
     end
     if k > 18
        Kmc(:,:,k,j) = T(:,:,i)'*Kc(:,:,4)*T(:,:,i);
     end
     if k > 24
        Kmc(:,:,k,j) = T(:,:,i)'*Kc(:,:,5)*T(:,:,i);
     end
    end
  end
     if i > 160
       Km(:,:,i) = T(:,:,i)'*Kb*T(:,:,i);
     end
     if i > 288
       Km(:,:,i) = T(:,:,i)'*Kd*T(:,:,i);
     end
       if i > 320
       Km(:,:,i) = T(:,:,i)'*Kd*T(:,:,i);
       end
    end

Many thanks

[Scott Banks]

All done, Matt.

[Matt J]

However, the Kmc matrices are just filled with zeros. There are no values in them at all. Hence, I am wondering why this is?

Because k is not changing throughout the 'loop'. It is set to a fixed value of 32, which means that the conditions k>6,..., k>24 are never met.

In any case, you should not be using loop structures like this. Everything you are currently doing can be done looplessly with pagemtimes, cat(), and repmat().

[Torsten]

@Scott Banks

If you add the lines

M = Kmc(:,:,32,5)
M = M(:)
norm(M)

at the end of your code, you will see that Kmc is not all zero.

[Scott Banks]

@Matt J

Ah, that's a typo. J should be 1:5 and k should be 1:32.

I am still getting all zeros, and @Torsten, that gives me a column vector of 245025 x 1.

[Torsten]

I am still getting all zeros, and @Torsten, that gives me a column vector of 245025 x 1.

Yes, and the fact that the norm of M is not zero shows that not all elements of Kmc(:,:,32,5) can be zero.

M is just Kmc(:,:,32,5), written as a column vector:

495*495

ans = 245025

So you don't get all zeros.

[Matt J]

@Scott Banks In addition to what Torsten said, lines like the following make no sense,

 Kmc(:,:,k,j) = T(:,:,i)'*Kc(:,:,1)*T(:,:,i);

Since the right-hand side depends on i, but the left-hand side does not, the computations for all i=1...251 are discarded.

This might be what you intended,

T = cat(3, Tcx, Tbx, Tdx1, Tdx2);   % n x m x 352
nT = size(T, 3);                     % 352

%% Construct the Kc page assigned to each k
% k =  1:6   -> Kc(:,:,1)
% k =  7:12  -> Kc(:,:,2)
% k = 13:18  -> Kc(:,:,3)
% k = 19:24  -> Kc(:,:,4)
% k = 25:32  -> Kc(:,:,5)

kcPage = repelem(1:5, [6 6 6 6 8]);

Kck = Kc(:,:,kcPage);                % n x n x 32

%% Compute T_i' * Kc_k * T_i for every (k,i)

% Arrange T as n x m x 1 x 352
Tp = reshape(T, size(T,1), size(T,2), 1, nT);

% Arrange Kc pages as n x n x 32 x 1
Kcp = reshape(Kck, size(Kck,1), size(Kck,2), 32, 1);

% MATLAB implicitly expands the singleton page dimensions:
%   Kcp: n x n x 32 x 1
%   Tp : n x m x  1 x 352
KT = pagemtimes(Kcp, Tp);            % n x m x 32 x 352

% 'ctranspose' applies to each matrix page.
Kmc0 = pagemtimes( ...
    Tp, "ctranspose", ...
    KT, "none");                     % m x m x 32 x 352

% Replicate across the unused j dimension
Kmc = repmat(reshape(Kmc0, ...
    size(Kmc0,1), size(Kmc0,2), 32, 1, nT), ...
    1, 1, 1, 5, 1);                 % m x m x 32 x 5 x 352



%%%%%%%%%%%%%%%
Km = zeros(size(T,2), size(T,2), nT, "like", T);

ib = 161:288;
id = 289:352;

Km(:,:,ib) = pagemtimes( ...
    T(:,:,ib), "ctranspose", ...
    pagemtimes(Kb, T(:,:,ib)), "none");

Km(:,:,id) = pagemtimes( ...
    T(:,:,id), "ctranspose", ...
    pagemtimes(Kd, T(:,:,id)), "none");