fix for out-of-place plans

JuliaApproximation/FastTransforms.jl#162 can't use an aliased pointer when the two arrays are genuinely different in the execution
MikaelSlevinsky · Feb 2, 2022 · 29d9e87 · 29d9e87
1 parent 4f0beb3
commit 29d9e87
Showing 1 changed file with 23 additions and 9 deletions.
diff --git 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;
 }