Speeding up code for lattice generation

Question

Below is a lattice point generation code I write

Clear[lattice];
lattice[basisvec_List, numofcell_List, base_List] := 
  Flatten[Transpose[# + Transpose[base]] & /@ (Tuples[Range /@ numofcell].basisvec), 1]

to use it for example,

a = 1.;
list = lattice[{{3 a, 0}, {0, Sqrt[3] a}}, {5, 
     5}, {{0, 0}, {a, 0}, {-a/2, Sqrt[3] a/2}, {3 a/2, 
      Sqrt[3] a/2}}]; 

generate a honeycomb lattice

To generate 4000000 points,

lattice[{{3 a, 0}, {0, Sqrt[3] a}}, {1000, 
    4000}, {{0, 0}, {a, 0}, {-a/2, Sqrt[3] a/2}, {3 a/2, 
     Sqrt[3] a/2}}]; // AbsoluteTiming

it takes 4.52389 sec

This is the fastest code I could write at the moment. I also tried compile as follows

latticecom = 
  Compile[{{basisvec, _Real, 2}, {numofcell, _Integer, 
     1}, {base, _Real, 2}}, Clear[lattice];
   Flatten[
    Transpose[# + Transpose[base]] & /@ (Tuples[
        Range /@ numofcell].basisvec), 1], RuntimeOptions -> "Speed", 
   CompilationTarget -> "C"];

But it will gives errors when using it

CompiledFunction::cflist: Nontensor object generated; proceeding with uncompiled evaluation. >>

Can someone explain this?

I also write a fortran version to test how efficiency can differ between fortran and mathematica. It turns out that fortran done the same thing almost instantly.

So I want to see how this code could be great improved to be comparable to fortran.

Answer 1

This appears to be faster, dealing with the x and y coordinate lists separately:

lattice2[basisvec_List, numofcell_List, base_List] :=
 Module[{x, y, basex, basey},
  {x, y} = Transpose[Tuples[Range /@ numofcell].basisvec];
  {basex, basey} = Transpose[base];
  Transpose[{Join @@ ((x + #) & /@ basex), Join @@ ((y + #) & /@ basey)}]]

Answer 2

I keep improving my version

**First version **, by separate Transpose[base] and Change the way Transpose act, time reduce from 4.52389 sec to 2.57898 sec

Clear[lattice];
lattice[basisvec_List, numofcell_List, base_List] :=
 (tmptmp = Transpose[base];
  Flatten[
   Transpose[(# + tmptmp) & /@ (Tuples[
        Range /@ numofcell].basisvec), {1, 3, 2}], 1])

second version, reverse vector adding order, time reduce from 4.52389 sec to 1.28168 sec

Clear[lattice];
lattice[basisvec_List, numofcell_List, base_List] :=
 (basetmp = N@base;
  tmptmp2 = Transpose[Tuples[Range /@ numofcell].basisvec];
  Flatten[Transpose[(# + tmptmp2) & /@ basetmp, {1, 3, 2}], 1]

Notice that I add basetmp = N@base, because I found that if this is not added, there will be unpacking process. This is because there is exact zero in base. But don't know why the "first version" has no unpacking problem. What is more, if we change Flatten[Transpose[(# + tmptmp2) & /@ basetmp, {1, 3, 2}], 1] to Flatten[Transpose[(# + tmptmp2)] & /@ basetmp, 1], it will get unpacking again, don't know why!!!

And finally, still cannot catch fortran.