easy-iso

Matrix3D.h (5273B)

---

 // File       : Matrix3D.h
 // Date       : Wed Nov 23 14:45:22 2016
 // Author     : Fabian Wermelinger
 // Description: 3D Matrix
 // Copyright 2016 ETH Zurich. All Rights Reserved.
 #ifndef MATRIX3D_H_L9MVKQGE
 #define MATRIX3D_H_L9MVKQGE
 
 #include <cassert>
 #include <cstdlib>
 #include <cstddef>
 #include <cstring>
 #include <iostream>
 
 template <typename T, int XS=0, int XE=0, int YS=0, int YE=0, int ZS=0, int ZE=0>
 class Matrix3D
 {
 public:
     static constexpr int XSTART = XS;
     static constexpr int XEND   = XE;
     static constexpr int YSTART = YS;
     static constexpr int YEND   = YE;
     static constexpr int ZSTART = ZS;
     static constexpr int ZEND   = ZE;
 
     Matrix3D() : m_Nx(0), m_Ny(0), m_Nz(0), m_allocated(false), m_data(nullptr) { }
 
     Matrix3D(const int Nx, const int Ny, const int Nz, const T* const data = nullptr) :
         m_Nx(0), m_Ny(0), m_Nz(0), m_allocated(false), m_data(nullptr)
     {
         alloc(Nx, Ny, Nz, data);
     }
 
     Matrix3D(const Matrix3D& rhs) :
         m_Nx(0), m_Ny(0), m_Nz(0), m_allocated(false), m_data(nullptr)
     {
         if (rhs.m_allocated)
         {
             m_data = _alloc(rhs.m_Nx, rhs.m_Ny, rhs.m_Nz);
             _copy(rhs.m_data);
         }
     }
 
     Matrix3D(Matrix3D&& rhs) :
         m_Nx(rhs.m_Nx), m_Ny(rhs.m_Ny), m_Nz(rhs.m_Nz),
         m_XS(rhs.m_XS), m_XE(rhs.m_XE), m_Xpitch(rhs.m_Xpitch), m_Npitched(rhs.m_Npitched),
         m_allocated(rhs.m_allocated), m_data(std::move(rhs.m_data))
     { }
 
     ~Matrix3D() { if (m_allocated) free(m_data); }
 
     Matrix3D& operator=(const Matrix3D& rhs)
     {
         if (this != &rhs)
         {
             _reset();
             if (rhs.m_allocated)
             {
                 m_data = _alloc(rhs.m_Nx, rhs.m_Ny, rhs.m_Nz);
                 _copy(rhs.m_data);
             }
         }
         return *this;
     }
 
     Matrix3D& operator=(Matrix3D&& rhs)
     {
         if (this != &rhs)
         {
             _reset();
             if (rhs.m_allocated)
             {
                 m_Nx = rhs.m_Nx;
                 m_Ny = rhs.m_Ny;
                 m_Nz = rhs.m_Nz;
                 m_XS = rhs.m_XS;
                 m_XE = rhs.m_XE;
                 m_Xpitch = rhs.m_Xpitch;
                 m_Npitched = rhs.m_Npitched;
                 m_allocated = rhs.m_allocated;
                 m_data = std::move(rhs.m_data);
                 rhs.m_allocated = false;
             }
         }
         return *this;
     }
 
     T operator()(const int ix, const int iy, const int iz) const
     {
         assert(m_allocated);
         assert(ix >= XS && ix < m_Nx+XE);
         assert(iy >= YS && iy < m_Ny+YE);
         assert(iz >= ZS && iz < m_Nz+ZE);
         return m_data[(ix-XS) + m_Xpitch*(iy-YS) + m_Xpitch*(m_Ny+YE-YS)*(iz-ZS)];
     }
 
     T& operator()(const int ix, const int iy, const int iz)
     {
         assert(m_allocated);
         assert(ix >= XS && ix < m_Nx+XE);
         assert(iy >= YS && iy < m_Ny+YE);
         assert(iz >= ZS && iz < m_Nz+ZE);
         return m_data[(ix-XS) + m_Xpitch*(iy-YS) + m_Xpitch*(m_Ny+YE-YS)*(iz-ZS)];
     }
 
     inline void alloc(const int Nx, const int Ny, const int Nz, const T* const data = nullptr)
     {
         if (m_allocated) free(m_data);
         m_data = _alloc(Nx, Ny, Nz);
 
         if (data)
         {
             for (int iz = 0; iz < m_Nz; ++iz)
                 for (int iy = 0; iy < m_Ny; ++iy)
                     {
                         const T* const src = data + m_Nx*iy + m_Nx*m_Ny*iz;
                         T* dst = &(this->operator()(0,iy,iz));
                         memcpy(dst, src, m_Nx*sizeof(T));
                     }
         }
     }
 
     inline int Nx() const { return m_Nx; }
     inline int Ny() const { return m_Ny; }
     inline int Nz() const { return m_Nz; }
 
 private:
     int m_Nx;
     int m_Ny;
     int m_Nz;
     int m_XS, m_XE, m_Xpitch;
     size_t m_Npitched;
     bool m_allocated;
     T* m_data;
 
     inline T* _alloc(const int Nx, const int Ny, const int Nz)
     {
         m_Nx = Nx;
         m_Ny = Ny;
         m_Nz = Nz;
         const int alignedElements = _ALIGN_/sizeof(T);
         m_XS = alignedElements * ((XS - (alignedElements-1))/alignedElements);
         m_XE = alignedElements * ((XE + (alignedElements-1))/alignedElements);
         m_Xpitch = m_Nx+m_XE-m_XS;
         m_Npitched = static_cast<size_t>(m_Xpitch) * static_cast<size_t>(m_Ny+YE-YS) * static_cast<size_t>(m_Nz+ZE-ZS);
 
         void* pmem;
         if (posix_memalign(&pmem, _ALIGN_, m_Npitched*sizeof(T)))
         {
             std::cerr << "ERROR: Matrix3D: can not allocate memory" << std::endl;
             abort();
         }
         memset(pmem, 0, m_Npitched*sizeof(T));
         m_allocated = true;
         return (T*)pmem;
     }
 
     inline void _copy(const T* const base)
     {
         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;
             T* dst = m_data + Nslice*iz;
             memcpy(dst, src, Nslice*sizeof(T));
         }
     }
 
     inline void _reset()
     {
         if (m_allocated) free(m_data);
         m_Nx = m_Ny = m_Nz = m_XS = m_XE = m_Xpitch = m_Npitched = 0;
         m_allocated = false;
         m_data = nullptr;
     }
 };
 
 #endif /* MATRIX3D_H_L9MVKQGE */