6
$\begingroup$

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

enter image description here

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.

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

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)}]]
Improve this answer
$\endgroup$
1
  • $\begingroup$ my test shows that yours and mine are almost the same efficiency. Thank you! $\endgroup$ – matheorem Nov 7 '15 at 11:56
1
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.