# Want advice on writing better code to do matrix calculation [closed]

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];
];
];

];

];

];


## closed as off-topic by m_goldberg, Yves Klett, José Antonio Díaz Navas, J. M. will be back soon♦Mar 22 '18 at 0:38

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question cannot be answered without additional information. Questions on problems in code must describe the specific problem and include valid code to reproduce it. Any data used for programming examples should be embedded in the question or code to generate the (fake) data must be included." – m_goldberg, Yves Klett, José Antonio Díaz Navas, J. M. will be back soon
If this question can be reworded to fit the rules in the help center, please edit the question.

• Don't use For-loops. They are very slow. Replace item by item operations with functional list manipulation operation or array manipulation operation. – m_goldberg Mar 19 '18 at 6:50
• Could you give me an example from the given code ? How could I rearrange, just for one loop, m_goldberg ? – user45055 Mar 19 '18 at 6:53
• @user45055: Very difficult to give an example from your code because we cannot evaluate your code. Try giving us a smaller example with all parameters defined. – Lotus Mar 19 '18 at 7:19
• No. Your code is incomplete. It will not run as posted. There is no way I could rewrite it to work with functional constructs without having definitions for all the quantities referenced. – m_goldberg Mar 19 '18 at 7:19
• @user45055, how would you write the matrix operations with pen and paper? If there is a compact way to express it on paper, it should be possible to do it almost directly like that in Mathematica. – Kiro Mar 19 '18 at 7:57

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.

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.
• Something like Outer[cf,Range[3],Range[3],Range[360]] should compute all needed values at once. Of course, you are free to replace 360 with 100. – Henrik Schumacher Mar 19 '18 at 12:52