diff --git a/src/finufft.cpp b/src/finufft.cpp deleted file mode 100644 index 758fcb723..000000000 --- a/src/finufft.cpp +++ /dev/null @@ -1,38 +0,0 @@ -// public header -#include - -// private headers for lib build -// (must come after finufft.h which clobbers FINUFFT* macros) -#include - -void FINUFFT_DEFAULT_OPTS(finufft_opts *o) { finufft_default_opts_t(o); } - -int FINUFFT_MAKEPLAN(int type, int dim, const BIGINT *n_modes, int iflag, int ntrans, - FLT tol, FINUFFT_PLAN *pp, finufft_opts *opts) { - return finufft_makeplan_t(type, dim, n_modes, iflag, ntrans, tol, - reinterpret_cast **>(pp), opts); -} - -int FINUFFT_SETPTS(FINUFFT_PLAN p, BIGINT nj, FLT *xj, FLT *yj, FLT *zj, BIGINT nk, - FLT *s, FLT *t, FLT *u) { - return finufft_setpts_t(reinterpret_cast *>(p), nj, xj, yj, zj, - nk, s, t, u); -} - -int FINUFFT_EXECUTE(FINUFFT_PLAN p, CPX *cj, CPX *fk) { - return finufft_execute_t(reinterpret_cast *>(p), cj, fk); -} - -int FINUFFT_DESTROY(FINUFFT_PLAN p) -// Free everything we allocated inside of finufft_plan pointed to by p. -// 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 (!p) // nullptr, so not a ptr to a plan, report error - return 1; - - delete reinterpret_cast *>(p); - p = nullptr; - return 0; // success -} diff --git a/src/finufft_core.cpp b/src/finufft_core.cpp index 834420bbe..60f62f1e5 100644 --- a/src/finufft_core.cpp +++ b/src/finufft_core.cpp @@ -12,7 +12,6 @@ #include #include -using namespace std; using namespace finufft; using namespace finufft::utils; using namespace finufft::spreadinterp; @@ -150,12 +149,12 @@ static void set_nhg_type3(T S, T X, finufft_opts opts, finufft_spread_opts spopt Xsafe = 1.0; Ssafe = 1.0; } else - Xsafe = max(Xsafe, 1 / S); + Xsafe = std::max(Xsafe, 1 / S); else - Ssafe = max(Ssafe, 1 / X); + Ssafe = std::max(Ssafe, 1 / X); // use the safe X and S... auto nfd = T(2.0 * opts.upsampfac * Ssafe * Xsafe / PI + nss); - if (!isfinite(nfd)) nfd = 0.0; // use T to catch inf + if (!std::isfinite(nfd)) nfd = 0.0; // use T to catch inf *nf = (BIGINT)nfd; // printf("initial nf=%lld, ns=%d\n",*nf,spopts.nspread); // catch too small nf, and nan or +-inf, otherwise spread fails... @@ -209,10 +208,10 @@ static void onedim_fseries_kernel(BIGINT nf, std::vector &fwkerhalf, a[n] = -exp(2 * PI * std::complex(0, 1) * z[n] / double(nf)); // phase winding // rates } - BIGINT nout = nf / 2 + 1; // how many values we're writing to - int nt = min(nout, (BIGINT)opts.nthreads); // how many chunks - std::vector brk(nt + 1); // start indices for each thread - for (int t = 0; t <= nt; ++t) // split nout mode indices btw threads + BIGINT nout = nf / 2 + 1; // how many values we're writing to + int nt = std::min(nout, (BIGINT)opts.nthreads); // how many chunks + std::vector brk(nt + 1); // start indices for each thread + for (int t = 0; t <= nt; ++t) // split nout mode indices btw threads brk[t] = (BIGINT)(0.5 + nout * t / (double)nt); #pragma omp parallel num_threads(nt) { // each thread gets own chunk to do @@ -618,7 +617,7 @@ int finufft_makeplan_t(int type, int dim, const BIGINT *n_modes, int iflag, int p->nbatch = 1 + (ntrans - 1) / nthr; // min # batches poss p->batchSize = 1 + (ntrans - 1) / p->nbatch; // then cut # thr in each b } else { // batchSize override by user - p->batchSize = min(p->opts.maxbatchsize, ntrans); + p->batchSize = std::min(p->opts.maxbatchsize, ntrans); p->nbatch = 1 + (ntrans - 1) / p->batchSize; // resulting # batches } if (p->opts.spread_thread == 0) p->opts.spread_thread = 2; // our auto choice @@ -907,8 +906,8 @@ int FINUFFT_PLAN_T::setpts(BIGINT nj, TF *xj, TF *yj, TF *zj, BIGINT nk, TF TF phase = t3P.D1 * xj[j]; if (d > 1) phase += t3P.D2 * yj[j]; if (d > 2) phase += t3P.D3 * zj[j]; - prephase[j] = cos(phase) + imasign * sin(phase); // Euler - // e^{+-i.phase} + prephase[j] = std::cos(phase) + imasign * std::sin(phase); // Euler + // e^{+-i.phase} } } else for (BIGINT j = 0; j < nj; ++j) @@ -939,12 +938,16 @@ int FINUFFT_PLAN_T::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 = isfinite(t3P.C1) && isfinite(t3P.C2) && 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 + 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 #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]; @@ -955,7 +958,7 @@ int FINUFFT_PLAN_T::setpts(BIGINT nj, TF *xj, TF *yj, TF *zj, BIGINT nk, TF TF phase = (s[k] - t3P.D1) * t3P.C1; if (d > 1) phase += (t[k] - t3P.D2) * t3P.C2; if (d > 2) phase += (u[k] - t3P.D3) * t3P.C3; - deconv[k] *= cos(phase) + imasign * sin(phase); // Euler e^{+-i.phase} + deconv[k] *= std::cos(phase) + imasign * std::sin(phase); // Euler e^{+-i.phase} } } if (opts.debug) @@ -971,12 +974,12 @@ int FINUFFT_PLAN_T::setpts(BIGINT nj, TF *xj, TF *yj, TF *zj, BIGINT nk, TF // Plan and setpts once, for the (repeated) inner type 2 finufft call... timer.restart(); - BIGINT t2nmodes[] = {nf1, nf2, nf3}; // t2 input is actually fw - finufft_opts t2opts = opts; // deep copy, since not ptrs - t2opts.modeord = 0; // needed for correct t3! - t2opts.debug = max(0, opts.debug - 1); // don't print as much detail - t2opts.spread_debug = max(0, opts.spread_debug - 1); - t2opts.showwarn = 0; // so don't see warnings 2x + BIGINT t2nmodes[] = {nf1, nf2, nf3}; // t2 input is actually fw + finufft_opts t2opts = opts; // deep copy, since not ptrs + t2opts.modeord = 0; // needed for correct t3! + t2opts.debug = std::max(0, opts.debug - 1); // don't print as much detail + 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; @@ -1049,7 +1052,7 @@ int FINUFFT_PLAN_T::execute(std::complex *cj, std::complex *fk) { for (int b = 0; b * batchSize < ntrans; b++) { // .....loop b over batches // current batch is either batchSize, or possibly truncated if last one - int thisBatchSize = min(ntrans - b * batchSize, batchSize); + int thisBatchSize = std::min(ntrans - b * batchSize, batchSize); int bB = b * batchSize; // index of vector, since batchsizes same std::complex *cjb = cj + bB * nj; // point to batch of weights std::complex *fkb = fk + bB * N; // point to batch of mode coeffs @@ -1110,7 +1113,7 @@ int FINUFFT_PLAN_T::execute(std::complex *cj, std::complex *fk) { for (int b = 0; b * batchSize < ntrans; b++) { // .....loop b over batches // batching and pointers to this batch, identical to t1,2 above... - int thisBatchSize = min(ntrans - b * batchSize, batchSize); + int thisBatchSize = std::min(ntrans - b * batchSize, batchSize); int bB = b * batchSize; std::complex *cjb = cj + bB * nj; // batch of input strengths std::complex *fkb = fk + bB * nk; // batch of output strengths