Problem explanation
I work with a symmetric matrix $M$ that consists of four single matrices. I calculate three of them such that $M$ = $\left(\begin{matrix}a & b \\ b^T & c\end{matrix}\right)$ and transpose $b$. The calculation of each matrix element is constant(!) in time, the Matrix scales according to $M\propto N^6$ where $N$ denotes an input parameter to the function that generates the matrix. The dimensions of the matrices are $a$: $N\times N$, $b$: $N \times N^3$, $c:$ $N^3 \times N^3$.
Expected as well as received scaling behaviour
My overall expectation is a scaling in time according to $N^6$ as this is the number of matrix elements that have to be calculated. I measured (expected):
- $N=13:$ 18.06 sec
- $N=14:$ 29.13 sec (28.17 sec)
- $N=20:$ 936.93 sec (239.46 sec)
As you can see for increasing $N$ the time noticeably deviates from the expectation. For small values of $N$ the approximation works quite well as can be seen with
exp[baseTime_, i_, n_] := baseTime*(i/6)^n;
baseTime = AbsoluteTiming[BuildMatrix[6]][[1]];
Table[{exp[baseTime, i, 2], exp[baseTime, i, 6], AbsoluteTiming[BuildMatrix[i]][[1]]}, {i, 7, 12}]
(* Out: {{0.590009, 1.09306, 1.09899}, {0.770624, 2.43555, 2.45416}, {0.975321, 4.93756, 4.989}, {1.2041, 9.2909, 9.27697}, {1.45696, 16.4594,16.3822}, {1.7339, 27.7425, 27.5057}}*)
and a quite perfect scaling behaviour according to a $\propto N^6$ expectation (compared to a wrongly assumed $\propto N^2$ behaviour here).
Code being used
The code I use is the following
t=1.;U=4.;
BuildMatrix[LatticeSize_] :=
Module[{m, n, o, i, j, k, matrix, aMatrix, bMatrix, bMatrixT, cMatrix},
aMatrix = ArrayReshape[Table[aElement[m, i, LatticeSize], {m, 0, LatticeSize - 1}, {i, 0, LatticeSize - 1}], {LatticeSize, LatticeSize}];
bMatrix = ArrayReshape[Table[bElement[i, j, k, p], {i, 0, LatticeSize - 1}, {j, 0, LatticeSize - 1}, {k, 0, LatticeSize - 1}, {p, 0, LatticeSize - 1}], {LatticeSize^3, LatticeSize}];
bMatrixT = Transpose[bMatrix];
cMatrix = ArrayReshape[Table[cElement[m, n, o, i, j, k, LatticeSize], {m, 0, LatticeSize - 1}, {n, 0, LatticeSize - 1}, {o, 0, LatticeSize - 1}, {i, 0, LatticeSize - 1}, {j, 0, LatticeSize - 1}, {k, 0, LatticeSize - 1}], {LatticeSize^3, LatticeSize^3}];
matrix = ArrayFlatten[{{aMatrix, bMatrixT}, {bMatrix, cMatrix}}];
matrix
];
where the methods defining the calculations are simple If-constructs of the form
aElement[m_, i_, LatticeSize_] := If[m == Mod[i - 1, LatticeSize], -t, 0.] + If[m == Mod[i + 1, LatticeSize], -t, 0.] + If[m == i, U/2, 0.]
bElement[i_, j_, k_, p_] := If[i == j == k == p, U/2, 0.];
cElement[m_, n_, o_, i_, j_, k_, LatticeSize_] := If[m == i && o == k && n == Mod[j + 1, LatticeSize] || n == j && o == k && m == Mod[i + 1, LatticeSize], -t, 0.] + If[m == i && o == k && n == Mod[j - 1, LatticeSize] || n == j && o == k && m == Mod[i - 1, LatticeSize], -t, 0.] + If[m == i && o == Mod[k + 1, LatticeSize] && n == j, t, 0.] + If[m == i && o == Mod[k - 1, LatticeSize] && n == j, t, 0.] + If[m == i && o == k && n == j, U/2, 0.];
where LatticeSize equals the formerly described input parameter $N$. Where is the bottleneck in my code? I am not able to see where I lose time in computation.
BuildMatrix[2]. Hard to give any ideas without working code. Very hard. $\endgroup$bElementandcElement. $\endgroup$