0
$\begingroup$

I have a function that behaves itself, but I want to make it much faster so I tried to put it into a compilable form After much fiddling I can't get it to work. Here's one of my failed attempts:

DiracPoint = Compile[{{angle, _Complex}},
   Module[{c1, c2, c3, Phi, cond1, cond2, cond3, sols, sol, theta, 
     phi},
    theta = Re[angle];
    phi = Im[angle];
    {c1, c2, c3} = 
     1 - 3 Sin[theta]^2 Cos[phi - 2 Pi (# - 1)/3]^2 & /@ {1, 2, 3};
    If[0 <= ((c2 + c3)^2 - c1^2)/(4 c2 c3) <= 1,
     Phi = Sqrt[((c2 + c3)^2 - c1^2)/(4 c2 c3)];
     cond1 = (c1 Cos[3 y/2] + (c2 + c3) Cos[Sqrt[3] x/2] == 0);
     cond2 = (-c1 Sin[3 y/2] + (-c2 + c3) Sin[Sqrt[3] x/2] == 0);
     cond3 = (Sin[Sqrt[3] x/2] == Phi);
     sols = Solve[cond1 && cond2 && cond3, {x, y}];
     {x, y} /. Refine[sols, {C[1] == 1, C[2] == 1}][[1]]
     , {$MaxMachineNumber, $MaxMachineNumber}
     ]
    ]
   ];

Why doesn't this work? I either get errors to do with incompatibility of rank sizes but I'm pretty confident the If statement is returning two objects of the same type and size. If I change {MaxM...,MaxM..} to just MaxM... I lose that error and get the error

"CompiledFunction::cfse: "Compiled expression y should be a !(\"machine-size real number\").""

If I then include x and y into my Module beginning declarations I get the error

Compile::initvar: The variable y has not been initialized or has been initialized to Null.

What is the proper way to initialise the variables x and y, or is my problem something else?

Here is how I apply it:

res = 4;
AngleArray = 
  Transpose[
   Table[theta + I phi, {theta, 0, Pi, Pi/res}, {phi, 0, 2 Pi, 
     Pi/res}]];
start = AbsoluteTime[];
array = Map[DiracPoint, AngleArray, {2}];
end = AbsoluteTime[] - start
$\endgroup$
  • $\begingroup$ @george2079 please don't feel obliged to respond but you were insightful on my previous, quite related post so thought you might have an idea on this one. Thanks. $\endgroup$ – Tom Mar 5 '15 at 0:33
  • 2
    $\begingroup$ relevant mathematica.stackexchange.com/questions/1096/… $\endgroup$ – george2079 Mar 5 '15 at 2:32
  • 1
    $\begingroup$ Compile isn't the catholicon for speed, the core part of your code is just not compilable. Generally Compile only speed up low level arithmetic calculation, see the link @george2079 gave for more details. $\endgroup$ – xzczd Mar 5 '15 at 3:26
  • $\begingroup$ Oh, so looks like Solve isn't even compatible. Bummer $\endgroup$ – Tom Mar 5 '15 at 10:32
  • 1
    $\begingroup$ As stated before you cannot compile Solve, but you might want to implement Newton's Method to solve your nonlinear system of equations. This iterative method can be compiled, no problem there. But then again it might be even faster or more accurate to use NSolve. $\endgroup$ – Wizard Mar 7 '15 at 17:57
2
$\begingroup$

In reply to @Wizard's comment. This is how I ended up solving/speeding up my code. I computed all the possible solutions to equation cond2 and cond3 and then picked whichever one first matched cond3. Probably an ugly way of doing it but seems to work.

DiracPoint[{theta_, phi_}] :=

  Module[{c1, c2, c3, Phi, x1, x2, y1a, y1b, y2a, y2b, Dx, Dy},
   {c1, c2, c3} = 
    1 - 3 Sin[theta]^2 Cos[phi - 2 Pi (# - 1)/3]^2 & /@ {1, 2, 3};
   {Dx, Dy} = If[0 <= ((c2 + c3)^2 - c1^2)/(4 c2 c3) <= 1,
     Phi = 
      Sqrt[((c2 + c3)^2 - c1^2)/(4 c2 c3)];(*Minus Phi gives K' point*)

          asx = ArcSin[Phi];
     asy1 = ArcSin[((c2 - c3)/c1) Sin[Sqrt[3] x1/2]];
     asy2 = ArcSin[((c2 - c3)/c1) Sin[Sqrt[3] x2/2]];
     x1 = (2/Sqrt[3]) asx;
     y1a = (2/3) asy1;
     y1b = (2/3) Sign[asy1] (Pi - Abs[asy1]);
     x2 = (2/Sqrt[3]) Sign[asx] (Pi - Abs[asx]);
     y2a = (2/3) asy2;
     y2b = (2/3) Sign[asy2] (Pi - Abs[asy2]);
     tol = 1*^-10;

     Which[
      Abs[(c2 + c3) Cos[Sqrt[3] x1/2] + c1 Cos[3 y1a/2]] < tol,
      {x1, y1a},
      Abs[(c2 + c3) Cos[Sqrt[3] x1/2] + c1 Cos[3 y1b/2]] < tol,
      {x1, y1b},
      Abs[(c2 + c3) Cos[Sqrt[3] x2/2] + c1 Cos[3 y2a/2]] < tol,
      {x2, y2a},
      Abs[(c2 + c3) Cos[Sqrt[3] x2/2] + c1 Cos[3 y2b/2]] < tol,
      {x2, y2b},
      1 == 1,
      {$MaxMachineNumber, $MaxMachineNumber}
      ]
     , {$MaxMachineNumber, $MaxMachineNumber}];
   {theta, phi, Dx, Dy}
   ];
$\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.