Skip to content

Commit

Permalink
more fixes suggested by Marco
Browse files Browse the repository at this point in the history
  • Loading branch information
mreineck committed Oct 22, 2024
1 parent ae4f75a commit 2b2dbc2
Show file tree
Hide file tree
Showing 3 changed files with 23 additions and 24 deletions.
9 changes: 9 additions & 0 deletions include/finufft/fft.h
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,9 @@ template<typename T> class Finufft_FFT_plan {
public:
[[maybe_unused]] Finufft_FFT_plan(void (*)(void *) = nullptr,
void (*)(void *) = nullptr, void * = nullptr) {}
// deleting these operations to be consistent with the FFTW plans (seel below)
Finufft_FFT_plan(const Finufft_FFT_plan &) = delete;
Finufft_FFT_plan &operator=(const Finufft_FFT_plan &) = delete;
[[maybe_unused]] void plan(const std::vector<int> & /*dims*/, size_t /*batchSize*/,
std::complex<T> * /*ptr*/, int /*sign*/, int /*options*/,
int /*nthreads*/) {}
Expand Down Expand Up @@ -63,11 +66,15 @@ template<> struct Finufft_FFT_plan<float> {
#endif
unlock();
}
// we have raw pointers in the object (the FFTW plan).
// If we allow copying those, we end up destroying the plans multiple times.
Finufft_FFT_plan(const Finufft_FFT_plan &) = delete;
[[maybe_unused]] ~Finufft_FFT_plan() {
lock();
fftwf_destroy_plan(plan_);
unlock();
}
Finufft_FFT_plan &operator=(const Finufft_FFT_plan &) = delete;

void plan
[[maybe_unused]] (const std::vector<int> &dims, size_t batchSize,
Expand Down Expand Up @@ -131,11 +138,13 @@ template<> struct Finufft_FFT_plan<double> {
#endif
unlock();
}
Finufft_FFT_plan(const Finufft_FFT_plan &) = delete;
[[maybe_unused]] ~Finufft_FFT_plan() {
lock();
fftw_destroy_plan(plan_);
unlock();
}
Finufft_FFT_plan &operator=(const Finufft_FFT_plan &) = delete;

void plan
[[maybe_unused]] (const std::vector<int> &dims, size_t batchSize,
Expand Down
3 changes: 2 additions & 1 deletion include/finufft/finufft_core.h
Original file line number Diff line number Diff line change
Expand Up @@ -201,7 +201,8 @@ template<typename TF> struct FINUFFT_PLAN_T { // the main plan object, fully C++
std::vector<TF> Sp, Tp, Up; // internal primed targs (s'_k, etc),
// allocated
type3params<TF> t3P; // groups together type 3 shift, scale, phase, parameters
FINUFFT_PLAN_T<TF> *innerT2plan = nullptr; // ptr used for type 2 in step 2 of type 3
std::unique_ptr<FINUFFT_PLAN_T<TF>> innerT2plan; // ptr used for type 2 in step 2 of
// type 3

// other internal structs
std::unique_ptr<Finufft_FFT_plan<TF>> fftPlan;
Expand Down
35 changes: 12 additions & 23 deletions src/finufft.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -713,7 +713,6 @@ FINUFFT_PLAN_T<TF>::FINUFFT_PLAN_T(int type_, int dim_, const BIGINT *n_modes, i
stderr,
"[%s] fwBatch would be bigger than MAX_NF, not attempting memory allocation!\n",
__func__);
// FIXME: this error causes memory leaks. We should free phiHat1, phiHat2, phiHat3
throw int(FINUFFT_ERR_MAXNALLOC);
}

Expand Down Expand Up @@ -831,7 +830,7 @@ int FINUFFT_PLAN_T<TF>::setpts(BIGINT nj, TF *xj, TF *yj, TF *zj, BIGINT nk, TF
U = u;

// pick x, s intervals & shifts & # fine grid pts (nf) in each dim...
TF S1, S2, S3; // get half-width X, center C, which contains {x_j}...
TF S1 = 0, S2 = 0, S3 = 0; // get half-width X, center C, which contains {x_j}...
arraywidcen(nj, xj, &(t3P.X1), &(t3P.C1));
arraywidcen(nk, s, &S1, &(t3P.D1)); // same D, S, but for {s_k}
set_nhg_type3(S1, t3P.X1, opts, spopts, &(nf1), &(t3P.h1),
Expand Down Expand Up @@ -928,9 +927,9 @@ int FINUFFT_PLAN_T<TF>::setpts(BIGINT nj, TF *xj, TF *yj, TF *zj, BIGINT nk, TF
}
} else
for (BIGINT j = 0; j < nj; ++j)
prephase[j] = (std::complex<TF>)1.0; // *** or keep flag so no mult in exec??
prephase[j] = {1, 0}; // *** or keep flag so no mult in exec??

// rescale the target s_k etc to s'_k etc...
// rescale the target s_k etc to s'_k etc...
#pragma omp parallel for num_threads(opts.nthreads) schedule(static)
for (BIGINT k = 0; k < nk; ++k) {
Sp[k] = t3P.h1 * t3P.gam1 * (s[k] - t3P.D1); // so |s'_k| < pi/R
Expand All @@ -955,16 +954,9 @@ int FINUFFT_PLAN_T<TF>::setpts(BIGINT nj, TF *xj, TF *yj, TF *zj, BIGINT nk, TF
phiHatk3.resize(nk);
onedim_nuft_kernel(nk, Up, phiHatk3, spopts); // fill phiHat3
}
int Cfinite =
std::isfinite(t3P.C1) && std::isfinite(t3P.C2) && std::isfinite(t3P.C3); // C can
// be nan
// or inf
// if
// M=0,
// no
// input
// NU pts
int Cnonzero = t3P.C1 != 0.0 || t3P.C2 != 0.0 || t3P.C3 != 0.0; // cen
// C can be nan or inf if M=0, no input NU pts
int Cfinite = std::isfinite(t3P.C1) && std::isfinite(t3P.C2) && std::isfinite(t3P.C3);
int Cnonzero = t3P.C1 != 0.0 || t3P.C2 != 0.0 || t3P.C3 != 0.0; // cen
#pragma omp parallel for num_threads(opts.nthreads) schedule(static)
for (BIGINT k = 0; k < nk; ++k) { // .... loop over NU targ freqs
TF phiHat = phiHatk1[k];
Expand Down Expand Up @@ -998,12 +990,11 @@ int FINUFFT_PLAN_T<TF>::setpts(BIGINT nj, TF *xj, TF *yj, TF *zj, BIGINT nk, TF
t2opts.spread_debug = std::max(0, opts.spread_debug - 1);
t2opts.showwarn = 0; // so don't see warnings 2x
// (...could vary other t2opts here?)
if (innerT2plan) {
delete innerT2plan;
innerT2plan = nullptr;
}
int ier = finufft_makeplan_t<TF>(2, d, t2nmodes, fftSign, batchSize, tol,
&innerT2plan, &t2opts);
// MR: temporary hack, until we have figured out the C++ interface.
FINUFFT_PLAN_T<TF> *tmpplan;
int ier = finufft_makeplan_t<TF>(2, d, t2nmodes, fftSign, batchSize, tol, &tmpplan,
&t2opts);
innerT2plan.reset(tmpplan);
if (ier > 1) { // if merely warning, still proceed
fprintf(stderr, "[%s t3]: inner type 2 plan creation failed with ier=%d!\n",
__func__, ier);
Expand Down Expand Up @@ -1181,10 +1172,8 @@ template<typename TF> FINUFFT_PLAN_T<TF>::~FINUFFT_PLAN_T() {
// Also must not crash if called immediately after finufft_makeplan.
// Thus either each thing free'd here is guaranteed to be nullptr or correctly
// allocated.
if (fftPlan) fftPlan->free(fwBatch); // free the big FFTW (or t3 spread) working array
if (fftPlan) fftPlan->free(fwBatch); // free the big FFT (or t3 spread) working array
if (type == 3) {
delete innerT2plan;
innerT2plan = nullptr;
delete[] X;
delete[] Y;
delete[] Z;
Expand Down

0 comments on commit 2b2dbc2

Please sign in to comment.