How to speed up my code using Sparse-array possibly?

Question

I created two loops using Parallel table and a table. The outer loop causes a (nstates by nstates) dimension matrix and the inner loop calculates the matrix elements based on their position in the outer matrix and also the initial predefined array (avec). Keeping in mind the avec I wrote here is just an example and always not so simple, and I do not usually end up with a final diagonal matrix, I want to know how I can make use of SparsArray and probably Associate to speed it up? I have tried to use SparseArray, but since I need to calculate the elements initially based on their position and I am not that familiar with the tricks to do that, I would appreciate it if someone can help.

here is the example code:

initial conditions:

ℓ0 = 8;
γ = 
Join[Table[{m, 1}, {m, -ℓ0, ℓ0}], 
Table[{m, -1}, {m, -ℓ0, ℓ0}]]
nstates = Length[γ]
ne = 2 ℓ0 + 2

simple avec:

 avec = Table[0, {ie, 1, ne}, {i, 1, nstates}];
 Do[avec[[1, 3]] = 1;
 avec[[i + 1, 2 + ne]] = 1, {i, 1, (ne - 1)}]

This is the function I am using in the loop, and probably its at its best speed, so I do not think I need to change anything about this part:

 ParallelEvaluate[Off[ClebschGordan::phy];
 ClearAll[j3s];
 j3s[a_, b_, c_] := j3s[a, b, c] = ThreeJSymbol[a, b, c];
 ClearAll[dfxn]; 
 dfxn[ℓ_, m1_, m2_, p1_, p2_] := 
 N@If[m1 + p1 == m2 + p2, 
 Sum[(2 ℓ + 1)^2 (2 ℓtemp + 1)/(4 π )
    Sum[If[m1 + p1 == mval && m2 + p2 == mval, 
     j3s[{ℓ, m1}, {ℓ, 
        p1}, {ℓtemp, -mval}] j3s[{ℓ, 
        m2}, {ℓ, 
        p2}, {ℓtemp, -mval}] j3s[{ℓ, 
        0}, {ℓ, 0}, {ℓtemp, 0}]^2, 
     0], {mval, -ℓtemp, ℓtemp}], \
 {ℓtemp, 0, 2 ℓ}], 0];];

and here is the loop which I want to modify the speed using SparseArray if possible:

   vex =(*(2 ℓ0 +1)^2*) ParallelTable[

  mpf = γ[[f, 1]];
  mk = γ[[k, 1]];
  μpf = γ[[f, 2]];
  μk = γ[[k, 2]];
  Chop[Total[
  Table[(* Here we loop over the HF states 
  *)Off[ClebschGordan::phy];
  pi = γ[[i, 1]]; 
  pj = γ[[j, 1]];
  μpi = γ[[i, 2]];
  μpj = γ[[j, 2]];
  If[μpi == μk && μpj == μpf, 
  N[Conjugate[avec[[ie, i]]] *avec[[ie, j]]* 
    dfxn[ ℓ0, pi, mk, mpf, pj]], 0]
 , {ie, 1, ne}, {i, 1, nstates}, {j, 1, nstates}]
 , Infinity]],
 {f, 1, nstates}, {k, 1, nstates}]
------------------------------------------------------------------

This was the simple version of my program, Thanks to you @Henrik Schumacher I got to learn about packed arrays. But I am still having problem with making it work for different avecs in which I may have more than one nonzero element in each row. for example if I define my avec as

 avec = SparseArray@ConstantArray[0., {ne, nstates}];

 stateList = Flatten[Table[stateA1 = {im - 1, -1}; stateA2 = {-im, 1};
 stateB1 = {-im, -1}; stateB2 = {im, 1};
 iA1 = Part[Position[\[Gamma], stateA1], 1, 1];
 iA2 = Part[Position[\[Gamma], stateA2], 1, 1];
 iB1 = Part[Position[\[Gamma], stateB1], 1, 1];
 iB2 = Part[Position[\[Gamma], stateB2], 1, 1];
 {{iA1, iA2}, {iB1, iB2}}, {im, 1, \[ScriptL]0}], 1];
 Do[avec[[ie, stateList[[ie, 1]]]] = 
 Sin[(ie \[Pi])/(2 (2 \[ScriptL]0 + 1))];
 avec[[ie, stateList[[ie, 2]]]] = 
 Cos[(ie \[Pi])/(2 (2 \[ScriptL]0 + 1))];
, {ie, 1, 2 \[ScriptL]0}]
 avec[[2 \[ScriptL]0 + 2, Part[Position[\[Gamma], {0, 1}], 1, 1]]] = 1;
 avec[[2 \[ScriptL]0 + 1, 
 Part[Position[\[Gamma], {\[ScriptL]0, -1}], 1, 1]]] = 1;

Then I don't know how to change this part of your code

 aa = ConjugateTranspose[avec].avec;
 {ilist, jlist} = Transpose[aa["NonzeroPositions"]];

Or if anyway I can use the same method still. Thank you again for your time and help.

Answer 1

Okay, here is my take on it, but I give no guarantees for correctness. The basic idea here is to replace summations by logic, to do as much logic as possible offline and to do remaining by filtering lists of summation indices.

ℓ0 = 8;
γ = Developer`ToPackedArray@ Join[Table[{m, 1}, {m, -ℓ0, ℓ0}], Table[{m, -1}, {m, -ℓ0, ℓ0}]];
{γ1, γ2} = Transpose[γ];
nstates = Length[γ];
ne = 2 ℓ0 + 2;

avec = SparseArray@ConstantArray[0., {ne, nstates}];
Do[avec[[1, 3]] = 1.; avec[[i + 1, 2 + ne]] = 1., {i, 1, (ne - 1)}];

aa = ConjugateTranspose[avec].avec;
{ilist, jlist} = Transpose[aa["NonzeroPositions"]];
vals = aa["NonzeroValues"];
γ1ilist = γ1[[ilist]];
γ1jlist = γ1[[jlist]];
γ2ilist = γ2[[ilist]];
γ2jlist = γ2[[jlist]];

ParallelEvaluate[Off[ClebschGordan::phy];
  ClearAll[j3s];
  j3s[a_, b_, c_] := j3s[a, b, c] = N@ThreeJSymbol[a, b, c];
  ClearAll[dfxn];
  dfxn[ℓ_, m1_, m2_, p1_, p2_] := 
   If[m1 + p1 == m2 + p2, 
    Sum[(2 ℓ + 1)^2 (2 ℓtemp + 1)/(4 π) Sum[
       If[m1 + p1 == mval && m2 + p2 == mval, 
        j3s[{ℓ, m1}, {ℓ, p1}, {ℓtemp, -mval}] j3s[{ℓ, m2}, {ℓ, p2}, {ℓtemp, -mval}] j3s[{ℓ, 0}, {ℓ, 0}, {ℓtemp, 0}]^2, 
        0.], {mval, -ℓtemp, ℓtemp}], 
{ℓtemp, 0, 2 ℓ}], 0.];];

And now

vex2 = SparseArray@ParallelTable[
    γ1f = γ1[[f]];
    γ1k = γ1[[k]];
    γ2f = γ2[[f]];
    γ2k = γ2[[k]];
    idx = Intersection[
      Random`Private`PositionsOf[γ2ilist, γ2k],
      Random`Private`PositionsOf[γ2jlist, γ2f]
      ];
    vals[[idx]].MapThread[
      dfxn[ℓ0, #1, γ1k, γ1f, #2] &,
      {γ1ilist[[idx]], γ1jlist[[idx]]}
      ],
    {f, 1, nstates}, {k, 1, nstates}];

This gets the job done in about 0.15 s (OP's code required about 100 s) on my machine. For the data provided by OP, vex and vex2 coincide. I am pretty sure that dfxn can be optimized further; it still contains summation of many zeroes. But I leave it as it is for the moment.