How to speed up the computation of a specific array

Question

I am doing some calculations on superconductivity and I REALLY need to speed up the way I calculate one of my arrays. I will put the input to the matrix (of course a tiny toy model of what I am doing), just in case somebody wants to time it, but the important part, and the part which I am asking for help about is delta, which is in the lower part.

Input to the matrix:

nn = 2; Nband = 2;
Nstates = 2*nn*Nband;

eigenvectors = Table[Range[Nstates] + i, {i, 1, Nstates}];

InverseFlatten[l_, dimensions_] := Fold[Partition, Flatten@l, Most[Reverse[dimensions]]];
uv = InverseFlatten[eigenvectors, {Nstates, Nband, 2, nn}];
u = uv[[1 ;; Nstates, 1 ;; Nband, 1]];
v = uv[[1 ;; Nstates, 1 ;; Nband, 2]];
f = Range[Nstates];
V = Table[Which[i == j + 1, 2], 
    {l, 1, Nband}, {m, 1, Nband}, {s, 1, Nband},{q, 1, Nband}, {i, 1, nn}, {j, 1, nn}];

Now, first I naively started by doing in it in a intuitive way, but it turned to be crazily slow!!

delta = Table[
         Table[Which[i == j + 1, 
           Sum[V[[μ, ν, q, s, i, j]]*Sum[u[[n, q, i]]*v[[n, s, j]]*f[[n]],
           {n, 1, Nstates}],{q, 1, Nband}, {s, 1, Nband}]],
         {i, 1, nn}, {j, 1, nn}],
        {μ, 1, Nband}, {ν, 1, Nband}];

Then thanks to this forum, I learned that one should avoid Sum and that it was smarter to do Dot products of arrays and to use Total, then I wrote:

delta = 
 Table[
   Total[Flatten[V[[μ, ν, ;; , ;; , ;; , ;;]]*
      Table[
         Table[
            Which[
                i == j + 1,
                f.(u[[;; , q, i]]*Conjugate[v[[;; , s, j]]])    ], 
              {i, 1, nn}, {j, 1, nn}],
          {q, 1, Nband}, {s, 1, Nband}], 1]], 
     {μ, 1, Nband}, {ν, 1, Nband}];

But this is still toooo slow.

Does any of the bright minds of this forum see a faster way to compute delta? If somebody wants to time it using Timing, can make it bigger by increasing nn and Nband, in my real calculations nn goes up to 600 and Nband up to 5.

(I do not care about the Null elements, just need the non Null elements)

Thanks

Answer 1

Use ParallelTable command...

nn = 100; Nband = 5;
Nstates = 2*nn*Nband;

eigenvectors = ParallelTable[Range[Nstates] + i, {i, 1, Nstates}];

InverseFlatten[l_, dimensions_] := 
 Fold[Partition, Flatten@l, Most[Reverse[dimensions]]];
uv = InverseFlatten[eigenvectors, {Nstates, Nband, 2, nn}];
u = uv[[1 ;; Nstates, 1 ;; Nband, 1]];
v = uv[[1 ;; Nstates, 1 ;; Nband, 2]];
f = Range[Nstates];
V = ParallelTable[
Which[i == j + 1, 2], {l, 1, Nband}, {m, 1, Nband}, {s, 1, 
Nband}, {q, 1, Nband}, {i, 1, nn}, {j, 1, nn}];

delta1 = ParallelTable[
 Total[Flatten[
   V[[μ, ν, ;; , ;; , ;; , ;;]]*
    Table[Table[
      Which[i == j + 1, 
       f.(u[[;; , q, i]]*Conjugate[v[[;; , s, j]]])], {i, 1, 
       nn}, {j, 1, nn}], {q, 1, Nband}, {s, 1, Nband}], 
   1]], {μ, 1, Nband}, {ν, 1, Nband}] /. {Null -> 
  0}; // Timing

For me it gives result {11.407000, Null}

Without parallel command...

delta = Table[
 Total[Flatten[
   V[[μ, ν, ;; , ;; , ;; , ;;]]*
    Table[Table[
      Which[i == j + 1, 
       f.(u[[;; , q, i]]*Conjugate[v[[;; , s, j]]])], {i, 1, 
       nn}, {j, 1, nn}], {q, 1, Nband}, {s, 1, Nband}], 
   1]], {μ, 1, Nband}, {ν, 1, Nband}] /. {Null -> 
  0}; // Timing

{23.768000, Null}