From 29d9e8760b4bf9cce32219a4c8a5dd76f8e2a6bc Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Wed, 2 Feb 2022 11:26:58 -0600 Subject: [PATCH] fix for out-of-place plans https://github.com/JuliaApproximation/FastTransforms.jl/issues/162 can't use an aliased pointer when the two arrays are genuinely different in the execution --- src/fftw.c | 32 +++++++++++++++++++++++--------- 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/src/fftw.c b/src/fftw.c index 064e79f8..e2c3ad34 100644 --- a/src/fftw.c +++ b/src/fftw.c @@ -99,10 +99,16 @@ ft_sphere_fftw_plan * ft_plan_sph_with_kind(const int N, const int M, const fftw idist = odist = 1; istride = ostride = N; howmany = N; - if (kind[2][0] == FFTW_HC2R) - P->planphi = fftw_plan_many_dft_c2r(rank, n, howmany, (fftw_complex *) P->Y, inembed, istride, idist, P->Y, onembed, ostride, odist, FT_FFTW_FLAGS); - else if (kind[2][0] == FFTW_R2HC) - P->planphi = fftw_plan_many_dft_r2c(rank, n, howmany, P->Y, inembed, istride, idist, (fftw_complex *) P->Y, onembed, ostride, odist, FT_FFTW_FLAGS); + if (kind[2][0] == FFTW_HC2R) { + double * Z = fftw_malloc(N*M*sizeof(double)); + P->planphi = fftw_plan_many_dft_c2r(rank, n, howmany, (fftw_complex *) P->Y, inembed, istride, idist, Z, onembed, ostride, odist, FT_FFTW_FLAGS); + fftw_free(Z); + } + else if (kind[2][0] == FFTW_R2HC) { + double * Z = fftw_malloc(N*M*sizeof(double)); + P->planphi = fftw_plan_many_dft_r2c(rank, n, howmany, Z, inembed, istride, idist, (fftw_complex *) P->Y, onembed, ostride, odist, FT_FFTW_FLAGS); + fftw_free(Z); + } return P; } @@ -436,10 +442,16 @@ ft_disk_fftw_plan * ft_plan_disk_with_kind(const int N, const int M, const fftw_ idist = odist = 1; istride = ostride = N; howmany = N; - if (kind[2][0] == FFTW_HC2R) - P->plantheta = fftw_plan_many_dft_c2r(rank, n, howmany, (fftw_complex *) P->Y, inembed, istride, idist, P->Y, onembed, ostride, odist, FT_FFTW_FLAGS); - else if (kind[2][0] == FFTW_R2HC) - P->plantheta = fftw_plan_many_dft_r2c(rank, n, howmany, P->Y, inembed, istride, idist, (fftw_complex *) P->Y, onembed, ostride, odist, FT_FFTW_FLAGS); + if (kind[2][0] == FFTW_HC2R) { + double * Z = fftw_malloc(N*M*sizeof(double)); + P->plantheta = fftw_plan_many_dft_c2r(rank, n, howmany, (fftw_complex *) P->Y, inembed, istride, idist, Z, onembed, ostride, odist, FT_FFTW_FLAGS); + fftw_free(Z); + } + else if (kind[2][0] == FFTW_R2HC) { + double * Z = fftw_malloc(N*M*sizeof(double)); + P->plantheta = fftw_plan_many_dft_r2c(rank, n, howmany, Z, inembed, istride, idist, (fftw_complex *) P->Y, onembed, ostride, odist, FT_FFTW_FLAGS); + fftw_free(Z); + } return P; } @@ -646,7 +658,9 @@ ft_spinsphere_fftw_plan * ft_plan_spinsph_with_kind(const int N, const int M, co idist = odist = 1; istride = ostride = N; howmany = N; - P->planphi = fftw_plan_many_dft(rank, n, howmany, (fftw_complex *) P->Y, inembed, istride, idist, (fftw_complex *) P->Y, onembed, ostride, odist, sign, FT_FFTW_FLAGS); + fftw_complex * Z = fftw_malloc(N*M*sizeof(fftw_complex)); + P->planphi = fftw_plan_many_dft(rank, n, howmany, (fftw_complex *) P->Y, inembed, istride, idist, Z, onembed, ostride, odist, sign, FT_FFTW_FLAGS); + fftw_free(Z); P->S = S; return P; }