Message Boards Message Boards


Parallelize code economically with Compile[]

Posted 10 years ago
0 Replies
6 Total Likes
The original idea of this article is from my primitive trial to solve Problem 143 of Project Euler. I am not presenting a solution to this problem, instead, I want to demonstrate how to easily launch multiple threads from compiled code without using extra licenses. 

Briefly speaking, I am looking for the triangle ABC with interger sides AB = c, BC= a and CA = b , provided that the distances between the first Fermat Point of this triangle and its vertices are integers (distances are denoted as p, q and r) and the sum p+q+r is preset with value s, as well. Finally I will sum up all distinct values of s under 120,000 which yields at least one integer tuple (a, b, c, p, q, r) . 

Here is my code: 
Because Compile[] only supports a limited number of Mathematica function, I need to define some internal function by myself. For instance, IntegerQ[Sqrt[<An Integer>]] will not work in the compiled code as the IntegerQ cannot be compiled. So this perfect square check function is rewrote in simple functions:
perfectSquareQ =
Compile[{{n, _Integer}},
  Module[{sq = Floor[Sqrt[Abs[n]]]}, If[sq^2 == n, True, False]]]
I can quickly check this function as Compile returns an functional object which is very similar to LinearModelFit 
In[1]:= perfectSquareQ[5]
Out[1]:= False
Then comes the main body of the procedure. I use "cf" to denote my compiled function. This "cf" basically checks if a given integer sum p+q+r can produce integer sides of a trianle. There may be multiple results or no result. "cf" terminates either when it finds the first solution or after it scans every possibilities in vain (a costly situation). 
 cf = Compile[{{limit, _Integer}},
   Module[{limP = Ceiling[limit/3], r = 1, a = 0, b = 0, c = 0,
     res = 0},
     If[res > 0, Break[]];
      r = limit - p - q;
      a = Abs[(limit - p)^2 - r*q];
      b = Abs[(limit - r)^2 - p*q];
     c = Abs[(limit - q)^2 - p*r];
     If[perfectSquareQ[a] && perfectSquareQ[b] && perfectSquareQ[c],
      res = limit; Break[]]
     , {q, Ceiling[(limit - p)/2], limit - p - 1}]
    , {p, limit - 2, limP, -1}];
  CompilationTarget -> "C",
  CompilationOptions -> {"InlineExternalDefinitions" -> True},
  RuntimeAttributes -> {Listable},
  Parallelization -> True
Besides the function code itself, you can see I have several options in the definition of "cf". 
  1. CompilationTarget: This option tells Mathematica to use the compiled C code instead of runing with "WVM" (wolfram virtual machine)
  2. CompilationOptions: I use this option because the "cf" calls another compile function "perfectSquareQ"
  3. RuntimeAttributes/Parallelization: These two options are critical to our parallel computation process. The symbol RuntimeAttributes is similar to SetAttributes, here it let the head "cf" take elements of a list automatically, ie. cf[{1,2,3}]  -> {cf[1],cf[2],cf[3]}. While Parallelization let cf[1] , cf[2] and cf[3] get executed on different cpu cores automatically. With these two options, cf[Range[1,100]] can be parallelized with no effort. 
As you run this code on your machine (with proper C compiler), you can see check the following statistics from the system monitor:

Of course this method still takes a rather long time to find the possibilites of p+q+r < 30,000 if you do not adopt a smart algorithm with number theory thinges. However, you may save a lot of time to do similar procedure programmings with many jumpy cases, which are not suitable for Map or Table function in Mathematica. 
POSTED BY: Shenghui Yang
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract