PieterjanRobbe / GaussianRandomFields.jl

A package for Gaussian random field generation in Julia
Other
66 stars 11 forks source link

Use in-place FFTW where possible #22

Closed giordano closed 4 years ago

giordano commented 4 years ago

I'm looking into speeding up _sample(grf::GaussianRandomField{CirculantEmbedding}, xi). A possibility would be to use an in-place FFT instead of the out-of-place version, reusing the y buffer that is already allocated and never used again. I tried this simple change:

diff --git a/src/generators/circulant_embedding.jl b/src/generators/circulant_embedding.jl
index 58904b9..e7c5e74 100644
--- a/src/generators/circulant_embedding.jl
+++ b/src/generators/circulant_embedding.jl
@@ -197,8 +197,8 @@ function _sample(grf::GaussianRandomField{CirculantEmbedding}, xi)
     v, P = grf.data

     # compute multiplication with square root of circulant embedding via FFT
-    y = v .* reshape(xi, size(v))
-    w = P * y
+    w = complex.(v .* reshape(xi, size(v)))
+    mul!(w, P, w)

     # extract realization of random field
     μ, σ = grf.mean, std(grf.cov)

but then tests fail with the error:

ERROR: LoadError: ArgumentError: FFTW out-of-place plan applied to in-place data

I didn't have the time to look into how the plan is created. It would be great to work this out, to cut down allocations and running time.

With this additional change, memory allocations can be brought to 0: define a function

_sample!(z, grf::GaussianRandomField{CirculantEmbedding}, xi)

which is equivalent to the current out-of-place _sample just without definition of z and then the out-of-place _sample would become

function _sample(grf::GaussianRandomField{CirculantEmbedding}, xi)
    z = Array{eltype(grf.cov)}(undef, length.(grf.pts))
    return _sample!(z, grf, xi)
end
PieterjanRobbe commented 4 years ago

That should be easy to fix: just use the in-place version of plan_fft, called plan_fft!? I'll try to make a commit later today.