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

Dear Mencia, after 25 question you sure can do better with regard to formatting your post? — Yves Klett, Dec 13 '13 at 11:14
@Yves Klett I tried to figure that out by putting in this forum, "how to format a quiestion", but did not find anything. Where can I find that? I definitely want to format, forsurwe — Mencia, Dec 13 '13 at 11:16
Please have a look at: http://mathematica.stackexchange.com/help/formatting. And/Or have a look at the edits to your questions (or any other Q) which will show you how the markdown works. — Yves Klett, Dec 13 '13 at 11:18
Most importantly, for code you simply indent by four spaces (look at my minimal edit). — Yves Klett, Dec 13 '13 at 11:25
Now, for that pro feeling you could even replace the fancy typesetting garble like \[Mu] using @halirutan´s answer here: http://meta.mathematica.stackexchange.com/q/1043/131 — Yves Klett, Dec 13 '13 at 11:27
At a quick glance, I think the delta calculation could be compiled (http://reference.wolfram.com/mathematica/ref/Compile.html). The naive version is actually a better starting point. With CompilationTarget->"C" you should obtain good performance. — Ajasja, Dec 13 '13 at 11:42
@Ajasja thanks so much, I see if I can conbine "Compile" with the naive way to calculate delta. — Mencia, Dec 13 '13 at 11:50
@Mencia are you sure that both of your versions give the correct result? The results do not look identical. — C. E., Dec 13 '13 at 11:52
@Anon thanks for the remark. Precisely they are not identical in the Null elements, but as I said I will only use the non Null elements, therefore, those are the ones that should be identical, both the number and the position in the final array. — Mencia, Dec 13 '13 at 11:54
@Mencia It's always good to work with PackedArrays http://mathematica.stackexchange.com/questions/3496/what-is-a-mathematica-packed-array, which means that Nulls can't be part of the array. — Ajasja, Dec 13 '13 at 12:03
@Mencia Of course they do. Do read about PackedArray. The bottom line is that you have to have an array (a list) where all elements are of the same type. — Ajasja, Dec 13 '13 at 12:46
As an example ``(Array of one type) na = ConstantArray[Null, 100]; Print@ByteCount@na; Print@Developer`PackedArrayQ@na;
(Array of floats) na = ConstantArray[1., 100]; Print@ByteCount@na; Print@Developer`PackedArrayQ@na;

(make mixed array) na[[1]] = Null; Print@ByteCount@na; Print@Developer`PackedArrayQ@na;`` — Ajasja, Dec 13 '13 at 13:02
@Ajasja how do I create a packed array, I tried for example: V = Table[ Which[i == j + 1, 2, True, 0], {l, 1, Nband}, {m, 1, Nband}, {s, 1, Nband}, {q, 1, Nband}, {i, 1, nn}, {j, 1, nn}];Developer`ToPackedArray@V; but stil not packed. — Mencia, Dec 14 '13 at 18:20

score 1 · Accepted Answer · edited Sep 21 '16 at 10:45

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}

How to speed up the computation of a specific array

1 Answers1