|
| 1 | +#include <finufft.h> |
| 2 | +#include <finufft_common/constants.h> |
| 3 | +#include <finufft_opts.h> |
| 4 | + |
| 5 | +#include <algorithm> |
| 6 | +#include <cmath> |
| 7 | +#include <complex> |
| 8 | +#include <cstdint> |
| 9 | +#include <cstdio> |
| 10 | +#include <thread> |
| 11 | +#include <vector> |
| 12 | + |
| 13 | +#include "utils/dirft1d.hpp" |
| 14 | +#include "utils/norms.hpp" |
| 15 | + |
| 16 | +int main() { |
| 17 | + constexpr int nthreads = 4; |
| 18 | + constexpr int nreps = 16; |
| 19 | + constexpr int M = 400; |
| 20 | + constexpr int64_t N1 = 2048; |
| 21 | + constexpr double tol = 1e-12; |
| 22 | + |
| 23 | + finufft_opts opts; |
| 24 | + finufft_default_opts(&opts); |
| 25 | + opts.nthreads = 1; // crucial: parallelism is across concurrent plan executes |
| 26 | + opts.debug = 0; |
| 27 | + |
| 28 | + std::vector<double> x(M); |
| 29 | + std::vector<std::complex<double>> c(M), ref(N1); |
| 30 | + for (int j = 0; j < M; ++j) { |
| 31 | + double t = static_cast<double>(j) / M; |
| 32 | + x[j] = -finufft::common::PI + 2.0 * finufft::common::PI * t; |
| 33 | + c[j] = std::complex<double>(0.5 * std::cos(13.0 * t) + 0.25 * std::sin(7.0 * t), |
| 34 | + 0.75 * std::sin(11.0 * t) - 0.2 * std::cos(5.0 * t)); |
| 35 | + } |
| 36 | + |
| 37 | + int64_t Ns[3] = {N1, 1, 1}; |
| 38 | + finufft_plan plan; |
| 39 | + int ier = finufft_makeplan(1, 1, Ns, +1, 1, tol, &plan, &opts); |
| 40 | + if (ier != 0) { |
| 41 | + std::fprintf(stderr, "finufft_makeplan failed: ier=%d\n", ier); |
| 42 | + return ier; |
| 43 | + } |
| 44 | + ier = finufft_setpts(plan, M, x.data(), nullptr, nullptr, 0, nullptr, nullptr, nullptr); |
| 45 | + if (ier != 0) { |
| 46 | + std::fprintf(stderr, "finufft_setpts failed: ier=%d\n", ier); |
| 47 | + finufft_destroy(plan); |
| 48 | + return ier; |
| 49 | + } |
| 50 | + |
| 51 | + dirft1d1<int64_t>(M, x, c, +1, N1, ref); |
| 52 | + |
| 53 | + std::vector<int> failures(nthreads, 0); |
| 54 | + |
| 55 | + std::vector<std::thread> workers; |
| 56 | + workers.reserve(nthreads); |
| 57 | + for (int tid = 0; tid < nthreads; ++tid) { |
| 58 | + workers.emplace_back([&, tid]() { |
| 59 | + std::vector<std::complex<double>> out(N1); |
| 60 | + for (int rep = 0; rep < nreps; ++rep) { |
| 61 | + int local_ier = finufft_execute(plan, c.data(), out.data()); |
| 62 | + double relerr = relerrtwonorm(N1, ref.data(), out.data()); |
| 63 | + if (local_ier != 0 || relerr > 10.0 * tol) { |
| 64 | + failures[tid] = 1; |
| 65 | + std::fprintf(stderr, "thread %d rep %d failed: ier=%d relerr=%.3g\n", tid, rep, |
| 66 | + local_ier, relerr); |
| 67 | + return; |
| 68 | + } |
| 69 | + } |
| 70 | + }); |
| 71 | + } |
| 72 | + |
| 73 | + for (auto &worker : workers) worker.join(); |
| 74 | + |
| 75 | + finufft_destroy(plan); |
| 76 | + return *std::max_element(failures.begin(), failures.end()); |
| 77 | +} |
0 commit comments