diff --git a/include/finufft/finufft_core.h b/include/finufft/finufft_core.h index cdfa4ab33..8f473dc86 100644 --- a/include/finufft/finufft_core.h +++ b/include/finufft/finufft_core.h @@ -171,10 +171,10 @@ template struct FINUFFT_PLAN_T { // the main plan object, fully C++ BIGINT mu; // number of modes in z (3) direction = N3 BIGINT N; // total # modes (prod of above three) - BIGINT nf1; // size of internal fine grid in x (1) direction - BIGINT nf2; // " y (2) - BIGINT nf3; // " z (3) - BIGINT nf; // total # fine grid points (product of the above three) + BIGINT nf1 = 1; // size of internal fine grid in x (1) direction + BIGINT nf2 = 1; // " y (2) + BIGINT nf3 = 1; // " z (3) + BIGINT nf = 1; // total # fine grid points (product of the above three) int fftSign; // sign in exponential for NUFFT defn, guaranteed to be +-1 @@ -182,15 +182,18 @@ template struct FINUFFT_PLAN_T { // the main plan object, fully C++ std::vector phiHat2; // " y-axis. std::vector phiHat3; // " z-axis. - TC *fwBatch = nullptr; // (batches of) fine grid(s) for FFTW to plan - // & act on. Usually the largest working array + TC *fwBatch = nullptr; // (batches of) fine grid(s) for the FFT to plan + // & act on. Usually the largest working array. + // Unfortunately this can't be a vector, since + // in the presence of FFTW it must be allocated + // by fftw_malloc. std::vector sortIndices; // precomputed NU pt permutation, speeds spread/interp bool didSort; // whether binsorting used (false: identity perm used) TF *X = nullptr, *Y = nullptr, *Z = nullptr; // for t1,2: ptr to user-supplied NU pts - // (no new allocs). for t3: allocated as - // "primed" (scaled) src pts x'_j, etc + // (no new allocs). for t3: points to + // "primed" (scaled) Xp, Yp, Zp // type 3 specific TF *S = nullptr, *T = nullptr, *U = nullptr; // pointers to user's target NU pts arrays @@ -198,8 +201,8 @@ template struct FINUFFT_PLAN_T { // the main plan object, fully C++ std::vector prephase; // pre-phase, for all input NU pts std::vector deconv; // reciprocal of kernel FT, phase, all output NU pts std::vector CpBatch; // working array of prephased strengths - std::vector Sp, Tp, Up; // internal primed targs (s'_k, etc), - // allocated + std::vector Xp, Yp, Zp; // internal primed NU points (x'_j, etc) + std::vector Sp, Tp, Up; // internal primed targs (s'_k, etc) type3params t3P; // groups together type 3 shift, scale, phase, parameters std::unique_ptr> innerT2plan; // ptr used for type 2 in step 2 of // type 3 diff --git a/src/finufft_core.cpp b/src/finufft_core.cpp index 2130cfd8a..464ac5446 100644 --- a/src/finufft_core.cpp +++ b/src/finufft_core.cpp @@ -641,18 +641,10 @@ FINUFFT_PLAN_T::FINUFFT_PLAN_T(int type_, int dim_, const BIGINT *n_modes, i if (ier > 1) // proceed if success or warning throw int(ier); - // set others as defaults (or unallocated for arrays)... - X = nullptr; - Y = nullptr; - Z = nullptr; - nf1 = 1; - nf2 = 1; - nf3 = 1; // crucial to leave as 1 for unused dims - // ------------------------ types 1,2: planning needed --------------------- if (type == 1 || type == 2) { - int nthr_fft = nthr; // give FFTW all threads (or use o.spread_thread?) + int nthr_fft = nthr; // give FFT all threads (or use o.spread_thread?) // Note: batchSize not used since might be only 1. spopts.spread_direction = type; @@ -885,18 +877,19 @@ int FINUFFT_PLAN_T::setpts(BIGINT nj, TF *xj, TF *yj, TF *zj, BIGINT nk, TF // printf("fwbatch, cpbatch ptrs: %llx %llx\n",fwBatch,CpBatch); // alloc rescaled NU src pts x'_j (in X etc), rescaled NU targ pts s'_k ... - // FIXME: should use realloc - if (X) delete[] X; - X = new TF[nj]; + // We do this by resizing Xp, Yp, and Zp, and pointing X, Y, Z to their data; + // this avoids any need for explicit cleanup. + Xp.resize(nj); + X = Xp.data(); Sp.resize(nk); if (d > 1) { - if (Y) delete[] Y; - Y = new TF[nj]; + Yp.resize(nj); + Y = Yp.data(); Tp.resize(nk); } if (d > 2) { - if (Z) delete[] Z; - Z = new TF[nj]; + Zp.resize(nj); + Z = Zp.data(); Up.resize(nk); } @@ -1173,11 +1166,6 @@ template FINUFFT_PLAN_T::~FINUFFT_PLAN_T() { // Thus either each thing free'd here is guaranteed to be nullptr or correctly // allocated. if (fftPlan) fftPlan->free(fwBatch); // free the big FFT (or t3 spread) working array - if (type == 3) { - delete[] X; - delete[] Y; - delete[] Z; - } } template FINUFFT_PLAN_T::~FINUFFT_PLAN_T(); template FINUFFT_PLAN_T::~FINUFFT_PLAN_T();