Commit 8a1d38e4 authored by Martin Reinecke's avatar Martin Reinecke

make code a bit more compact

parent fb991029
......@@ -595,7 +595,7 @@ template<typename T0> class cfftp
template<bool fwd, typename T> void pass2 (size_t ido, size_t l1,
const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
const cmplx<T0> * POCKETFFT_RESTRICT wa)
const cmplx<T0> * POCKETFFT_RESTRICT wa) const
{
constexpr size_t cdim=2;
......@@ -644,7 +644,7 @@ template<bool fwd, typename T> void pass2 (size_t ido, size_t l1,
}
template<bool fwd, typename T> void pass3 (size_t ido, size_t l1,
const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
const cmplx<T0> * POCKETFFT_RESTRICT wa)
const cmplx<T0> * POCKETFFT_RESTRICT wa) const
{
constexpr size_t cdim=3;
constexpr T0 tw1r=-0.5,
......@@ -684,7 +684,7 @@ template<bool fwd, typename T> void pass3 (size_t ido, size_t l1,
template<bool fwd, typename T> void pass4 (size_t ido, size_t l1,
const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
const cmplx<T0> * POCKETFFT_RESTRICT wa)
const cmplx<T0> * POCKETFFT_RESTRICT wa) const
{
constexpr size_t cdim=4;
......@@ -760,7 +760,7 @@ template<bool fwd, typename T> void pass4 (size_t ido, size_t l1,
}
template<bool fwd, typename T> void pass5 (size_t ido, size_t l1,
const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
const cmplx<T0> * POCKETFFT_RESTRICT wa)
const cmplx<T0> * POCKETFFT_RESTRICT wa) const
{
constexpr size_t cdim=5;
constexpr T0 tw1r= T0(0.3090169943749474241022934171828191L),
......@@ -832,7 +832,7 @@ template<bool fwd, typename T> void pass5 (size_t ido, size_t l1,
template<bool fwd, typename T> void pass7(size_t ido, size_t l1,
const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
const cmplx<T0> * POCKETFFT_RESTRICT wa)
const cmplx<T0> * POCKETFFT_RESTRICT wa) const
{
constexpr size_t cdim=7;
constexpr T0 tw1r= T0(0.6234898018587335305250048840042398L),
......@@ -881,7 +881,7 @@ template<bool fwd, typename T> void pass7(size_t ido, size_t l1,
#undef POCKETFFT_PARTSTEP7a
#undef POCKETFFT_PREP7
template <bool fwd, typename T> void ROTX45(T &a)
template <bool fwd, typename T> void ROTX45(T &a) const
{
constexpr T0 hsqt2=T0(0.707106781186547524400844362104849L);
if (fwd)
......@@ -889,7 +889,7 @@ template <bool fwd, typename T> void ROTX45(T &a)
else
{ auto tmp_=a.r; a.r=hsqt2*(a.r-a.i); a.i=hsqt2*(a.i+tmp_); }
}
template <bool fwd, typename T> void ROTX135(T &a)
template <bool fwd, typename T> void ROTX135(T &a) const
{
constexpr T0 hsqt2=T0(0.707106781186547524400844362104849L);
if (fwd)
......@@ -900,7 +900,7 @@ template <bool fwd, typename T> void ROTX135(T &a)
template<bool fwd, typename T> void pass8 (size_t ido, size_t l1,
const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
const cmplx<T0> * POCKETFFT_RESTRICT wa)
const cmplx<T0> * POCKETFFT_RESTRICT wa) const
{
constexpr size_t cdim=8;
......@@ -1015,7 +1015,7 @@ template<bool fwd, typename T> void pass8 (size_t ido, size_t l1,
template<bool fwd, typename T> void pass11 (size_t ido, size_t l1,
const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
const cmplx<T0> * POCKETFFT_RESTRICT wa)
const cmplx<T0> * POCKETFFT_RESTRICT wa) const
{
constexpr size_t cdim=11;
constexpr T0 tw1r= T0(0.8412535328311811688618116489193677L),
......@@ -1077,7 +1077,7 @@ template<bool fwd, typename T> void pass11 (size_t ido, size_t l1,
template<bool fwd, typename T> void passg (size_t ido, size_t ip,
size_t l1, T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
const cmplx<T0> * POCKETFFT_RESTRICT wa,
const cmplx<T0> * POCKETFFT_RESTRICT csarr)
const cmplx<T0> * POCKETFFT_RESTRICT csarr) const
{
const size_t cdim=ip;
size_t ipph = (ip+1)/2;
......@@ -1183,7 +1183,7 @@ template<bool fwd, typename T> void passg (size_t ido, size_t ip,
}
}
template<bool fwd, typename T> void pass_all(T c[], T0 fct)
template<bool fwd, typename T> void pass_all(T c[], T0 fct) const
{
if (length==1) { c[0]*=fct; return; }
size_t l1=1;
......@@ -1232,11 +1232,8 @@ template<bool fwd, typename T> void pass_all(T c[], T0 fct)
}
public:
template<typename T> void forward(T c[], T0 fct)
{ pass_all<true>(c, fct); }
template<typename T> void backward(T c[], T0 fct)
{ pass_all<false>(c, fct); }
template<typename T> void exec(T c[], T0 fct, bool fwd) const
{ fwd ? pass_all<true>(c, fct) : pass_all<false>(c, fct); }
private:
POCKETFFT_NOINLINE void factorize()
......@@ -1335,12 +1332,12 @@ template<typename T0> class rfftp
/* (a+ib) = conj(c+id) * (e+if) */
template<typename T1, typename T2, typename T3> inline void MULPM
(T1 &a, T1 &b, T2 c, T2 d, T3 e, T3 f)
(T1 &a, T1 &b, T2 c, T2 d, T3 e, T3 f) const
{ a=c*e+d*f; b=c*f-d*e; }
template<typename T> void radf2 (size_t ido, size_t l1,
const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
const T0 * POCKETFFT_RESTRICT wa)
const T0 * POCKETFFT_RESTRICT wa) const
{
constexpr size_t cdim=2;
......@@ -1379,7 +1376,7 @@ template<typename T> void radf2 (size_t ido, size_t l1,
template<typename T> void radf3(size_t ido, size_t l1,
const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
const T0 * POCKETFFT_RESTRICT wa)
const T0 * POCKETFFT_RESTRICT wa) const
{
constexpr size_t cdim=3;
constexpr T0 taur=-0.5, taui=T0(0.8660254037844386467637231707529362L);
......@@ -1419,7 +1416,7 @@ template<typename T> void radf3(size_t ido, size_t l1,
template<typename T> void radf4(size_t ido, size_t l1,
const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
const T0 * POCKETFFT_RESTRICT wa)
const T0 * POCKETFFT_RESTRICT wa) const
{
constexpr size_t cdim=4;
constexpr T0 hsqt2=T0(0.707106781186547524400844362104849L);
......@@ -1467,7 +1464,7 @@ template<typename T> void radf4(size_t ido, size_t l1,
template<typename T> void radf5(size_t ido, size_t l1,
const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
const T0 * POCKETFFT_RESTRICT wa)
const T0 * POCKETFFT_RESTRICT wa) const
{
constexpr size_t cdim=5;
constexpr T0 tr11= T0(0.3090169943749474241022934171828191L),
......@@ -1524,7 +1521,7 @@ template<typename T> void radf5(size_t ido, size_t l1,
template<typename T> void radfg(size_t ido, size_t ip, size_t l1,
T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
const T0 * POCKETFFT_RESTRICT wa, const T0 * POCKETFFT_RESTRICT csarr)
const T0 * POCKETFFT_RESTRICT wa, const T0 * POCKETFFT_RESTRICT csarr) const
{
const size_t cdim=ip;
size_t ipph=(ip+1)/2;
......@@ -1666,7 +1663,7 @@ template<typename T> void radfg(size_t ido, size_t ip, size_t l1,
template<typename T> void radb2(size_t ido, size_t l1,
const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
const T0 * POCKETFFT_RESTRICT wa)
const T0 * POCKETFFT_RESTRICT wa) const
{
constexpr size_t cdim=2;
......@@ -1698,7 +1695,7 @@ template<typename T> void radb2(size_t ido, size_t l1,
template<typename T> void radb3(size_t ido, size_t l1,
const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
const T0 * POCKETFFT_RESTRICT wa)
const T0 * POCKETFFT_RESTRICT wa) const
{
constexpr size_t cdim=3;
constexpr T0 taur=-0.5, taui=T0(0.8660254037844386467637231707529362L);
......@@ -1739,7 +1736,7 @@ template<typename T> void radb3(size_t ido, size_t l1,
template<typename T> void radb4(size_t ido, size_t l1,
const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
const T0 * POCKETFFT_RESTRICT wa)
const T0 * POCKETFFT_RESTRICT wa) const
{
constexpr size_t cdim=4;
constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
......@@ -1792,7 +1789,7 @@ template<typename T> void radb4(size_t ido, size_t l1,
template<typename T> void radb5(size_t ido, size_t l1,
const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
const T0 * POCKETFFT_RESTRICT wa)
const T0 * POCKETFFT_RESTRICT wa) const
{
constexpr size_t cdim=5;
constexpr T0 tr11= T0(0.3090169943749474241022934171828191L),
......@@ -1852,7 +1849,7 @@ template<typename T> void radb5(size_t ido, size_t l1,
template<typename T> void radbg(size_t ido, size_t ip, size_t l1,
T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
const T0 * POCKETFFT_RESTRICT wa, const T0 * POCKETFFT_RESTRICT csarr)
const T0 * POCKETFFT_RESTRICT wa, const T0 * POCKETFFT_RESTRICT csarr) const
{
const size_t cdim=ip;
size_t ipph=(ip+1)/ 2;
......@@ -1985,7 +1982,7 @@ template<typename T> void radbg(size_t ido, size_t ip, size_t l1,
}
}
template<typename T> void copy_and_norm(T *c, T *p1, size_t n, T0 fct)
template<typename T> void copy_and_norm(T *c, T *p1, size_t n, T0 fct) const
{
if (p1!=c)
{
......@@ -2002,15 +1999,15 @@ template<typename T> void radbg(size_t ido, size_t ip, size_t l1,
}
public:
template<typename T> void forward(T c[], T0 fct)
template<typename T> void exec(T c[], T0 fct, bool r2hc) const
{
if (length==1) { c[0]*=fct; return; }
size_t n=length;
size_t l1=n, nf=fact.size();
size_t n=length, nf=fact.size();
arr<T> ch(n);
T *p1=c, *p2=ch.data();
for(size_t k1=0; k1<nf;++k1)
if (r2hc)
for(size_t k1=0, l1=n; k1<nf;++k1)
{
size_t k=nf-k1-1;
size_t ip=fact[k].fct;
......@@ -2028,18 +2025,8 @@ template<typename T> void radbg(size_t ido, size_t ip, size_t l1,
{ radfg(ido, ip, l1, p1, p2, fact[k].tw, fact[k].tws); swap (p1,p2); }
swap (p1,p2);
}
copy_and_norm(c,p1,n,fct);
}
template<typename T> void backward(T c[], T0 fct)
{
if (length==1) { c[0]*=fct; return; }
size_t n=length;
size_t l1=1, nf=fact.size();
arr<T> ch(n);
T *p1=c, *p2=ch.data();
for(size_t k=0; k<nf; k++)
else
for(size_t k=0, l1=1; k<nf; k++)
{
size_t ip = fact[k].fct,
ido= n/(ip*l1);
......@@ -2056,6 +2043,7 @@ template<typename T> void radbg(size_t ido, size_t ip, size_t l1,
swap (p1,p2);
l1*=ip;
}
copy_and_norm(c,p1,n,fct);
}
......@@ -2153,7 +2141,7 @@ template<typename T0> class fftblue
arr<cmplx<T0>> mem;
cmplx<T0> *bk, *bkf;
template<bool fwd, typename T> void fft(cmplx<T> c[], T0 fct)
template<bool fwd, typename T> void fft(cmplx<T> c[], T0 fct) const
{
arr<cmplx<T>> akf(n2);
......@@ -2164,14 +2152,14 @@ template<typename T0> class fftblue
for (size_t m=n; m<n2; ++m)
akf[m]=zero;
plan.forward (akf.data(),1.);
plan.exec (akf.data(),1.,true);
/* do the convolution */
for (size_t m=0; m<n2; ++m)
akf[m] = akf[m].template special_mul<!fwd>(bkf[m]);
/* inverse FFT */
plan.backward (akf.data(),1.);
plan.exec (akf.data(),1.,false);
/* multiply by b_k */
for (size_t m=0; m<n; ++m)
......@@ -2203,18 +2191,26 @@ template<typename T0> class fftblue
bkf[m] = bkf[n2-m] = bk[m]*xn2;
for (size_t m=n;m<=(n2-n);++m)
bkf[m].Set(0.,0.);
plan.forward(bkf,1.);
plan.exec(bkf,1.,true);
}
template<typename T> void backward(cmplx<T> c[], T0 fct)
{ fft<false>(c,fct); }
template<typename T> void forward(cmplx<T> c[], T0 fct)
{ fft<true>(c,fct); }
template<typename T> void exec(cmplx<T> c[], T0 fct, bool fwd) const
{ fwd ? fft<true>(c,fct) : fft<false>(c,fct); }
template<typename T> void backward_r(T c[], T0 fct)
template<typename T> void exec_r(T c[], T0 fct, bool fwd)
{
arr<cmplx<T>> tmp(n);
if (fwd)
{
auto zero = T0(0)*c[0];
for (size_t m=0; m<n; ++m)
tmp[m].Set(c[m], zero);
fft<true>(tmp.data(),fct);
c[0] = tmp[0].r;
memcpy (c+1, tmp.data()+1, (n-1)*sizeof(T));
}
else
{
tmp[0].Set(c[0],c[0]*0);
memcpy (reinterpret_cast<void *>(tmp.data()+1),
reinterpret_cast<void *>(c+1), (n-1)*sizeof(T));
......@@ -2225,16 +2221,6 @@ template<typename T0> class fftblue
for (size_t m=0; m<n; ++m)
c[m] = tmp[m].r;
}
template<typename T> void forward_r(T c[], T0 fct)
{
arr<cmplx<T>> tmp(n);
auto zero = T0(0)*c[0];
for (size_t m=0; m<n; ++m)
tmp[m].Set(c[m], zero);
fft<true>(tmp.data(),fct);
c[0] = tmp[0].r;
memcpy (c+1, tmp.data()+1, (n-1)*sizeof(T));
}
};
......@@ -2269,11 +2255,8 @@ template<typename T0> class pocketfft_c
packplan=unique_ptr<cfftp<T0>>(new cfftp<T0>(length));
}
template<typename T> POCKETFFT_NOINLINE void backward(cmplx<T> c[], T0 fct) const
{ packplan ? packplan->backward(c,fct) : blueplan->backward(c,fct); }
template<typename T> POCKETFFT_NOINLINE void forward(cmplx<T> c[], T0 fct) const
{ packplan ? packplan->forward(c,fct) : blueplan->forward(c,fct); }
template<typename T> POCKETFFT_NOINLINE void exec(cmplx<T> c[], T0 fct, bool fwd) const
{ packplan ? packplan->exec(c,fct,fwd) : blueplan->exec(c,fct,fwd); }
size_t length() const { return len; }
};
......@@ -2309,17 +2292,8 @@ template<typename T0> class pocketfft_r
packplan=unique_ptr<rfftp<T0>>(new rfftp<T0>(length));
}
template<typename T> POCKETFFT_NOINLINE void backward(T c[], T0 fct) const
{
packplan ? packplan->backward(c,fct)
: blueplan->backward_r(c,fct);
}
template<typename T> POCKETFFT_NOINLINE void forward(T c[], T0 fct) const
{
packplan ? packplan->forward(c,fct)
: blueplan->forward_r(c,fct);
}
template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool fwd) const
{ packplan ? packplan->exec(c,fct,fwd) : blueplan->exec_r(c,fct,fwd); }
size_t length() const { return len; }
};
......@@ -2349,7 +2323,7 @@ template<typename T0> class T_dct1
tmp[0] = c[0];
for (size_t i=1; i<n; ++i)
tmp[i] = tmp[N-i] = c[i];
fftplan.forward(tmp.data(), fct);
fftplan.exec(tmp.data(), fct, true);
c[0] = tmp[0];
for (size_t i=1; i<n; ++i)
c[i] = tmp[2*i-1];
......@@ -2377,7 +2351,7 @@ template<typename T0> class T_dst1
tmp[0] = tmp[n+1] = c[0]*0;
for (size_t i=0; i<n; ++i)
{ tmp[i+1]=c[i]; tmp[N-1-i]=-c[i]; }
fftplan.forward(tmp.data(), fct);
fftplan.exec(tmp.data(), fct, true);
for (size_t i=0; i<n; ++i)
c[i] = -tmp[2*i+2];
}
......@@ -2415,7 +2389,7 @@ template<typename T0> class T_dcst23
if ((N&1)==0) c[N-1]*=2;
for (size_t k=1; k<N-1; k+=2)
MPINPLACE(c[k+1], c[k]);
fftplan.backward(c, fct);
fftplan.exec(c, fct, false);
for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
{
T t1 = twiddle[k-1]*c[kc]+twiddle[kc-1]*c[k];
......@@ -2443,7 +2417,7 @@ template<typename T0> class T_dcst23
}
if ((N&1)==0)
c[NS2] *= 2*twiddle[NS2-1];
fftplan.forward(c, fct);
fftplan.exec(c, fct, true);
for (size_t k=1; k<N-1; k+=2)
MPINPLACE(c[k], c[k+1]);
if (!cosine)
......@@ -2511,7 +2485,7 @@ template<typename T0> class T_dcst4
for (; i<N; ++i, m+=4)
y[i] = c[m];
}
rfft->forward(y.data(), fct);
rfft->exec(y.data(), fct, true);
{
size_t i=0;
for (; i+i+1<n2; ++i)
......@@ -2544,7 +2518,7 @@ template<typename T0> class T_dcst4
y[i].Set(c[2*i],c[N-1-2*i]);
y[i] *= C2[i];
}
fft->forward(y.data(), fct);
fft->exec(y.data(), fct, true);
for(size_t i=0, ic=n2-1; i<n2; ++i, --ic)
{
c[2*i ] = 2*(y[i ].r*C2[i ].r-y[i ].i*C2[i ].i);
......@@ -2901,7 +2875,7 @@ template<typename T> POCKETFFT_NOINLINE void general_c(
tdatav[i].r[j] = tin[it.iofs(j,i)].r;
tdatav[i].i[j] = tin[it.iofs(j,i)].i;
}
forward ? plan->forward (tdatav, fct) : plan->backward(tdatav, fct);
plan->exec (tdatav, fct, forward);
for (size_t i=0; i<len; ++i)
for (size_t j=0; j<vlen; ++j)
out[it.oofs(j,i)].Set(tdatav[i].r[j],tdatav[i].i[j]);
......@@ -2916,7 +2890,7 @@ template<typename T> POCKETFFT_NOINLINE void general_c(
if (buf != &tin[it.iofs(0)])
for (size_t i=0; i<len; ++i)
buf[i] = tin[it.iofs(i)];
forward ? plan->forward (buf, fct) : plan->backward(buf, fct);
plan->exec (buf, fct, forward);
if (buf != &out[it.oofs(0)])
for (size_t i=0; i<len; ++i)
out[it.oofs(i)] = buf[i];
......@@ -2956,7 +2930,7 @@ template<typename T> POCKETFFT_NOINLINE void general_hartley(
for (size_t i=0; i<len; ++i)
for (size_t j=0; j<vlen; ++j)
tdatav[i][j] = tin[it.iofs(j,i)];
plan->forward(tdatav, fct);
plan->exec(tdatav, fct, true);
for (size_t j=0; j<vlen; ++j)
out[it.oofs(j,0)] = tdatav[0][j];
size_t i=1, i1=1, i2=len-1;
......@@ -2977,7 +2951,7 @@ template<typename T> POCKETFFT_NOINLINE void general_hartley(
auto tdata = reinterpret_cast<T *>(storage.data());
for (size_t i=0; i<len; ++i)
tdata[i] = tin[it.iofs(i)];
plan->forward(tdata, fct);
plan->exec(tdata, fct, true);
// Hartley order
out[it.oofs(0)] = tdata[0];
size_t i=1, i1=1, i2=len-1;
......@@ -3072,7 +3046,7 @@ template<typename T> POCKETFFT_NOINLINE void general_r2c(
for (size_t i=0; i<len; ++i)
for (size_t j=0; j<vlen; ++j)
tdatav[i][j] = in[it.iofs(j,i)];
plan->forward(tdatav, fct);
plan->exec(tdatav, fct, true);
for (size_t j=0; j<vlen; ++j)
out[it.oofs(j,0)].Set(tdatav[0][j]);
size_t i=1, ii=1;
......@@ -3095,7 +3069,7 @@ template<typename T> POCKETFFT_NOINLINE void general_r2c(
auto tdata = reinterpret_cast<T *>(storage.data());
for (size_t i=0; i<len; ++i)
tdata[i] = in[it.iofs(i)];
plan->forward(tdata, fct);
plan->exec(tdata, fct, true);
out[it.oofs(0)].Set(tdata[0]);
size_t i=1, ii=1;
if (forward)
......@@ -3145,7 +3119,7 @@ template<typename T> POCKETFFT_NOINLINE void general_c2r(
for (size_t j=0; j<vlen; ++j)
tdatav[i][j] = in[it.iofs(j,ii)].r;
}
plan->backward(tdatav, fct);
plan->exec(tdatav, fct, false);
for (size_t i=0; i<len; ++i)
for (size_t j=0; j<vlen; ++j)
out[it.oofs(j,i)] = tdatav[i][j];
......@@ -3167,7 +3141,7 @@ template<typename T> POCKETFFT_NOINLINE void general_c2r(
if (i<len)
tdata[i] = in[it.iofs(ii)].r;
}
plan->backward(tdata, fct);
plan->exec(tdata, fct, false);
for (size_t i=0; i<len; ++i)
out[it.oofs(i)] = tdata[i];
}
......@@ -3207,8 +3181,7 @@ template<typename T> POCKETFFT_NOINLINE void general_r(
if ((!r2c) && forward)
for (size_t i=2; i<len; i+=2)
tdatav[i] = -tdatav[i];
forward ? plan->forward (tdatav, fct)
: plan->backward(tdatav, fct);
plan->exec (tdatav, fct, forward);
if (r2c && (!forward))
for (size_t i=2; i<len; i+=2)
tdatav[i] = -tdatav[i];
......@@ -3229,7 +3202,7 @@ template<typename T> POCKETFFT_NOINLINE void general_r(
if ((!r2c) && forward)
for (size_t i=2; i<len; i+=2)
buf[i] = -buf[i];
forward ? plan->forward(buf, fct) : plan->backward(buf, fct);
plan->exec(buf, fct, forward);
if (r2c && (!forward))
for (size_t i=2; i<len; i+=2)
buf[i] = -buf[i];
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment