Commit 7db2f080 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

cleanup

parent c3ac9bf3
......@@ -97,9 +97,11 @@ class GL_Integrator
if (n&1) x[0] = 0.; // set to exact zero
}
template<typename Func> double integrate(Func f)
template<typename Func> auto integrate(Func f) -> decltype(f(0.))
{
double res=0, istart=0;
using T = decltype(f(0.));
T res=0;
size_t istart=0;
if (n_&1)
{
res = f(x[0])*w[0];
......
......@@ -79,12 +79,6 @@ template<typename T, size_t ndim> class mav
{ return d; }
bool last_contiguous() const
{ return (str[ndim-1]==1) || (str[ndim-1]==0); }
void check_storage(const char *name) const
{
if (!last_contiguous())
cout << "Array '" << name << "': last dimension is not contiguous.\n"
"This may slow down computation significantly!\n";
}
bool contiguous() const
{
ptrdiff_t stride=1;
......
......@@ -224,7 +224,8 @@ class Distribution
nthreads_ = 1;
thread_map(move(f));
}
void execStatic(size_t nwork, size_t nthreads, size_t chunksize, std::function<void(Scheduler &)> f)
void execStatic(size_t nwork, size_t nthreads, size_t chunksize,
std::function<void(Scheduler &)> f)
{
mode = STATIC;
nthreads_ = (nthreads==0) ? get_default_nthreads() : nthreads;
......@@ -237,8 +238,8 @@ class Distribution
nextstart[i] = i*chunksize_;
thread_map(move(f));
}
void execDynamic(size_t nwork,
size_t nthreads, size_t chunksize_min, double fact_max, std::function<void(Scheduler &)> f)
void execDynamic(size_t nwork, size_t nthreads, size_t chunksize_min,
double fact_max, std::function<void(Scheduler &)> f)
{
mode = DYNAMIC;
nthreads_ = (nthreads==0) ? get_default_nthreads() : nthreads;
......@@ -392,7 +393,8 @@ void execSingle(size_t nwork, std::function<void(Scheduler &)> func)
MyScheduler sched(nwork);
func(sched);
}
void execStatic(size_t nwork, size_t, size_t, std::function<void(Scheduler &)> func)
void execStatic(size_t nwork, size_t, size_t,
std::function<void(Scheduler &)> func)
{
MyScheduler sched(nwork);
func(sched);
......
......@@ -16,17 +16,15 @@
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
/* Copyright (C) 2019 Max-Planck-Society
/* Copyright (C) 2019-2020 Max-Planck-Society
Author: Martin Reinecke */
#ifndef MRUTIL_UNITY_ROOTS_H
#define MRUTIL_UNITY_ROOTS_H
#include <complex>
#include <cmath>
#include <type_traits>
#include <vector>
#include <mr_util/cmplx.h>
namespace mr {
......@@ -34,28 +32,29 @@ namespace detail_unity_roots {
using namespace std;
template<typename T, typename Tc=complex<T>> class UnityRoots
template<typename T, typename Tc> class UnityRoots
{
private:
using Thigh = typename conditional<(sizeof(T)>sizeof(double)), T, double>::type;
struct cmplx_ { Thigh r, i; };
size_t N, mask, shift;
vector<Cmplx<Thigh>> v1, v2;
vector<cmplx_> v1, v2;
static Cmplx<Thigh> calc(size_t x, size_t n, Thigh ang)
static cmplx_ calc(size_t x, size_t n, Thigh ang)
{
x<<=3;
if (x<4*n) // first half
{
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 {cos(Thigh(x)*ang), sin(Thigh(x)*ang)};
return {sin(Thigh(2*n-x)*ang), 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 {-sin(Thigh(x)*ang), cos(Thigh(x)*ang)};
return {-cos(Thigh(2*n-x)*ang), sin(Thigh(2*n-x)*ang)};
}
}
else
......@@ -63,14 +62,14 @@ template<typename T, typename Tc=complex<T>> class UnityRoots
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 {cos(Thigh(x)*ang), -sin(Thigh(x)*ang)};
return {sin(Thigh(2*n-x)*ang), -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 {-sin(Thigh(x)*ang), -cos(Thigh(x)*ang)};
return {-cos(Thigh(2*n-x)*ang), -sin(Thigh(2*n-x)*ang)};
}
}
}
......@@ -86,11 +85,11 @@ template<typename T, typename Tc=complex<T>> class UnityRoots
while((size_t(1)<<shift)*(size_t(1)<<shift) < nval) ++shift;
mask = (size_t(1)<<shift)-1;
v1.resize(mask+1);
v1[0].Set(Thigh(1), Thigh(0));
v1[0]={Thigh(1), Thigh(0)};
for (size_t i=1; i<v1.size(); ++i)
v1[i]=calc(i,n,ang);
v2.resize((nval+mask)/(mask+1));
v2[0].Set(Thigh(1), Thigh(0));
v2[0]={Thigh(1), Thigh(0)};
for (size_t i=1; i<v2.size(); ++i)
v2[i]=calc(i*(mask+1),n,ang);
}
......
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