title: Can you parallelise inexact Jacobian calculation in Julia without special arrays?
body: <p>In Julia, you wanted to calculate an inexact Jacobian based on a vector function, f(x), that requires a lot of calculation to evaluate. Evaluation of the Jacobian is obviously quite neatly parallelisable in concept. My question is, can this be done in Julia without resorting to DistributedArray, SharedArray, etc?</p>
  
  <p>For example, suppose you have the code:</p>
  
  <pre><code>function Jacob(f::Function,x)
    eps=1e-7
    delta=eps*eye(length(x))
    J=zeros(length(x),length(x))
    for i=1:length(x)
      J[:,i]=(f(x+delta[:,i])-f(x-delta[:,i]))/2/eps
    end
    J
  end
  </code></pre>
  
  <p>Is it possible to parallelise this in the same way that you might parallelise the summation of 200000000 random coin flips, as per the manual? That is, something equivalent to</p>
  
  <pre><code>nheads = @parallel (+) for i=1:200000000
    int(randbool())
  end
  </code></pre>
  
  <p>I have tried this:</p>
  
  <pre><code>function Jacob(f::Function,x)
    require("testfunc.jl");
    eps=1e-7
    delta=eps*eye(length(x))
    J=zeros(length(x),length(x))
    J=@parallel (+) for i=1:length(x)
      J[:,i]=(f(x+delta[:,i])-f(x-delta[:,i]))/2/eps
      J
    end
    J
  end
  </code></pre>
  
  <p>where "testfunc.jl" is the name of the file where this code, and the definition of f itself, is found. When I tried this, with f simply evaluating to x.^2+cos(x), I was able to get a proper (diagonal) matrix out, but the values did not match those given by the non-parallel code (which I can confirm to be the correct values). Further investigation suggests that the resulting Jacobian has some values being multiplied by 2 or 3, when using julia -p 4.</p>
  
  <p>Is the approach I have described a plausible one (and simply requires tweaking to prevent duplication of evaluations)? If not, is there another method by which I can evaluate the Jacobian without using the more complicated special Array types?</p>
  
  <p>It appears that adding "J=zeros(n,n)" as the first operation inside the parallel for loop corrects this duplication issue. Can the same thing be done without resorting to such brute force clearing of the J array?</p>