easy-iso

commit b85945a0b4bbab8569a22a797112583eab9091b8
parent ca56afe397349b569b4bc8b01fdf6865597eacde
Author: Fabian Wermelinger <fabianw@mavt.ethz.ch>
Date:   Mon, 28 Nov 2016 19:36:53 +0100

fixed integer overflow for VERY large datasets

Diffstat:
M
Makefile.config
 | 
6
++++--
M
include/Matrix3D.h
 | 
7
++++---
M
src/Interpolator.h
 | 
39
++++++++++++---------------------------

3 files changed, 20 insertions(+), 32 deletions(-)
diff --git a/Makefile.config b/Makefile.config
@@ -47,8 +47,10 @@ endif
 
 ifeq "$(hdf)" "1"
 ifneq "$(findstring kilo,$(shell hostname))" ""
-	INC += -I${HOME}/local/hdf5-1.8.17/build/include
-	LIB += -L${HOME}/local/hdf5-1.8.17/build/lib -lhdf5
+	# INC += -I${HOME}/local/hdf5-1.8.17/build/include
+	# LIB += -L${HOME}/local/hdf5-1.8.17/build/lib -lhdf5
+	INC += -I${HOME}/local/hdf5-1.8.17/build-parallel/include
+	LIB += -L${HOME}/local/hdf5-1.8.17/build-parallel/lib -lhdf5
 else
 	INC += -I/usr/local/include
 	LIB += -lhdf5
diff --git a/include/Matrix3D.h b/include/Matrix3D.h
@@ -121,7 +121,8 @@ private:
     int m_Nx;
     int m_Ny;
     int m_Nz;
-    int m_XS, m_XE, m_Xpitch, m_Npitched;
+    int m_XS, m_XE, m_Xpitch;
+    size_t m_Npitched;
     bool m_allocated;
     T* m_data;
 
@@ -134,7 +135,7 @@ private:
         m_XS = alignedElements * (XS - (alignedElements-1))/alignedElements;
         m_XE = alignedElements * (XE + (alignedElements-1))/alignedElements;
         m_Xpitch = (m_Nx+m_XE) - m_XS;
-        m_Npitched = m_Xpitch * (m_Ny+YE-YS) * (m_Nz+ZE-ZS);
+        m_Npitched = static_cast<size_t>(m_Xpitch) * (m_Ny+YE-YS) * (m_Nz+ZE-ZS);
 
         void* pmem;
         if (posix_memalign(&pmem, _ALIGN_, m_Npitched*sizeof(T)))
@@ -149,7 +150,7 @@ private:
 
     inline void _copy(const T* const base)
     {
-        const int Nslice = m_Xpitch * (m_Ny+YE-YS);
+        const size_t Nslice = static_cast<size_t>(m_Xpitch) * (m_Ny+YE-YS);
         for (int iz = 0; iz < m_Nz+ZE-ZS; ++iz)
         {
             const T* const src = base + Nslice*iz;
diff --git a/src/Interpolator.h b/src/Interpolator.h
@@ -17,9 +17,6 @@
 #include <hdf5.h>
 #endif /* _USE_HDF_ */
 
-// TODO: (fabianw@mavt.ethz.ch; Thu Nov 24 08:21:14 2016) testing
-#include <fstream>
-
 #include "common.h"
 #include "ArgumentParser.h"
 #ifdef _USE_CUBISMZ_
@@ -31,6 +28,7 @@ template <typename TKernel>
 class Interpolator
 {
 public:
+    Interpolator() = default;
     Interpolator(ArgumentParser& p) :
         m_extent(p("extent").asDouble(1.0))
     {
@@ -105,24 +103,8 @@ public:
     int  getNz() const { return m_data.Nz(); }
     int  isNodal() const { return m_isNodal; }
 
-    // TODO: (fabianw@mavt.ethz.ch; Thu Nov 24 16:38:31 2016) trash this
-    // void dump_data(const Real x, const Real y, const Real z) const
-    // {
-    //     const Real xo = x - m_xorigin_off;
-    //     const Real yo = y - m_yorigin_off;
-    //     const Real zo = z - m_zorigin_off;
-
-    //     const int iy = _idx(yo);
-    //     const int iz = _idx(zo);
-
-    //     ofstream idata("raw_data.dat");
-    //     for (int ix = 0; ix < m_data.Nx(); ++ix)
-    //         idata << _pos(ix)+m_xorigin_off << "\t" << m_data(ix,iy,iz) << std::endl;
-    // }
-
-
-private:
-    const Real m_extent;
+protected:
+    Real m_extent;
     Real m_xorigin_off;
     Real m_yorigin_off;
     Real m_zorigin_off;
@@ -132,14 +114,15 @@ private:
     Matrix_t m_data;
     TKernel m_kernel;
 
-    // loader
-    void _load_h5(ArgumentParser& p);
-    void _load_wavelet(ArgumentParser& p);
-
     // helper
     std::function<Real(const Real, const int)> _getH;
     inline int _idx(const Real x) const { return x*m_invh; }
     inline Real _pos(const int i) const { return i*m_h; }
+
+private:
+    // loader
+    void _load_h5(ArgumentParser& p);
+    void _load_wavelet(ArgumentParser& p);
 };
 
 
@@ -233,7 +216,8 @@ void Interpolator<TKernel>::_load_wavelet(ArgumentParser& p)
     const int NBY = myreader.yblocks();
     const int NBZ = myreader.zblocks();
     const int maxDim[3] = {NBX*_BLOCKSIZE_, NBY*_BLOCKSIZE_, NBZ*_BLOCKSIZE_};
-    Real* const data = new Real[maxDim[0]*maxDim[1]*maxDim[2]];
+    const size_t num_elem = static_cast<size_t>(maxDim[0]) * maxDim[1] * maxDim[2];
+    Real* const data = new Real[num_elem];
 
     Real blockdata[_BLOCKSIZE_][_BLOCKSIZE_][_BLOCKSIZE_];
     for (int iz=0; iz < NBZ; ++iz)
@@ -247,7 +231,8 @@ void Interpolator<TKernel>::_load_wavelet(ArgumentParser& p)
                     {
                         assert(iy*_BLOCKSIZE_+y < maxDim[1]);
                         assert(iz*_BLOCKSIZE_+z < maxDim[2]);
-                        Real* const dst = data + _BLOCKSIZE_*(ix + NBX*(y+iy*_BLOCKSIZE_ + (z+iz*_BLOCKSIZE_)*NBY*_BLOCKSIZE_));
+                        const size_t offset = _BLOCKSIZE_*(static_cast<size_t>(ix) + NBX*(y+iy*_BLOCKSIZE_ + (z+static_cast<size_t>(iz)*_BLOCKSIZE_)*NBY*_BLOCKSIZE_));
+                        Real* const dst = data + offset;
                         memcpy(dst, &blockdata[z][y][0], _BLOCKSIZE_*sizeof(Real));
                     }
             }