Want advice on writing better code to do matrix calculation

Question

I want to perform the computation shown in the code given below in a much faster way.

I have three loops for e, i and j. The matrices B1e, B2e, ...., B6e are functions of xi and eta. There are two ReplaceAll expressions which will replace the new value of xi and eta in the B1e, B2e, ...., B6e matrices, which are calculated repeatedly for the each i and j in the loops. But it is too slow.

Is there a way to make it faster? Maybe parallel computation? I'm open to advice.

\[Alpha] = MatrixForm[Table[i^2, {i, 5}]];

D1 = Table[i*j, {i, 3}, {j, 3}];

B1 = Table[\[Xi]^(i + e)*\[Eta]^(j + e), {i, 3}, {j, 3}, {e, 360}];
Ke1 = ConstantArray[0, {3, 3, 360}];

For[e = 1, e < 361, e++,
  For[i = 1, i < 6, i++,
    For[j = 1, j < 6, j++,
      Ke1[[All, All, e]] = 
        ReplaceAll[
         ReplaceAll[(Transpose[
             B1[[All, All, e]]].D1.B1[[All, All, 
              e]]), \[Xi] -> \[Alpha][[1, i]]], \[Eta] -> \[Alpha][[1,
             j]]];
      ];
    ];
  ];

---

@Henrik Schumacher the problem is that my ξ and η are functions of i,j and e but they are defined in matrices such as Alpha, U and V in which their dimensions are already known. Here Compile[{{i, _Integer}, {j, _Integer}, {e, _Integer}},if I dont mention the limits of i, j and e while compiling the function f, there will be a problem. Lets say na[[1,e]] matrix is defined for e=6 and for higher values of e>6 can not be found in na[[1,6]],na[[1,7]],..., and compilation can give an error. How could we fix it ? Below, you could see the code which is revised according to your solution strategy.

 Block[{i, j, e}, 
  f = {i, j, e} \[Function] 
    Evaluate[
     With[{ξ = ((U[[1, na[[1, e]] + 1]] - 
              U[[1, na[[1, e]]]])*\[Alpha][[NGaussxitilda[[1, 2]], 
             i]] + (U[[1, na[[1, e]] + 1]] + 
             U[[1, na[[1, e]]]]))*0.5, η = ((V[[1, 
               nb[[1, e]] + 1]] - V[[1, nb[[1, e]]]])*\[Alpha][[
             NGaussetatilda[[1, 2]], j]] + (V[[1, nb[[1, e]] + 1]] + 
             V[[1, nb[[1, e]]]]))*0.5}, (Transpose[
           B1e [[All, All, e]]].D1hat.B1e[[All, All, e]] + 
         Transpose[B1e[[All, All, e]]].D2hat.B2e[[All, All, e]] + 
         Transpose[B2e[[All, All, e]]].D3hat.B1e[[All, All, e]] + 
         Transpose[B2e[[All, All, e]]].D4hat.B2e[[All, All, e]] + 
         Transpose[B3e[[All, All, e]]].D5hat.B3e[[All, All, e]] + 
         Transpose[B3e[[All, All, e]]].D6hat.B4e[[All, All, e]] + 
         Transpose[B4e[[All, All, e]]].D7hat.B3e[[All, All, e]] + 
         Transpose[B4e[[All, All, e]]].D8hat.B4e[[All, All, e]] + 
         Transpose[B5e[[All, All, e]]].D9hat.B5e[[All, All, e]] + 
         Transpose[B6e[[All, All, e]]].D10hat.B6e[[All, All, e]])*
       detJe[[1, e]]
      ]]];


cf = With[{code = f[i, j, e]}, 
   Compile[{{i, _Integer}, {j, _Integer}, {e, _Integer}}, code, 
    CompilationTarget -> "C", RuntimeAttributes -> {Listable}, 
    Parallelization -> True, RuntimeOptions -> "Speed"]];


For[e = 1, e < nel + 1, e++,

  For[i = 1, i < NGaussxitilda[[1, 2]] + 1, i++,

    For[j = 1, j < NGaussetatilda[[1, 2]] + 1, j++,

      For[m = 1, m < nen*ndf + 1, m++,
        For[n = 1, n < nen*ndf + 1, n++,
          Kel[[m, n, e]] = Kel[[m, n, e]] + f[i, j, e];
          ];
        ];

      ];

    ];

  ];

Answer 1

A major issue with your code is that you do symbolic calculations over and over again within a threefold nested loop. Symbolic computations are slows. In the and you want to have a numerical array, so you get tid of symbolic computations as early as possible. Compute the entries of the resulting array once outside any loops and create a function depending on the loop indices. This could, e.g., like this:

Block[{i, j, e},

  f = {i, j, e} \[Function] 
    Evaluate[With[{ξ = i^2, η = j^2},

      (ξ η)^(2 e) Dot[Range[3],ξ^Range[3]]^2 TensorProduct[η^Range[3],η^Range[3]]]

     ]];

So, f[i,j,e] returns what the monstrosity

Ke1[[All, All, e]] = 
        ReplaceAll[
         ReplaceAll[(Transpose[
             B1[[All, All, e]]].D1.B1[[All, All, 
              e]]), ξ -> α[[1, i]]], η -> α[[1,
             j]]]

is supposed to compute. This can even be sped up by compiling the function cf:

cf = With[{code = f[i, j, e]},
   Compile[{{i, _Integer}, {j, _Integer}, {e, _Integer}},
    code,
    CompilationTarget -> "C",
    RuntimeAttributes -> {Listable},
    Parallelization -> True,
    RuntimeOptions -> "Speed"
    ]
   ];

The full array can be obtained, e.g., with

Table[cf[i,j,e],{i,1,3},{j,1,3},{e,1,360}]

I'll stop here since it is not exactly clear to me, what you are going to do with it. Be sure to get rid of the For loops: They are slow and do not scope their rinnung indices, so it is easy to write buggy code with it. Use Do or Table instead.

Addendum

Better not compile the code or at least not with RuntimeOptions -> "Speed". I have the very bad habit to use this option which deactivates---amongst other things---checks for overflow. And integer overflow is precisely what is happening here: The numbers that occur in the result are simply too large to be stored in machine integer type.