# Speeding up code for lattice generation

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 a}}, {5,
5}, {{0, 0}, {a, 0}, {-a/2, Sqrt a/2}, {3 a/2,
Sqrt a/2}}];


generate a honeycomb lattice To generate 4000000 points,

lattice[{{3 a, 0}, {0, Sqrt a}}, {1000,
4000}, {{0, 0}, {a, 0}, {-a/2, Sqrt a/2}, {3 a/2,
Sqrt 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.

• Tuples[] is not compilable, no? – J. M.'s torpor Nov 7 '15 at 8:50
• @J.M. Yeah I know, but Tuples is not in a loop, so I think it is harmless – matheorem Nov 7 '15 at 8:57
• Anyway, try this: lattice[basisvec_?MatrixQ, numofcell_?VectorQ, base_?MatrixQ] := Flatten[Through[(TranslationTransform /@ (Tuples[Range /@ numofcell].basisvec))[base]], 1]. – J. M.'s torpor Nov 7 '15 at 9:09
• @J.M. Nope, much slower than mine :） – matheorem Nov 7 '15 at 10:13
• This should be pretty fast, then: lattice[basisvec_List, numofcell_List, base_List] := Flatten[Outer[Plus, Tuples[Range /@ numofcell].basisvec, base, 1], 1]. – J. M.'s torpor Nov 7 '15 at 10:28

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)}]]

• my test shows that yours and mine are almost the same efficiency. Thank you! – matheorem Nov 7 '15 at 11:56

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.