Skip to content

GitLab

P
pypocketfft
Commit b9f99ec0 authored by Martin Reinecke's avatar Martin Reinecke
Browse files
Options

Merge branch 'thread-qual' into 'master' 

Use qualified std::thread to avoid name collisions on AIX

See merge request !31
parents 2c80bd37 56bcb939
Hide whitespace changes
Inline Side-by-side
Showing with 123 additions and 117 deletions
pocketfft_hdronly.h
View file @ b9f99ec0
......@@ -90,11 +90,16 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
namespace pocketfft {
namespace detail {
using std::size_t;
using std::ptrdiff_t;
using namespace std;
// Always use std:: for <cmath> functions
template <typename T> T cos(T) = delete;
template <typename T> T sin(T) = delete;
template <typename T> T sqrt(T) = delete;
using shape_t = vector<size_t>;
using stride_t = vector<ptrdiff_t>;
using shape_t = std::vector<size_t>;
using stride_t = std::vector<ptrdiff_t>;
constexpr bool FORWARD = true,
BACKWARD = false;
......@@ -146,7 +151,7 @@ template<typename T> class arr
{
if (num==0) return nullptr;
void *res = malloc(num*sizeof(T));
if (!res) throw bad_alloc();
if (!res) throw std::bad_alloc();
return reinterpret_cast<T *>(res);
}
static void dealloc(T *ptr)
......@@ -156,7 +161,7 @@ template<typename T> class arr
{
if (num==0) return nullptr;
void *res = aligned_alloc(64,num*sizeof(T));
if (!res) throw bad_alloc();
if (!res) throw std::bad_alloc();
return reinterpret_cast<T *>(res);
}
static void dealloc(T *ptr)
......@@ -166,7 +171,7 @@ template<typename T> class arr
{
if (num==0) return nullptr;
void *ptr = malloc(num*sizeof(T)+64);
if (!ptr) throw bad_alloc();
if (!ptr) throw std::bad_alloc();
T *res = reinterpret_cast<T *>
((reinterpret_cast<size_t>(ptr) & ~(size_t(63))) + 64);
(reinterpret_cast<void**>(res))[-1] = ptr;
......@@ -269,7 +274,7 @@ template<bool fwd, typename T> void ROTX90(cmplx<T> &a)
template<typename T> class sincos_2pibyn
{
private:
using Thigh = typename conditional<(sizeof(T)>sizeof(double)), T, double>::type;
using Thigh = typename std::conditional<(sizeof(T)>sizeof(double)), T, double>::type;
size_t N, mask, shift;
arr<cmplx<Thigh>> v1, v2;
......@@ -280,14 +285,14 @@ template<typename T> class sincos_2pibyn
{
if (x<2*n) // first quadrant
{
if (x<n) return cmplx<Thigh>(cos(Thigh(x)*ang), sin(Thigh(x)*ang));
return cmplx<Thigh>(sin(Thigh(2*n-x)*ang), cos(Thigh(2*n-x)*ang));
if (x<n) return cmplx<Thigh>(std::cos(Thigh(x)*ang), std::sin(Thigh(x)*ang));
return cmplx<Thigh>(std::sin(Thigh(2*n-x)*ang), std::cos(Thigh(2*n-x)*ang));
}
else // second quadrant
{
x-=2*n;
if (x<n) return cmplx<Thigh>(-sin(Thigh(x)*ang), cos(Thigh(x)*ang));
return cmplx<Thigh>(-cos(Thigh(2*n-x)*ang), sin(Thigh(2*n-x)*ang));
if (x<n) return cmplx<Thigh>(-std::sin(Thigh(x)*ang), std::cos(Thigh(x)*ang));
return cmplx<Thigh>(-std::cos(Thigh(2*n-x)*ang), std::sin(Thigh(2*n-x)*ang));
}
}
else
......@@ -295,14 +300,14 @@ template<typename T> class sincos_2pibyn
x=8*n-x;
if (x<2*n) // third quadrant
{
if (x<n) return cmplx<Thigh>(cos(Thigh(x)*ang), -sin(Thigh(x)*ang));
return cmplx<Thigh>(sin(Thigh(2*n-x)*ang), -cos(Thigh(2*n-x)*ang));
if (x<n) return cmplx<Thigh>(std::cos(Thigh(x)*ang), -std::sin(Thigh(x)*ang));
return cmplx<Thigh>(std::sin(Thigh(2*n-x)*ang), -std::cos(Thigh(2*n-x)*ang));
}
else // fourth quadrant
{
x-=2*n;
if (x<n) return cmplx<Thigh>(-sin(Thigh(x)*ang), -cos(Thigh(x)*ang));
return cmplx<Thigh>(-cos(Thigh(2*n-x)*ang), -sin(Thigh(2*n-x)*ang));
if (x<n) return cmplx<Thigh>(-std::sin(Thigh(x)*ang), -std::cos(Thigh(x)*ang));
return cmplx<Thigh>(-std::cos(Thigh(2*n-x)*ang), -std::sin(Thigh(2*n-x)*ang));
}
}
}
......@@ -439,11 +444,11 @@ struct util // hack to avoid duplicate symbols
const stride_t &stride_in, const stride_t &stride_out, bool inplace)
{
auto ndim = shape.size();
if (ndim<1) throw runtime_error("ndim must be >= 1");
if (ndim<1) throw std::runtime_error("ndim must be >= 1");
if ((stride_in.size()!=ndim) || (stride_out.size()!=ndim))
throw runtime_error("stride dimension mismatch");
throw std::runtime_error("stride dimension mismatch");
if (inplace && (stride_in!=stride_out))
throw runtime_error("stride mismatch");
throw std::runtime_error("stride mismatch");
}
static POCKETFFT_NOINLINE void sanity_check(const shape_t &shape,
......@@ -455,8 +460,8 @@ struct util // hack to avoid duplicate symbols
shape_t tmp(ndim,0);
for (auto ax : axes)
{
if (ax>=ndim) throw invalid_argument("bad axis number");
if (++tmp[ax]>1) throw invalid_argument("axis specified repeatedly");
if (ax>=ndim) throw std::invalid_argument("bad axis number");
if (++tmp[ax]>1) throw std::invalid_argument("axis specified repeatedly");
}
}
......@@ -465,7 +470,7 @@ struct util // hack to avoid duplicate symbols
size_t axis)
{
sanity_check(shape, stride_in, stride_out, inplace);
if (axis>=shape.size()) throw invalid_argument("bad axis number");
if (axis>=shape.size()) throw std::invalid_argument("bad axis number");
}
#ifdef POCKETFFT_NO_MULTITHREADING
......@@ -482,8 +487,8 @@ struct util // hack to avoid duplicate symbols
if (shape[axis] < 1000)
parallel /= 4;
size_t max_threads = nthreads == 0 ?
thread::hardware_concurrency() : nthreads;
return max(size_t(1), min(parallel, max_threads));
std::thread::hardware_concurrency() : nthreads;
return std::max(size_t(1), std::min(parallel, max_threads));
}
#endif
};
......@@ -504,14 +509,14 @@ void thread_map(size_t /* nthreads */, Func f)
thread_local size_t thread_id = 0;
thread_local size_t num_threads = 1;
static const size_t max_threads = max(1u, thread::hardware_concurrency());
static const size_t max_threads = std::max(1u, std::thread::hardware_concurrency());
class latch
{
atomic<size_t> num_left_;
mutex mut_;
condition_variable completed_;
using lock_t = unique_lock<mutex>;
std::atomic<size_t> num_left_;
std::mutex mut_;
std::condition_variable completed_;
using lock_t = std::unique_lock<std::mutex>;
public:
latch(size_t n): num_left_(n) {}
......@@ -534,11 +539,11 @@ class latch
template <typename T> class concurrent_queue
{
queue<T> q_;
mutex mut_;
condition_variable item_added_;
std::queue<T> q_;
std::mutex mut_;
std::condition_variable item_added_;
bool shutdown_;
using lock_t = unique_lock<mutex>;
using lock_t = std::unique_lock<std::mutex>;
public:
concurrent_queue(): shutdown_(false) {}
......@@ -548,7 +553,7 @@ template <typename T> class concurrent_queue
{
lock_t lock(mut_);
if (shutdown_)
throw runtime_error("Item added to queue after shutdown");
throw std::runtime_error("Item added to queue after shutdown");
q_.push(move(val));
}
item_added_.notify_one();
......@@ -580,12 +585,12 @@ template <typename T> class concurrent_queue
class thread_pool
{
concurrent_queue<function<void()>> work_queue_;
vector<thread> threads_;
concurrent_queue<std::function<void()>> work_queue_;
std::vector<std::thread> threads_;
void worker_main()
{
function<void()> work;
std::function<void()> work;
while (work_queue_.pop(work))
work();
}
......@@ -595,7 +600,7 @@ class thread_pool
size_t nthreads = threads_.size();
for (size_t i=0; i<nthreads; ++i)
{
try { threads_[i] = thread([this]{ worker_main(); }); }
try { threads_[i] = std::thread([this]{ worker_main(); }); }
catch (...)
{
shutdown();
......@@ -613,7 +618,7 @@ class thread_pool
~thread_pool() { shutdown(); }
void submit(function<void()> work)
void submit(std::function<void()> work)
{
work_queue_.push(move(work));
}
......@@ -637,8 +642,8 @@ thread_pool & get_pool()
{
static thread_pool pool;
#ifdef POCKETFFT_PTHREADS
static once_flag f;
call_once(f,
static std::once_flag f;
std::call_once(f,
[]{
pthread_atfork(
+[]{ get_pool().shutdown(); }, // prepare
......@@ -663,8 +668,8 @@ void thread_map(size_t nthreads, Func f)
auto & pool = get_pool();
latch counter(nthreads);
exception_ptr ex;
mutex ex_mut;
std::exception_ptr ex;
std::mutex ex_mut;
for (size_t i=0; i<nthreads; ++i)
{
pool.submit(
......@@ -674,15 +679,15 @@ void thread_map(size_t nthreads, Func f)
try { f(); }
catch (...)
{
lock_guard<mutex> lock(ex_mut);
ex = current_exception();
std::lock_guard<std::mutex> lock(ex_mut);
ex = std::current_exception();
}
counter.count_down();
});
}
counter.wait();
if (ex)
rethrow_exception(ex);
std::rethrow_exception(ex);
}
#endif
......@@ -704,7 +709,7 @@ template<typename T0> class cfftp
size_t length;
arr<cmplx<T0>> mem;
vector<fctdata> fact;
std::vector<fctdata> fact;
void add_factor(size_t factor)
{ fact.push_back({factor, nullptr, nullptr}); }
......@@ -1318,9 +1323,9 @@ template<bool fwd, typename T> void pass_all(T c[], T0 fct) const
else
{
passg<fwd>(ido, ip, l1, p1, p2, fact[k1].tw, fact[k1].tws);
swap(p1,p2);
std::swap(p1,p2);
}
swap(p1,p2);
std::swap(p1,p2);
l1=l2;
}
if (p1!=c)
......@@ -1354,7 +1359,7 @@ template<bool fwd, typename T> void pass_all(T c[], T0 fct) const
len>>=1;
// factor 2 should be at the front of the factor list
add_factor(2);
swap(fact[0].fct, fact.back().fct);
std::swap(fact[0].fct, fact.back().fct);
}
for (size_t divisor=3; divisor*divisor<=len; divisor+=2)
while ((len%divisor)==0)
......@@ -1407,7 +1412,7 @@ template<bool fwd, typename T> void pass_all(T c[], T0 fct) const
POCKETFFT_NOINLINE cfftp(size_t length_)
: length(length_)
{
if (length==0) throw runtime_error("zero-length FFT requested");
if (length==0) throw std::runtime_error("zero-length FFT requested");
if (length==1) return;
factorize();
mem.resize(twsize());
......@@ -1430,7 +1435,7 @@ template<typename T0> class rfftp
size_t length;
arr<T0> mem;
vector<fctdata> fact;
std::vector<fctdata> fact;
void add_factor(size_t factor)
{ fact.push_back({factor, nullptr, nullptr}); }
......@@ -2117,8 +2122,8 @@ template<typename T> void radbg(size_t ido, size_t ip, size_t l1,
else if(ip==5)
radf5(ido, l1, p1, p2, fact[k].tw);
else
{ radfg(ido, ip, l1, p1, p2, fact[k].tw, fact[k].tws); swap (p1,p2); }
swap (p1,p2);
{ radfg(ido, ip, l1, p1, p2, fact[k].tw, fact[k].tws); std::swap (p1,p2); }
std::swap (p1,p2);
}
else
for(size_t k=0, l1=1; k<nf; k++)
......@@ -2135,7 +2140,7 @@ template<typename T> void radbg(size_t ido, size_t ip, size_t l1,
radb5(ido, l1, p1, p2, fact[k].tw);
else
radbg(ido, ip, l1, p1, p2, fact[k].tw, fact[k].tws);
swap (p1,p2);
std::swap (p1,p2);
l1*=ip;
}
......@@ -2153,7 +2158,7 @@ template<typename T> void radbg(size_t ido, size_t ip, size_t l1,
len>>=1;
// factor 2 should be at the front of the factor list
add_factor(2);
swap(fact[0].fct, fact.back().fct);
std::swap(fact[0].fct, fact.back().fct);
}
for (size_t divisor=3; divisor*divisor<=len; divisor+=2)
while ((len%divisor)==0)
......@@ -2216,7 +2221,7 @@ template<typename T> void radbg(size_t ido, size_t ip, size_t l1,
POCKETFFT_NOINLINE rfftp(size_t length_)
: length(length_)
{
if (length==0) throw runtime_error("zero-length FFT requested");
if (length==0) throw std::runtime_error("zero-length FFT requested");
if (length==1) return;
factorize();
mem.resize(twsize());
......@@ -2334,28 +2339,28 @@ template<typename T0> class fftblue
template<typename T0> class pocketfft_c
{
private:
unique_ptr<cfftp<T0>> packplan;
unique_ptr<fftblue<T0>> blueplan;
std::unique_ptr<cfftp<T0>> packplan;
std::unique_ptr<fftblue<T0>> blueplan;
size_t len;
public:
POCKETFFT_NOINLINE pocketfft_c(size_t length)
: len(length)
{
if (length==0) throw runtime_error("zero-length FFT requested");
if (length==0) throw std::runtime_error("zero-length FFT requested");
size_t tmp = (length<50) ? 0 : util::largest_prime_factor(length);
if (tmp*tmp <= length)
{
packplan=unique_ptr<cfftp<T0>>(new cfftp<T0>(length));
packplan=std::unique_ptr<cfftp<T0>>(new cfftp<T0>(length));
return;
}
double comp1 = util::cost_guess(length);
double comp2 = 2*util::cost_guess(util::good_size_cmplx(2*length-1));
comp2*=1.5; /* fudge factor that appears to give good overall performance */
if (comp2<comp1) // use Bluestein
blueplan=unique_ptr<fftblue<T0>>(new fftblue<T0>(length));
blueplan=std::unique_ptr<fftblue<T0>>(new fftblue<T0>(length));
else
packplan=unique_ptr<cfftp<T0>>(new cfftp<T0>(length));
packplan=std::unique_ptr<cfftp<T0>>(new cfftp<T0>(length));
}
template<typename T> POCKETFFT_NOINLINE void exec(cmplx<T> c[], T0 fct, bool fwd) const
......@@ -2371,28 +2376,28 @@ template<typename T0> class pocketfft_c
template<typename T0> class pocketfft_r
{
private:
unique_ptr<rfftp<T0>> packplan;
unique_ptr<fftblue<T0>> blueplan;
std::unique_ptr<rfftp<T0>> packplan;
std::unique_ptr<fftblue<T0>> blueplan;
size_t len;
public:
POCKETFFT_NOINLINE pocketfft_r(size_t length)
: len(length)
{
if (length==0) throw runtime_error("zero-length FFT requested");
if (length==0) throw std::runtime_error("zero-length FFT requested");
size_t tmp = (length<50) ? 0 : util::largest_prime_factor(length);
if (tmp*tmp <= length)
{
packplan=unique_ptr<rfftp<T0>>(new rfftp<T0>(length));
packplan=std::unique_ptr<rfftp<T0>>(new rfftp<T0>(length));
return;
}
double comp1 = 0.5*util::cost_guess(length);
double comp2 = 2*util::cost_guess(util::good_size_cmplx(2*length-1));
comp2*=1.5; /* fudge factor that appears to give good overall performance */
if (comp2<comp1) // use Bluestein
blueplan=unique_ptr<fftblue<T0>>(new fftblue<T0>(length));
blueplan=std::unique_ptr<fftblue<T0>>(new fftblue<T0>(length));
else
packplan=unique_ptr<rfftp<T0>>(new rfftp<T0>(length));
packplan=std::unique_ptr<rfftp<T0>>(new rfftp<T0>(length));
}
template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool fwd) const
......@@ -2466,7 +2471,7 @@ template<typename T0> class T_dcst23
{
private:
pocketfft_r<T0> fftplan;
vector<T0> twiddle;
std::vector<T0> twiddle;
public:
POCKETFFT_NOINLINE T_dcst23(size_t length)
......@@ -2503,7 +2508,7 @@ template<typename T0> class T_dcst23
c[NS2] *= twiddle[NS2-1];
if (!cosine)
for (size_t k=0, kc=N-1; k<kc; ++k, --kc)
swap(c[k], c[kc]);
std::swap(c[k], c[kc]);
if (ortho) c[0]*=sqrt2*T0(0.5);
}
else
......@@ -2511,7 +2516,7 @@ template<typename T0> class T_dcst23
if (ortho) c[0]*=sqrt2;
if (!cosine)
for (size_t k=0, kc=N-1; k<NS2; ++k, --kc)
swap(c[k], c[kc]);
std::swap(c[k], c[kc]);
for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
{
T t1=c[k]+c[kc], t2=c[k]-c[kc];
......@@ -2536,8 +2541,8 @@ template<typename T0> class T_dcst4
{
private:
size_t N;
unique_ptr<pocketfft_c<T0>> fft;
unique_ptr<pocketfft_r<T0>> rfft;
std::unique_ptr<pocketfft_c<T0>> fft;
std::unique_ptr<pocketfft_r<T0>> rfft;
arr<cmplx<T0>> C2;
public:
......@@ -2561,7 +2566,7 @@ template<typename T0> class T_dcst4
size_t n2 = N/2;
if (!cosine)
for (size_t k=0, kc=N-1; k<n2; ++k, --kc)
swap(c[k], c[kc]);
std::swap(c[k], c[kc]);
if (N&1)
{
// The following code is derived from the FFTW3 function apply_re11()
......@@ -2637,18 +2642,18 @@ template<typename T0> class T_dcst4
// multi-D infrastructure
//
template<typename T> shared_ptr<T> get_plan(size_t length)
template<typename T> std::shared_ptr<T> get_plan(size_t length)
{
#if POCKETFFT_CACHE_SIZE==0
return make_shared<T>(length);
return std::make_shared<T>(length);
#else
constexpr size_t nmax=POCKETFFT_CACHE_SIZE;
static array<shared_ptr<T>, nmax> cache;
static array<size_t, nmax> last_access{{0}};
static std::array<std::shared_ptr<T>, nmax> cache;
static std::array<size_t, nmax> last_access{{0}};
static size_t access_counter = 0;
static mutex mut;
static std::mutex mut;
auto find_in_cache = [&]() -> shared_ptr<T>
auto find_in_cache = [&]() -> std::shared_ptr<T>
{
for (size_t i=0; i<nmax; ++i)
if (cache[i] && (cache[i]->length()==length))
......@@ -2668,13 +2673,13 @@ template<typename T> shared_ptr<T> get_plan(size_t length)
};
{
lock_guard<mutex> lock(mut);
std::lock_guard<std::mutex> lock(mut);
auto p = find_in_cache();
if (p) return p;
}
auto plan = make_shared<T>(length);
auto plan = std::make_shared<T>(length);
{
lock_guard<mutex> lock(mut);
std::lock_guard<std::mutex> lock(mut);
auto p = find_in_cache();
if (p) return p;
......@@ -2762,9 +2767,9 @@ template<size_t N> class multi_iter
{
auto nshares = threading::num_threads;
if (nshares==1) return;
if (nshares==0) throw runtime_error("can't run with zero threads");
if (nshares==0) throw std::runtime_error("can't run with zero threads");
auto myshare = threading::thread_id;
if (myshare>=nshares) throw runtime_error("impossible share requested");
if (myshare>=nshares) throw std::runtime_error("impossible share requested");
size_t nbase = rem/nshares;
size_t additional = rem%nshares;
size_t lo = myshare*nbase + ((myshare<additional) ? myshare : additional);
......@@ -2786,7 +2791,7 @@ template<size_t N> class multi_iter
}
void advance(size_t n)
{
if (rem<n) throw runtime_error("underrun");
if (rem<n) throw std::runtime_error("underrun");
for (size_t i=0; i<n; ++i)
{
p_i[i] = p_ii;
......@@ -2839,8 +2844,8 @@ class rev_iter
private:
shape_t pos;
const arr_info &arr;
vector<char> rev_axis;
vector<char> rev_jump;
std::vector<char> rev_axis;
std::vector<char> rev_jump;
size_t last_axis, last_size;
shape_t shp;
ptrdiff_t p, rp;
......@@ -2998,7 +3003,7 @@ POCKETFFT_NOINLINE void general_nd(const cndarr<T> &in, ndarr<T> &out,
const shape_t &axes, T0 fct, size_t nthreads, const Exec & exec,
const bool allow_inplace=true)
{
shared_ptr<Tplan> plan;
std::shared_ptr<Tplan> plan;
for (size_t iax=0; iax<axes.size(); ++iax)
{
......@@ -3244,7 +3249,7 @@ struct ExecR2R
template<typename T> void c2c(const shape_t &shape, const stride_t &stride_in,
const stride_t &stride_out, const shape_t &axes, bool forward,
const complex<T> *data_in, complex<T> *data_out, T fct,
const std::complex<T> *data_in, std::complex<T> *data_out, T fct,
size_t nthreads=1)
{
if (util::prod(shape)==0) return;
......@@ -3258,7 +3263,7 @@ template<typename T> void dct(const shape_t &shape,
const stride_t &stride_in, const stride_t &stride_out, const shape_t &axes,
int type, const T *data_in, T *data_out, T fct, bool ortho, size_t nthreads=1)
{
if ((type<1) || (type>4)) throw invalid_argument("invalid DCT type");
if ((type<1) || (type>4)) throw std::invalid_argument("invalid DCT type");
if (util::prod(shape)==0) return;
util::sanity_check(shape, stride_in, stride_out, data_in==data_out, axes);
cndarr<T> ain(data_in, shape, stride_in);
......@@ -3276,7 +3281,7 @@ template<typename T> void dst(const shape_t &shape,
const stride_t &stride_in, const stride_t &stride_out, const shape_t &axes,
int type, const T *data_in, T *data_out, T fct, bool ortho, size_t nthreads=1)
{
if ((type<1) || (type>4)) throw invalid_argument("invalid DST type");
if ((type<1) || (type>4)) throw std::invalid_argument("invalid DST type");
if (util::prod(shape)==0) return;
util::sanity_check(shape, stride_in, stride_out, data_in==data_out, axes);
cndarr<T> ain(data_in, shape, stride_in);
......@@ -3292,7 +3297,7 @@ template<typename T> void dst(const shape_t &shape,
template<typename T> void r2c(const shape_t &shape_in,
const stride_t &stride_in, const stride_t &stride_out, size_t axis,
bool forward, const T *data_in, complex<T> *data_out, T fct,
bool forward, const T *data_in, std::complex<T> *data_out, T fct,
size_t nthreads=1)
{
if (util::prod(shape_in)==0) return;
......@@ -3306,7 +3311,7 @@ template<typename T> void r2c(const shape_t &shape_in,
template<typename T> void r2c(const shape_t &shape_in,
const stride_t &stride_in, const stride_t &stride_out, const shape_t &axes,
bool forward, const T *data_in, complex<T> *data_out, T fct,
bool forward, const T *data_in, std::complex<T> *data_out, T fct,
size_t nthreads=1)
{
if (util::prod(shape_in)==0) return;
......@@ -3324,7 +3329,7 @@ template<typename T> void r2c(const shape_t &shape_in,
template<typename T> void c2r(const shape_t &shape_out,
const stride_t &stride_in, const stride_t &stride_out, size_t axis,
bool forward, const complex<T> *data_in, T *data_out, T fct,
bool forward, const std::complex<T> *data_in, T *data_out, T fct,
size_t nthreads=1)
{
if (util::prod(shape_out)==0) return;
......@@ -3338,7 +3343,7 @@ template<typename T> void c2r(const shape_t &shape_out,
template<typename T> void c2r(const shape_t &shape_out,
const stride_t &stride_in, const stride_t &stride_out, const shape_t &axes,
bool forward, const complex<T> *data_in, T *data_out, T fct,
bool forward, const std::complex<T> *data_in, T *data_out, T fct,
size_t nthreads=1)
{
if (util::prod(shape_out)==0) return;
......@@ -3354,7 +3359,7 @@ template<typename T> void c2r(const shape_t &shape_out,
for (int i=int(shape_in.size())-2; i>=0; --i)
stride_inter[size_t(i)] =
stride_inter[size_t(i+1)]*ptrdiff_t(shape_in[size_t(i+1)]);
arr<complex<T>> tmp(nval);
arr<std::complex<T>> tmp(nval);
auto newaxes = shape_t({axes.begin(), --axes.end()});
c2c(shape_in, stride_in, stride_inter, newaxes, forward, data_in, tmp.data(),
T(1), nthreads);
......@@ -3398,9 +3403,9 @@ template<typename T> void r2r_genuine_hartley(const shape_t &shape,
util::sanity_check(shape, stride_in, stride_out, data_in==data_out, axes);
shape_t tshp(shape);
tshp[axes.back()] = tshp[axes.back()]/2+1;
arr<complex<T>> tdata(util::prod(tshp));
arr<std::complex<T>> tdata(util::prod(tshp));
stride_t tstride(shape.size());
tstride.back()=sizeof(complex<T>);
tstride.back()=sizeof(std::complex<T>);
for (size_t i=tstride.size()-1; i>0; --i)
tstride[i-1]=tstride[i]*ptrdiff_t(tshp[i]);
r2c(shape, stride_in, tstride, axes, true, data_in, tdata.data(), fct, nthreads);
......
pypocketfft.cc
View file @ b9f99ec0
......@@ -20,9 +20,10 @@
namespace {
using namespace std;
using pocketfft::shape_t;
using pocketfft::stride_t;
using std::size_t;
using std::ptrdiff_t;
namespace py = pybind11;
......@@ -66,13 +67,13 @@ shape_t makeaxes(const py::array &in, const py::object &axes)
auto tmp=axes.cast<std::vector<ptrdiff_t>>();
auto ndim = in.ndim();
if ((tmp.size()>size_t(ndim)) || (tmp.size()==0))
throw runtime_error("bad axes argument");
throw std::runtime_error("bad axes argument");
for (auto& sz: tmp)
{
if (sz<0)
sz += ndim;
if ((sz>=ndim) || (sz<0))
throw invalid_argument("axes exceeds dimensionality of output");
throw std::invalid_argument("axes exceeds dimensionality of output");
}
return shape_t(tmp.begin(), tmp.end());
}
......@@ -82,7 +83,7 @@ shape_t makeaxes(const py::array &in, const py::object &axes)
if (py::isinstance<py::array_t<T1>>(arr)) return func<double> args; \
if (py::isinstance<py::array_t<T2>>(arr)) return func<float> args; \
if (py::isinstance<py::array_t<T3>>(arr)) return func<ldbl_t> args; \
throw runtime_error("unsupported data type"); \
throw std::runtime_error("unsupported data type"); \
}
template<typename T> T norm_fct(int inorm, size_t N)
......@@ -90,7 +91,7 @@ template<typename T> T norm_fct(int inorm, size_t N)
if (inorm==0) return T(1);
if (inorm==2) return T(1/ldbl_t(N));
if (inorm==1) return T(1/sqrt(ldbl_t(N)));
throw invalid_argument("invalid value for inorm (must be 0, 1, or 2)");
throw std::invalid_argument("invalid value for inorm (must be 0, 1, or 2)");
}
template<typename T> T norm_fct(int inorm, const shape_t &shape,
......@@ -109,7 +110,7 @@ template<typename T> py::array_t<T> prepare_output(py::object &out_,
if (out_.is_none()) return py::array_t<T>(dims);
auto tmp = out_.cast<py::array_t<T>>();
if (!tmp.is(out_)) // a new object was created during casting
throw runtime_error("unexpected data type for output array");
throw std::runtime_error("unexpected data type for output array");
return tmp;
}
......@@ -119,11 +120,11 @@ template<typename T> py::array c2c_internal(const py::array &in,
{
auto axes = makeaxes(in, axes_);