bug in parallel_scan ?

[ghennequin]

I'm having issues with parallel_scan which behaves oddly for a particular associative operator -- in particular, it doesn't return the same thing as a fold_left.

open Base
open Domainslib

(* associative operator corresponding to the linear dynamics x_(t+1) = 0.5*x_t + u_t *)
let op (a, u) (b, v) = Float.(a * b, b*u + v)

(* input timeseries that will give us the impulse response *) 
let inputs = Array.init 4 ~f:(function 0 -> 1. | _ -> 0.)

(* assume time invariant decay factor of 0.5 *)
let c = Array.map inputs ~f:(fun x -> (0.5, x))

(* compute the impulse response with a parallel_scan *)
let pool = Task.setup_pool ~num_domains:3 ()
let y = Task.run pool (fun () -> Task.parallel_scan pool op c) |> Array.map ~f:snd
(* this returns [|1.; 1.; 0.5; 1.|], which is wrong; 
    setting num_domains=0 yields the correct result though *)

(* compute the impulse response with a fold *)
let y' = 
     let identity = (1., 0.) in
     Array.fold c ~init:([], identity) ~f:(fun (accu, c) x -> 
          let x' = op c x in 
          x' :: accu, x')
    |> fst |> List.rev |> Array.of_list |> Array.map ~f:snd
(* this returns the correct impulse response [|1.; 0.5; 0.25; 0.125|] *)

parallel_scan seems to be operating correctly on symmetric operators, e.g. sum. I'll try and look into the code for parallel_scan to see if perhaps the arguments of op got swapped somewhere they shouldn't.

[polytypic]

Apologies, I haven’t tried the code sample or looked at it deeply, but I noticed the operation uses floating point and with floating point operations generally cannot be reordered without potentially changing the results slightly due to limited floating point precision. Could that be the case here?