DiagOp with a repeated FFTOp returns wrong products when run with multiple threads #44

DiagOp((op, op)) returns wrong forward and adjoint products when op is an FFTOp and Julia runs with more than one thread. The same operator gives correct results with one thread, and also when the second block is copy(op).

Reproducer

# JULIA_NUM_THREADS=2 julia repro.jl
using LinearOperatorCollection, FFTW, LinearAlgebra, Random

N = 128
op = FFTOp(ComplexF64; shape = (N, N))
D = DiagOp((op, op))          # DiagOp((op, copy(op))) gives correct results
n = N * N

Random.seed!(1)
bad = 0
for _ in 1:200
  x = rand(ComplexF64, 2n)
  ref = vcat(op * x[1:n], op * x[n+1:end])   # same blocks, applied sequentially
  global bad += norm(D * x - ref) / norm(ref) > 1e-8
end
println("nthreads = $(Threads.nthreads()), mismatches = $bad/200")

With JULIA_NUM_THREADS=2, 20 out of 20 runs had mismatches: 195–199 of 200 forward products and 196–199 of 200 adjoint products (adjoint(D) * x) were wrong. In each run, the largest relative error was between 6e10 and 7e12. With 2 threads and DiagOp((op, copy(op))), 20 out of 20 runs had no mismatches. With 1 thread and DiagOp((op, op)), 20 out of 20 runs also had no mismatches.

Environment: Julia 1.10.12 (julia:1.10-bookworm Docker image, Linux x86_64). LinearOperatorCollection main at 17b5893, FFTW 1.10.0, LinearOperators 2.14.3, OhMyThreads 0.8.9.

Cause

diagOpProd, diagOpTProd and diagOpCTProd run each block's mul! in its own task (@tasks with the default DynamicScheduler()). The FFTOp prod!/ctprod! closures all capture the operator's single tmpVec and apply in-place plans to it. When two blocks are the same FFTOpImpl, two tasks write the same buffer at once.

DiagOp(op, N) copies op for each block (copyOpsFn = copy), so it is not affected. DiagOp(ops) and DiagOp(ops...) keep the children as given. The same applies to any child operator that owns a mutable work buffer.

Possible directions

I'm not sure which behaviour you'd prefer, so here are the options I can see:

1. At construction, detect children that are the same object (===) and copy the repeats, as DiagOp(op, N) does.
2. Run the blocks sequentially when children share an object.
3. Document that DiagOp children must not share stateful operators when the scheduler is parallel. Users can then pass distinct copies or scheduler = SerialScheduler() (from OhMyThreads). With 2 threads, SerialScheduler() also gave 0/200 mismatches on the snippet above.

I'm happy to open a PR for whichever you prefer.