How can I make this program work faster?

Previous question: Can any one help me make my program work faster?

This question is an extension to the question I asked previously and referred to in the link above. I want to add spin to my system, so I changed it to the following code. the problem is that although I am using Compile, it is still very slow. I didn't change it much, so I don't understand why it is becoming so slow?

This is the code after including the spin by using the function kdfxn which acts as a delta function. gives 1 in case I have equal spin for instance if sz1 and sz2 are equal and gives 0 otherwise. The only changed I made is in chmat, the rest is the same as previous.

 nvec = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
19, 20, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 19, 20};
svec = {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, \
-1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1};
ne = 5;
nμ = 40;
δ = -150;
β = 1;
kdfxn[i_, j_] := If[i == j, 1, 0]
avec = {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0,
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0}};

cfxn = Block[{n1a, n1b, n2a, n2b},
With[{code = Which[
n1a == n1b && n2a == n2b,
Evaluate[N[1/6 (1 - 3/(n1a^2 π^2) - 3/(n2a^2 π^2))]],

n1a == n1b && n2a != n2b,
Evaluate[N[(4 (1 + (-1)^(n2a + n2b)) n2a n2b)/((n2a^2 - n2b^2)^2 π^2)]],

n1a != n1b && n2a == n2b,
Evaluate[N[(4 (1 + (-1)^(n1a + n1b)) n1a n1b)/((n1a^2 - n1b^2)^2 π^2)]],

True,
Evaluate[N[-((32 (-1 + (-1)^(n1a + n1b)) (-1 + (-1)^(n2a + n2b)) n1a n1b
n2a n2b)/((n1a^2 - n1b^2)^2 (n2a^2 - n2b^2)^2 π^4))]]
]},
Compile[{{n1a, _Integer}, {n1b, _Integer}, {n2a, _Integer}, {n2b,
_Integer}},
code,
CompilationTarget -> "C"
]
]];


This is the same function as previous question, which I did not change, and it will be use later to calculate a matrix:

   chmat = With[{ccfxn = cfxn, kkdfxn = kdfxn},
Compile[{{nm, _Integer}, {ne, _Integer}, {b, _Real}, {d, _Real}, \
{avec, _Real, 2}, {nvec, _Real, 1}, {svec, _Real, 1}},
Block[{sz0, sz1, sz2, sz3, n0, n1, n2, n3, h1, h2, tmp, tmp2,
tmp21, kf01, kf23, kf13, kf02},
Table[n0 = CompileGetElement[nvec, nm0];
n1 = CompileGetElement[nvec, nm1];
sz0 = CompileGetElement[svec, nm0];
sz1 = CompileGetElement[svec, nm1];
tmp = 0.;
Do[sz2 = CompileGetElement[svec, nm2];
sz3 = CompileGetElement[svec, nm3];
n2 = CompileGetElement[nvec, nm2];
n3 = CompileGetElement[nvec, nm3];
tmp2 = ccfxn[n1, n0, n3, n2];
tmp21 = ccfxn[n1, n3, n0, n2];

kf01 = kkdfxn[sz0, sz1];
kf23 = kkdfxn[sz2, sz3];
kf13 = kkdfxn[sz1, sz3];
kf02 = kkdfxn[sz0, sz2];

Do[
tmp += (tmp2 kf23 kf01 - tmp21 kf13 kf02) CompileGetElement[
avec, j, nm3] CompileGetElement[avec, j, nm2], {j, 1,
ne}], {nm2, 1, nm}, {nm3, 1, nm}];
d tmp +
If[nm0 == nm1, (n0^2 Pi^2 + b CompileGetElement[svec, nm0]),
0.], {nm0, 1, nm}, {nm1, 1, nm}]], CompilationTarget -> "C",
CompilationOptions -> {"InlineCompiledFunctions" -> True},
RuntimeOptions -> "Speed"]];


and then I use the iteration method to converge the energies, which is the same as previous question and the speed is fine,:

  Table[
hmat = chmat[nμ, ne, δ, β, avec, nvec, svec];
{evals, evecs} = Eigensystem[hmat];
pos = Ordering[evals][[1 ;; ne]];
bvec = Map[x \[Function] If[Total[x] < 0, -x, x], evecs[[pos]]];
residual = Max[Abs[avec - bvec]];
avec = bvec;
{residual, Total[evals[[pos]]]},
{j, 1, 30}]


I wish some one could tell me why does it become so slow as I have included spin, and is there any way I can make it work faster?

• @b3m2a1 Ok, Thank you. Sorry for the Junk!!! I am beginner, maybe I couldn't find a way to make it shorter. Anyway, Thanks for your help Commented Oct 23, 2017 at 21:05
• "junk" was too harsh a word. I just meant it's a lot. Generally we advocate breaking things down into the most minimal piece you can find. In this case that might be using a chmat adapted so that it used a smaller cfxn that had fewer parameters being fed into it. That would also have helped you find the issue, which was the kdfxn. Commented Oct 23, 2017 at 21:07
• @b3m2a1 sure, Thank you Commented Oct 23, 2017 at 21:09

Despite telling myself I wouldn't fix this, here you are... Your problem is that you can't Compile DownValues.
Replace kdfxn with (If[# == #2, 1, 0] &) and it'll work fine.
It might help to read up on Compile a bit so this won't be an issue in the future
P.S. KroneckerDelta already exists in the language and if you were to do your own version you'd want to use === instead of ==
• @DelaramNematollahi understandable. It's a very different way to think about the code. Just keep in mind that DownValues` can't be compiled and you'll be much better off. Commented Oct 23, 2017 at 21:08