From 6efe49a1efc0d0253e61973876d88fe33baa02bc Mon Sep 17 00:00:00 2001 From: Your Name Date: Thu, 29 Mar 2012 15:54:34 -0400 Subject: [PATCH] Machine precision detection --- src/cic.cpp | 205 +++++++++++++++++++++++++++++++++++++++++++++++++++ src/cic.hpp | 35 +++++++++ src/mach.hpp | 18 +++++ 3 files changed, 258 insertions(+) create mode 100644 src/cic.cpp create mode 100644 src/cic.hpp create mode 100644 src/mach.hpp diff --git a/src/cic.cpp b/src/cic.cpp new file mode 100644 index 0000000..71a106d --- /dev/null +++ b/src/cic.cpp @@ -0,0 +1,205 @@ +#include +#include +#include +#include "cic.hpp" + +CICFilter::CICFilter(uint32_t N, double len) +{ + spatialLen = len; + szGrid = N; + totalSize = N*N*N; + densityGrid = new CICType[totalSize]; + resetMesh(); +} + +CICFilter::~CICFilter() +{ + delete[] densityGrid; +} + +void CICFilter::resetMesh() +{ + for (uint32_t i = 0; i < totalSize; i++) + densityGrid[i] = 0; +} + +void CICFilter::putParticles(CICParticles *particles, uint32_t N) +{ +#if 0 + uint32_t numCorners = 1 << NUMDIMS; + + for (uint32_t i = 0; i < N; i++) + { + Coordinates xyz; + int32_t ixyz[NUMDIMS]; + int32_t rxyz[NUMDIMS]; + CICType alpha[NUMDIMS]; + CICType beta[NUMDIMS]; + for (int j = 0; j < NUMDIMS; j++) + { + xyz[j] = (particles[i].coords[j] / spatialLen * szGrid); + ixyz[j] = (int32_t)floor(xyz[j] - 0.5); + beta[j] = xyz[j] - ixyz[j] - 0.5; + alpha[j] = 1 - beta[j]; + if (ixyz[j] < 0) + ixyz[j] = szGrid-1; + } + + CICType tot_mass = 0; + for (int j = 0; j < numCorners; j++) + { + CICType rel_mass = 1; + uint32_t idx = 0; + uint32_t mul = 1; + uint32_t mul2 = 1; + + for (int k = 0; k < NUMDIMS; k++) + { + uint32_t ipos = ((j & mul2) != 0); + + if (ipos == 1) + { + rel_mass *= beta[k]; + } + else + { + rel_mass *= alpha[k]; + } + + rxyz[k] = ixyz[k] + ipos; + + if (rxyz[k] >= szGrid) + idx += (rxyz[k] - szGrid) * mul; + else + idx += rxyz[k] * mul; + + mul2 *= 2; + mul *= szGrid; + } + + assert(rel_mass > 0); + assert(rel_mass < 1); + assert(idx < totalSize); + densityGrid[idx] += rel_mass * particles[i].mass; + tot_mass += rel_mass; + } + assert(tot_mass < 1.1); + assert(tot_mass > 0.9); + } +#endif +#if 0 + for (uint32_t i = 0; i < N; i++) + { + Coordinates xyz; + int32_t ixyz[NUMDIMS]; + for (int j = 0; j < NUMDIMS; j++) + { + xyz[j] = (particles[i].coords[j] / spatialLen * szGrid); + ixyz[j] = (int32_t)round(xyz[j] - 0.5); + if (ixyz[j] < 0) + ixyz[j] = szGrid-1; + else if (ixyz[j] >= szGrid) + ixyz[j] = 0; + } + + uint32_t idx = ixyz[0] + ixyz[1] * szGrid + ixyz[2] * szGrid * szGrid; + densityGrid[idx] += particles[i].mass; + } + +#endif + + for (uint32_t i = 0; i < N; i++) + { + CICType x, y, z; + int32_t ix, iy, iz; + int32_t ix2, iy2, iz2; + + x = particles[i].coords[0] / spatialLen * szGrid + 0.5; + y = particles[i].coords[1] / spatialLen * szGrid + 0.5; + z = particles[i].coords[2] / spatialLen * szGrid + 0.5; + + if (x < 0) + x += szGrid; + if (y < 0) + y += szGrid; + if (z < 0) + z += szGrid; + + ix = ((int32_t)floor(x)); + iy = ((int32_t)floor(y)); + iz = ((int32_t)floor(z)); + + ix2 = (ix + 1) % szGrid; + iy2 = (iy + 1) % szGrid; + iz2 = (iz + 1) % szGrid; + + CICType alpha_x = x - ix; + CICType alpha_y = y - iy; + CICType alpha_z = z - iz; + + ix %= szGrid; + iy %= szGrid; + iz %= szGrid; + + assert(alpha_x >= 0); + assert(alpha_y >= 0); + assert(alpha_z >= 0); + + CICType beta_x = 1 - alpha_x; + CICType beta_y = 1 - alpha_y; + CICType beta_z = 1 - alpha_z; + + assert(beta_x >= 0); + assert(beta_y >= 0); + assert(beta_z >= 0); + + CICType mass = particles[i].mass; + uint32_t idx; + + // 000 + idx = ix + (iy + iz * szGrid) * szGrid; + densityGrid[idx] += + mass * beta_x * beta_y * beta_z; + + // 100 + idx = ix2 + (iy + iz * szGrid) * szGrid; + densityGrid[idx] += + mass * alpha_x * beta_y * beta_z; + + // 010 + idx = ix + (iy2 + iz * szGrid) * szGrid; + densityGrid[idx] += + mass * beta_x * alpha_y * beta_z; + + // 110 + idx = ix2 + (iy2 + iz * szGrid) * szGrid; + densityGrid[idx] += + mass * alpha_x * alpha_y * beta_z; + + // 001 + idx = ix + (iy + iz2 * szGrid) * szGrid; + densityGrid[idx] += + mass * beta_x * beta_y * alpha_z; + + // 101 + idx = ix2 + (iy + iz2 * szGrid) * szGrid; + densityGrid[idx] += + mass * alpha_x * beta_y * alpha_z; + + // 011 + idx = ix + (iy2 + iz2 * szGrid) * szGrid; + densityGrid[idx] += + mass * beta_x * alpha_y * alpha_z; + + // 111 + idx = ix2 + (iy2 + iz2 * szGrid) * szGrid; + densityGrid[idx] += + mass * alpha_x * alpha_y * alpha_z; + } +} + +void CICFilter::getDensityField(CICType*& field, uint32_t& res) +{ + field = densityGrid; + res = totalSize; +} diff --git a/src/cic.hpp b/src/cic.hpp new file mode 100644 index 0000000..d522d00 --- /dev/null +++ b/src/cic.hpp @@ -0,0 +1,35 @@ +#ifndef __CICFILTER_HPP +#define __CICFILTER_HPP + +#include "CosmoTool/config.hpp" +#include + +using namespace CosmoTool; + +typedef float CICType; + + typedef struct + { + float mass; + Coordinates coords; + } CICParticles; + + class CICFilter + { + public: + CICFilter(uint32_t resolution, double spatialLen); + ~CICFilter(); + + void resetMesh(); + void putParticles(CICParticles *particles, uint32_t N); + + void getDensityField(CICType*& field, uint32_t& res); + + protected: + CICType *densityGrid; + double spatialLen; + uint32_t totalSize; + uint32_t szGrid; + }; + +#endif diff --git a/src/mach.hpp b/src/mach.hpp new file mode 100644 index 0000000..60c98bf --- /dev/null +++ b/src/mach.hpp @@ -0,0 +1,18 @@ +#ifndef __MACHINE_INFORMATION_HPP +#define __MACHINE_INFORMATION_HPP + +template +T mach_epsilon() +{ + T eps = (T)1; + + do + { + eps /= 2; + } + while ((T)(1 + (eps/2)) != (T)1); + + return eps; +} + +#endif