Skip to content

Commit

Permalink
remove finufft.cpp; remove blanket use of namspace std
Browse files Browse the repository at this point in the history
  • Loading branch information
@mreineck
mreineck committed Oct 18, 2024
1 parent 3d102df commit 9374600
Show file tree
Hide file tree
Showing 2 changed files with 29 additions and 64 deletions.
38 changes: 0 additions & 38 deletions src/finufft.cpp
View file

This file was deleted.

55 changes: 29 additions & 26 deletions src/finufft_core.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@
#include <memory>
#include <vector>

using namespace std;
using namespace finufft;
using namespace finufft::utils;
using namespace finufft::spreadinterp;
Expand Down Expand Up @@ -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...
Expand Down Expand Up @@ -209,10 +208,10 @@ static void onedim_fseries_kernel(BIGINT nf, std::vector<T> &fwkerhalf,
a[n] = -exp(2 * PI * std::complex<double>(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<BIGINT> 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<BIGINT> 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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -907,8 +906,8 @@ int FINUFFT_PLAN_T<TF>::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)
Expand Down Expand Up @@ -939,12 +938,16 @@ 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 = 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];
Expand All @@ -955,7 +958,7 @@ int FINUFFT_PLAN_T<TF>::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)
Expand All @@ -971,12 +974,12 @@ int FINUFFT_PLAN_T<TF>::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;
Expand Down Expand Up @@ -1049,7 +1052,7 @@ int FINUFFT_PLAN_T<TF>::execute(std::complex<TF> *cj, std::complex<TF> *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<TF> *cjb = cj + bB * nj; // point to batch of weights
std::complex<TF> *fkb = fk + bB * N; // point to batch of mode coeffs
Expand Down Expand Up @@ -1110,7 +1113,7 @@ int FINUFFT_PLAN_T<TF>::execute(std::complex<TF> *cj, std::complex<TF> *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<TF> *cjb = cj + bB * nj; // batch of input strengths
std::complex<TF> *fkb = fk + bB * nk; // batch of output strengths
Expand Down

0 comments on commit 9374600

Please sign in to comment.