finufft_plan.cpp

#include <finufft/plan.hpp>
#include <finufft/finufft_core.h>
#include <stdexcept>

namespace finufft {

struct plan_nu2u::impl : public FINUFFT_PLAN_T<float> {
    using FINUFFT_PLAN_T<float>::FINUFFT_PLAN_T;
};

plan_nu2u::plan_nu2u(std::span<const int64_t> modes, const plan_opts& opts) {
    finufft_opts iopts;
    finufft_default_opts_t(&iopts);
    iopts.modeord = opts.modeord;
    iopts.spreadinterponly = opts.spreadinterponly;
    iopts.debug = opts.debug;
    iopts.spread_debug = opts.spread_debug;
    iopts.showwarn = opts.showwarn;
    iopts.nthreads = opts.nthreads;
    iopts.spread_sort = opts.spread_sort;
    iopts.spread_kerevalmeth = opts.spread_kerevalmeth;
    iopts.spread_kerpad = opts.spread_kerpad;
    iopts.upsampfac = opts.upsampfac;
    iopts.spread_thread = opts.spread_thread;
    iopts.maxbatchsize = opts.maxbatchsize;
    iopts.spread_nthr_atomic = opts.spread_nthr_atomic;
    iopts.spread_max_sp_size = opts.spread_max_sp_size;

    int ier;
    m_impl = std::make_unique<impl>(1, int(modes.size()), modes.data(), opts.iexp,
        opts.num_trans, float(opts.tol), &iopts, ier);

    if (ier) {
        throw std::runtime_error("finufft::plan_nu2u: error " + std::to_string(ier) +
            " in plan creation");
    }
}

plan_nu2u::~plan_nu2u() = default;

void plan_nu2u::setpts(int64_t nusize, const float* x, const float* y, const float* z) {
    m_impl->setpts(nusize, const_cast<float*>(x), const_cast<float*>(y), const_cast<float*>(z),
        0, nullptr, nullptr, nullptr // unused for type 1
    );
}

void plan_nu2u::execute(const std::complex<float>* nudata, std::complex<float>* out_udata) {
    m_impl->execute(const_cast<std::complex<float>*>(nudata), out_udata);
}

void plan_nu2u::execute_adjoint(std::complex<float>* out_nudata, const std::complex<float>* udata) {
    m_impl->execute_adjoint(out_nudata, const_cast<std::complex<float>*>(udata));
}

} // namespace finufft

finufft_plan.hpp

#pragma once
#include "../finufft_common/defines.h"
#include <cstdint>
#include <span>
#include <memory>
#include <complex>
// custom extension to finufft to make it more C++ friendly

namespace finufft {

struct plan_opts {
    // execution opts...
    int iexp = 1;      // if >=0, uses +i in Fourier exponential, otherwise -i
    int num_trans = 1; // number of transforms to do (default 1)
    double tol = 1e-3; // desired tolerance (default 1e-3)

    // data handling opts...
    int modeord = 0;           // (type 1,2 only): 0 CMCL-style increasing mode order
                               //                  1 FFT-style mode order
    int spreadinterponly = 0;  // (type 1,2 only): 0 do actual NUFFT
                               // 1 only spread (if type 1) or interpolate (type 2)

    // diagnostic opts...
    int debug = 0;        // 0 silent, 1 some timing/debug, or 2 more
    int spread_debug = 0; // spreader: 0 silent, 1 some timing/debug, or 2 tonnes
    int showwarn = 1;     // 0 don't print warnings to stderr, 1 do

    // algorithm performance opts...
    int nthreads = 0;           // number of threads to use, or 0 uses all available
    int spread_sort = 2;        // spreader: 0 don't sort, 1 do, or 2 heuristic choice
    int spread_kerevalmeth = 1; // spreader: 0 exp(sqrt()), 1 Horner piecewise poly (faster)
    int spread_kerpad = 1;      // (exp(sqrt()) only): 0 don't pad kernel to 4n, 1 do
    double upsampfac = 0.0;     // upsampling ratio sigma: 2.0 std, 1.25 small FFT, 0.0 auto
    int spread_thread = 0;      // (vectorized ntr>1 only): 0 auto, 1 seq multithreaded,
                                //                          2 parallel single-thread spread
    int maxbatchsize = 0;       // (vectorized ntr>1 only): max transform batch, 0 auto
    int spread_nthr_atomic = -1;// if >=0, threads above which spreader OMP critical goes
                                // atomic
    int spread_max_sp_size = 0; // if >0, overrides spreader (dir=1) max subproblem size
                                // sphinx tag (don't remove): @opts_end

};

// non-uniform to uniform (type 1) plan
class FINUFFT_EXPORT plan_nu2u {
public:
    plan_nu2u(std::span<const int64_t> modes, const plan_opts& opts = {});
    ~plan_nu2u();

    void setpts(
        int64_t nusize,
        // x, y, and z, if present, must point to bufs of size nusize
        const float* x,
        const float* y = nullptr,
        const float* z = nullptr
    );

    void setpts(
        std::span<const float> x, // determines nusize
        // y and z, if present, must point to bufs of the same size as x
        const float* y = nullptr,
        const float* z = nullptr
    ) {
        setpts(x.size(), x.data(), y, z);
    }

    void setpts(
        std::span<const float> x, // determines nusize
        // y and z, if present, must be same size as x
        std::span<const float> y,
        std::span<const float> z = {}
    ) {
        setpts(x.size(), x.data(), y.data(), z.data());
    }

    // nudata must point to a buf of size ntrans*nusize,
    // udata must point to a buf of size ntrans*prod(modes)
    void execute(const std::complex<float>* nudata, std::complex<float>* out_udata);
    void execute_adjoint(std::complex<float>* out_nudata, const std::complex<float>* udata);

private:
    struct impl;
    std::unique_ptr<impl> m_impl;
};

} // namespace finufft