title: Parallelism on divide & conquer algorithm
body: <p>I'm facing problems to make my code runs in parallel. It is a 3D Delaunay generator using a divide &amp; conquer algorithm named <a href="http://www.google.be/url?sa=t&amp;source=web&amp;cd=5&amp;ved=0CD8QFjAE&amp;url=http%3A%2F%2Fce.sharif.edu%2F~ghodsi%2Farchive%2FCourse-Materials%2FComputational%2520Geometry%2Fpapers%2FDeWall-A%2520fast%2520Divide%2520and%2520Conquer%2520DT%2520Alg%2520in%2520Ed.ps&amp;ei=bPZOTaTSBMyFhQf3npjGDg&amp;usg=AFQjCNFdJGybWug6CQNDmOddg116iJP80Q&amp;sig2=owvTJyG_lLqNKXt3SleOLQ" rel="nofollow">DeWall</a>.</p>
  
  <p>The main function is:</p>
  
  <pre><code>deWall::[SimplexPointer] -&gt; SetSimplexFace -&gt; Box -&gt; StateT DeWallSets IO ([Simplex], [Edge])
  deWall p afl box = do
     ...
     ...
     get &gt;&gt;= recursion box1 box2 p1 p2 sigma edges
     ...
     ...
  </code></pre>
  
  <p>It calls the "recursion" function that might call the dewall function back. And it is here where the parallization opportunity appears. The following code shows the sequential solution.  </p>
  
  <pre><code>recursion::Box -&gt; Box -&gt; [SimplexPointer] -&gt; [SimplexPointer] -&gt; [Simplex] -&gt; [Edge] -&gt; DeWallSets -&gt; StateT DeWallSets IO ([Simplex], [Edge])    
  recursion box1 box2 p1 p2 sigma edges deWallSet
          | null afl1 &amp;&amp; null afl2 = return (sigma, edges)
          | (null) afl1 = do
              (s, e) &lt;- deWall p2 afl2 box2
              return (s ++ sigma, e ++ edges)
          | (null) afl2 = do
              (s,e) &lt;- deWall p1 afl1 box1
              return (s ++ sigma, e ++ edges)
          | otherwise   = do
              x &lt;- get
              liftIO $ do
                  (s1, e1) &lt;- evalStateT (deWall p1 afl1 box1) x
                  (s2, e2) &lt;- evalStateT (deWall p2 afl2 box2) x
                  return (s1 ++ s2 ++ sigma, e1 ++ e2 ++ edges)
  
          where   afl1 = aflBox1 deWallSet
                  afl2 = aflBox2 deWallSet
  </code></pre>
  
  <p>State and IO monads are used to pipe the state and to generate UID for each tetrahedron found using MVar's. My first attempt was to add a forkIO but it doesn't work. It gives a wrong output due a lack of control during the merge part that doesn't wait for both threads to finish. I don't know how to make it wait for them.</p>
  
  <pre><code>            liftIO $ do
                  let 
                      s1 = evalStateT (deWall p1 afl1 box1) x
                      s2 = evalStateT (deWall p2 afl2 box2) x
                      concatThread var (a1, b1) = takeMVar var &gt;&gt;= \(a2, b2) -&gt; putMVar var (a1 ++ a2, b1 ++ b2)
                  mv &lt;- newMVar ([],[])
                  forkIO (s1 &gt;&gt;= concatThread mv)
                  forkIO (s2 &gt;&gt;= concatThread mv)
                  takeMVar mv &gt;&gt;= \(s, e) -&gt; return (s ++ sigma, e ++ edges)
  </code></pre>
  
  <p>So, my next attempt was to use a better parallel strategy "par" and "pseq" which gives the right result but no parallel execution according to threadScope.  </p>
  
  <pre><code>        liftIO $ do
              let
                  s1 = evalStateT (deWall p1 afl1 box1) x
                  s2 = evalStateT (deWall p2 afl2 box2) x
                  conc = liftM2 (\(a1, b1) (a2, b2) -&gt; (a1 ++ a2, b1 ++ b2))
              (stotal, etotal) = s1 `par` (s2 `pseq` (s1 `conc` s2))
              return (stotal ++ sigma, etotal ++ edges)
  </code></pre>
  
  <p>What am I doing wrong?</p>
  
  <p><strong>UPDATE</strong>: Somehow this problem seems to be related with the presence of IO monads. In an other (old) version with no IO monad, only State monad, the parallel execution runs with <code>'par'</code> and <code>'pseq'</code>. The GHC -sstderr gives <code>SPARKS: 1160 (69 converted, 1069 pruned)</code>.</p>
  
  <pre><code>recursion::Box -&gt; Box -&gt; [SimplexPointer] -&gt; [SimplexPointer] -&gt; [Simplex] -&gt; [Edge] -&gt; DeWallSets -&gt; State DeWallSets ([Simplex], [Edge])  
  recursion p1 p2 sigma deWallSet
       | null afl1 &amp;&amp; null afl2 = return sigma
       | (null) afl1 = do
           s &lt;- deWall p2 afl2 box2
           return (s ++ sigma)
       | (null) afl2 = do
           s &lt;- deWall p1 afl1 box1
           return (s ++ sigma)
       | otherwise   = do
                       x &lt;- get
                       let s1 = evalState (deWall p1 afl1 box1) x
                       let s2 = evalState (deWall p2 afl2 box2) x
                       return $ s1 `par` (s2 `pseq` (s1 ++ s2 ++ sigma))
       where   afl1 = aflBox1 deWallSet
               afl2 = aflBox2 deWallSet
  </code></pre>
  
  <p>Cloud someone explain that?</p>