VecGeom/base/SOA3D.h

0001 /// \file SOA3D.h
0002 /// \author Johannes de Fine Licht (johannes.definelicht@cern.ch)
0003
0004 #ifndef VECGEOM_BASE_SOA3D_H_
0005 #define VECGEOM_BASE_SOA3D_H_
0006
0007 #include "VecGeom/base/Cuda.h"
0008 #include "VecGeom/base/Global.h"
0009
0010 #include "VecGeom/base/Container3D.h"
0011 #include "VecGeom/backend/scalar/Backend.h"
0012 #ifdef VECGEOM_CUDA_INTERFACE
0013 #include "VecGeom/backend/cuda/Interface.h"
0014 #endif
0015
0016 #include <fstream>
0017
0018 namespace vecgeom {
0019
0020 VECGEOM_DEVICE_FORWARD_DECLARE(template <typename Type> class SOA3D;);
0021
0022 inline namespace VECGEOM_IMPL_NAMESPACE {
0023
0024 // gcc 4.8.2's -Wnon-virtual-dtor is broken and turned on by -Weffc++, we
0025 // need to disable it for SOA3D
0026
0027 #if __GNUC__ < 3 || (__GNUC__ == 4 && __GNUC_MINOR__ <= 8)
0028
0029 #pragma GCC diagnostic push
0030 #pragma GCC diagnostic ignored "-Wnon-virtual-dtor"
0031 #pragma GCC diagnostic ignored "-Weffc++"
0032 #define GCC_DIAG_POP_NEEDED
0033 #endif
0034
0035 template <typename T>
0036 class SOA3D : Container3D<SOA3D<T>> {
0037
0038 private:
0039   bool fAllocated = false;
0040   size_t fSize = 0, fCapacity = 0;
0041   T *fX = nullptr, *fY = nullptr, *fZ = nullptr;
0042
0043 public:
0044   typedef T value_type;
0045
0046   VECCORE_ATT_HOST_DEVICE
0047   SOA3D(T *x, T *y, T *z, size_t size);
0048
0049   VECCORE_ATT_HOST_DEVICE
0050   SOA3D(size_t size);
0051
0052   /// @brief Construct aligned data using specialized allocator.
0053   /// @details To allocate the full object, call this constructor with placement new.
0054   /// @param size Size of the internal arrays
0055   /// @param a Aligned allocator, pre-initialized to fit the content
0056   VECCORE_ATT_HOST_DEVICE
0057   SOA3D(size_t size, AlignedAllocator &a);
0058
0059   SOA3D(SOA3D<T> const &other);
0060
0061   SOA3D() = default;
0062
0063   VECCORE_ATT_HOST_DEVICE
0064   SOA3D &operator=(SOA3D<T> const &other);
0065
0066   VECCORE_ATT_HOST_DEVICE
0067   ~SOA3D();
0068
0069   /// @brief Compute size of a buffer to hold the aligned arrays fX, fY and fZ
0070   /// @param initSize Number of elements
0071   /// @return Size to allocate
0072   VECCORE_ATT_HOST_DEVICE
0073   VECGEOM_FORCE_INLINE
0074   static size_t aligned_sizeof_data(size_t initSize);
0075
0076   VECCORE_ATT_HOST_DEVICE
0077   VECGEOM_FORCE_INLINE
0078   size_t size() const;
0079
0080   VECCORE_ATT_HOST_DEVICE
0081   VECGEOM_FORCE_INLINE
0082   size_t capacity() const;
0083
0084   VECCORE_ATT_HOST_DEVICE
0085   VECGEOM_FORCE_INLINE
0086   void resize(size_t newSize);
0087
0088   VECCORE_ATT_HOST_DEVICE
0089   VECGEOM_FORCE_INLINE
0090   void reserve(size_t newCapacity);
0091
0092   VECGEOM_FORCE_INLINE
0093   void clear();
0094
0095   // Element access methods
0096
0097   VECCORE_ATT_HOST_DEVICE
0098   VECGEOM_FORCE_INLINE
0099   Vector3D<T> operator[](size_t index) const;
0100
0101   VECCORE_ATT_HOST_DEVICE
0102   VECGEOM_FORCE_INLINE
0103   T x(size_t index) const;
0104
0105   VECCORE_ATT_HOST_DEVICE
0106   VECGEOM_FORCE_INLINE
0107   T &x(size_t index);
0108
0109   VECCORE_ATT_HOST_DEVICE
0110   VECGEOM_FORCE_INLINE
0111   T *x();
0112
0113   VECCORE_ATT_HOST_DEVICE
0114   VECGEOM_FORCE_INLINE
0115   T const *x() const;
0116
0117   VECCORE_ATT_HOST_DEVICE
0118   VECGEOM_FORCE_INLINE
0119   T y(size_t index) const;
0120
0121   VECCORE_ATT_HOST_DEVICE
0122   VECGEOM_FORCE_INLINE
0123   T &y(size_t index);
0124
0125   VECCORE_ATT_HOST_DEVICE
0126   VECGEOM_FORCE_INLINE
0127   T *y();
0128
0129   VECCORE_ATT_HOST_DEVICE
0130   VECGEOM_FORCE_INLINE
0131   T const *y() const;
0132
0133   VECCORE_ATT_HOST_DEVICE
0134   VECGEOM_FORCE_INLINE
0135   T z(size_t index) const;
0136
0137   VECCORE_ATT_HOST_DEVICE
0138   VECGEOM_FORCE_INLINE
0139   T &z(size_t index);
0140
0141   VECCORE_ATT_HOST_DEVICE
0142   VECGEOM_FORCE_INLINE
0143   T *z();
0144
0145   VECCORE_ATT_HOST_DEVICE
0146   VECGEOM_FORCE_INLINE
0147   T const *z() const;
0148
0149   // Element manipulation methods
0150
0151   VECCORE_ATT_HOST_DEVICE
0152   VECGEOM_FORCE_INLINE
0153   void set(size_t index, T x, T y, T z);
0154
0155   VECCORE_ATT_HOST_DEVICE
0156   VECGEOM_FORCE_INLINE
0157   void set(size_t index, Vector3D<T> const &vec);
0158
0159   VECCORE_ATT_HOST_DEVICE
0160   VECGEOM_FORCE_INLINE
0161   void push_back(T x, T y, T z);
0162
0163   VECCORE_ATT_HOST_DEVICE
0164   VECGEOM_FORCE_INLINE
0165   void push_back(Vector3D<T> const &vec);
0166
0167 #ifdef VECGEOM_CUDA_INTERFACE
0168   DevicePtr<cuda::SOA3D<T>> CopyToGpu(DevicePtr<T> xGpu, DevicePtr<T> yGpu, DevicePtr<T> zGpu) const;
0169   DevicePtr<cuda::SOA3D<T>> CopyToGpu(DevicePtr<T> xGpu, DevicePtr<T> yGpu, DevicePtr<T> zGpu, size_t size) const;
0170 #endif // VECGEOM_CUDA_INTERFACE
0171
0172   void ToFile(std::string /*filename*/) const;
0173   int FromFile(std::string /*filename*/);
0174
0175 private:
0176   VECCORE_ATT_HOST_DEVICE
0177   void Allocate();
0178
0179   VECCORE_ATT_HOST_DEVICE
0180   void Deallocate();
0181 };
0182
0183 #if defined(GCC_DIAG_POP_NEEDED)
0184
0185 #pragma GCC diagnostic pop
0186 #undef GCC_DIAG_POP_NEEDED
0187
0188 #endif
0189
0190 template <typename T>
0191 VECCORE_ATT_HOST_DEVICE SOA3D<T>::SOA3D(T *xval, T *yval, T *zval, size_t sz)
0192     : fAllocated(false), fSize(sz), fCapacity(fSize), fX(xval), fY(yval), fZ(zval)
0193 {
0194 }
0195
0196 template <typename T>
0197 VECCORE_ATT_HOST_DEVICE SOA3D<T>::SOA3D(size_t initSize, AlignedAllocator &a)
0198     : fAllocated(false), fSize(initSize), fCapacity(initSize)
0199 {
0200   fX = a.aligned_alloc<T>(initSize, kAlignmentBoundary);
0201   fY = a.aligned_alloc<T>(initSize, kAlignmentBoundary);
0202   fZ = a.aligned_alloc<T>(initSize, kAlignmentBoundary);
0203 }
0204
0205 template <typename T>
0206 VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE size_t SOA3D<T>::aligned_sizeof_data(size_t initSize)
0207 {
0208   // Size for separate aligned arrays fX, fY and fZ
0209   return 3 * AlignedAllocator::aligned_sizeof<T>(initSize, kAlignmentBoundary);
0210 }
0211
0212 template <typename T>
0213 VECCORE_ATT_HOST_DEVICE SOA3D<T>::SOA3D(size_t sz) : fSize(sz), fCapacity(sz)
0214 {
0215   Allocate();
0216 }
0217
0218 template <typename T>
0219 SOA3D<T>::SOA3D(SOA3D<T> const &rhs) : fAllocated(false), fSize(rhs.fSize), fCapacity(rhs.fCapacity)
0220 {
0221   if (rhs.fAllocated) {
0222     Allocate();
0223     copy(rhs.fX, rhs.fX + rhs.fSize, fX);
0224     copy(rhs.fY, rhs.fY + rhs.fSize, fY);
0225     copy(rhs.fZ, rhs.fZ + rhs.fSize, fZ);
0226   } else {
0227     fX = rhs.fX;
0228     fY = rhs.fY;
0229     fZ = rhs.fZ;
0230   }
0231 }
0232
0233 template <typename T>
0234 VECCORE_ATT_HOST_DEVICE SOA3D<T> &SOA3D<T>::operator=(SOA3D<T> const &rhs)
0235 {
0236 #ifndef VECCORE_CUDA_DEVICE_COMPILATION
0237   fSize     = rhs.fSize;
0238   fCapacity = rhs.fCapacity;
0239   Deallocate();
0240   if (rhs.fAllocated && rhs.fSize > 0) {
0241     Allocate();
0242     copy(rhs.fX, rhs.fX + rhs.fSize, fX);
0243     copy(rhs.fY, rhs.fY + rhs.fSize, fY);
0244     copy(rhs.fZ, rhs.fZ + rhs.fSize, fZ);
0245   } else {
0246     fX = rhs.fX;
0247     fY = rhs.fY;
0248     fZ = rhs.fZ;
0249   }
0250 #else
0251   fAllocated = false;
0252   fSize      = rhs.fSize;
0253   fCapacity  = rhs.fCapacity;
0254   fX         = rhs.fX;
0255   fY         = rhs.fY;
0256   fZ         = rhs.fZ;
0257 #endif
0258   return *this;
0259 }
0260
0261 template <typename T>
0262 VECCORE_ATT_HOST_DEVICE SOA3D<T>::~SOA3D()
0263 {
0264 #ifndef VECCORE_CUDA_DEVICE_COMPILATION
0265   Deallocate();
0266 #endif
0267 }
0268
0269 template <typename T>
0270 VECCORE_ATT_HOST_DEVICE size_t SOA3D<T>::size() const
0271 {
0272   return fSize;
0273 }
0274
0275 template <typename T>
0276 VECCORE_ATT_HOST_DEVICE size_t SOA3D<T>::capacity() const
0277 {
0278   return fCapacity;
0279 }
0280
0281 template <typename T>
0282 VECCORE_ATT_HOST_DEVICE void SOA3D<T>::resize(size_t newSize)
0283 {
0284   VECGEOM_ASSERT(newSize <= fCapacity);
0285   fSize = newSize;
0286 }
0287
0288 template <typename T>
0289 VECCORE_ATT_HOST_DEVICE void SOA3D<T>::reserve(size_t newCapacity)
0290 {
0291   fCapacity = newCapacity;
0292   T *xNew, *yNew, *zNew;
0293   xNew  = AlignedAllocate<T>(fCapacity);
0294   yNew  = AlignedAllocate<T>(fCapacity);
0295   zNew  = AlignedAllocate<T>(fCapacity);
0296   fSize = (fSize > fCapacity) ? fCapacity : fSize;
0297   if (fX && fY && fZ) {
0298     copy(fX, fX + fSize, xNew);
0299     copy(fY, fY + fSize, yNew);
0300     copy(fZ, fZ + fSize, zNew);
0301   }
0302   Deallocate();
0303   fX         = xNew;
0304   fY         = yNew;
0305   fZ         = zNew;
0306   fAllocated = true;
0307 }
0308
0309 template <typename T>
0310 void SOA3D<T>::clear()
0311 {
0312   Deallocate();
0313   fSize     = 0;
0314   fCapacity = 0;
0315 }
0316
0317 template <typename T>
0318 VECCORE_ATT_HOST_DEVICE void SOA3D<T>::Allocate()
0319 {
0320   if (fCapacity == 0) return;
0321
0322   fX         = AlignedAllocate<T>(fCapacity);
0323   fY         = AlignedAllocate<T>(fCapacity);
0324   fZ         = AlignedAllocate<T>(fCapacity);
0325   fAllocated = true;
0326 }
0327
0328 template <typename T>
0329 VECCORE_ATT_HOST_DEVICE void SOA3D<T>::Deallocate()
0330 {
0331   if (fAllocated) {
0332     AlignedFree(fX);
0333     AlignedFree(fY);
0334     AlignedFree(fZ);
0335   }
0336   fAllocated = false;
0337 }
0338
0339 template <typename T>
0340 VECCORE_ATT_HOST_DEVICE Vector3D<T> SOA3D<T>::operator[](size_t index) const
0341 {
0342   return Vector3D<T>(fX[index], fY[index], fZ[index]);
0343 }
0344
0345 template <typename T>
0346 VECCORE_ATT_HOST_DEVICE T SOA3D<T>::x(size_t index) const
0347 {
0348   return fX[index];
0349 }
0350
0351 template <typename T>
0352 VECCORE_ATT_HOST_DEVICE T &SOA3D<T>::x(size_t index)
0353 {
0354   return fX[index];
0355 }
0356
0357 template <typename T>
0358 VECCORE_ATT_HOST_DEVICE T *SOA3D<T>::x()
0359 {
0360   return fX;
0361 }
0362
0363 template <typename T>
0364 VECCORE_ATT_HOST_DEVICE T const *SOA3D<T>::x() const
0365 {
0366   return fX;
0367 }
0368
0369 template <typename T>
0370 VECCORE_ATT_HOST_DEVICE T SOA3D<T>::y(size_t index) const
0371 {
0372   return fY[index];
0373 }
0374
0375 template <typename T>
0376 VECCORE_ATT_HOST_DEVICE T &SOA3D<T>::y(size_t index)
0377 {
0378   return fY[index];
0379 }
0380
0381 template <typename T>
0382 VECCORE_ATT_HOST_DEVICE T *SOA3D<T>::y()
0383 {
0384   return fY;
0385 }
0386
0387 template <typename T>
0388 VECCORE_ATT_HOST_DEVICE T const *SOA3D<T>::y() const
0389 {
0390   return fY;
0391 }
0392
0393 template <typename T>
0394 VECCORE_ATT_HOST_DEVICE T SOA3D<T>::z(size_t index) const
0395 {
0396   return fZ[index];
0397 }
0398
0399 template <typename T>
0400 VECCORE_ATT_HOST_DEVICE T &SOA3D<T>::z(size_t index)
0401 {
0402   return fZ[index];
0403 }
0404
0405 template <typename T>
0406 VECCORE_ATT_HOST_DEVICE T *SOA3D<T>::z()
0407 {
0408   return fZ;
0409 }
0410
0411 template <typename T>
0412 VECCORE_ATT_HOST_DEVICE T const *SOA3D<T>::z() const
0413 {
0414   return fZ;
0415 }
0416
0417 template <typename T>
0418 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void SOA3D<T>::set(size_t index, T xval, T yval, T zval)
0419 {
0420 // not asserting in case of NVCC -- still getting annoying
0421 // errors on CUDA < 8.0
0422 #ifndef VECCORE_CUDA
0423   VECGEOM_ASSERT(index < fCapacity);
0424 #endif
0425   fX[index] = xval;
0426   fY[index] = yval;
0427   fZ[index] = zval;
0428 }
0429
0430 template <typename T>
0431 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void SOA3D<T>::set(size_t index, Vector3D<T> const &vec)
0432 {
0433 // not asserting in case of NVCC -- still getting annoying
0434 // errors on CUDA < 8.0
0435 #ifndef VECCORE_CUDA
0436   VECGEOM_ASSERT(index < fCapacity);
0437 #endif
0438   fX[index] = vec[0];
0439   fY[index] = vec[1];
0440   fZ[index] = vec[2];
0441 }
0442
0443 template <typename T>
0444 VECCORE_ATT_HOST_DEVICE void SOA3D<T>::push_back(T xval, T yval, T zval)
0445 {
0446   fX[fSize] = xval;
0447   fY[fSize] = yval;
0448   fZ[fSize] = zval;
0449   ++fSize;
0450 }
0451
0452 template <typename T>
0453 VECCORE_ATT_HOST_DEVICE void SOA3D<T>::push_back(Vector3D<T> const &vec)
0454 {
0455   push_back(vec[0], vec[1], vec[2]);
0456 }
0457
0458 template <typename T>
0459 void SOA3D<T>::ToFile(std::string filename) const
0460 {
0461   std::ofstream outfile(filename, std::ios::binary);
0462   outfile.write(reinterpret_cast<const char *>(&fSize), sizeof(fSize));
0463   outfile.write(reinterpret_cast<const char *>(&fCapacity), sizeof(fCapacity));
0464   outfile.write(reinterpret_cast<char *>(fX), fCapacity * sizeof(T));
0465   outfile.write(reinterpret_cast<char *>(fY), fCapacity * sizeof(T));
0466   outfile.write(reinterpret_cast<char *>(fZ), fCapacity * sizeof(T));
0467 }
0468
0469 // read SOA3D from file serialized with ToFile
0470 // returns number of elements read or -1 if failure
0471 template <typename T>
0472 int SOA3D<T>::FromFile(std::string filename)
0473 {
0474   // should be called on an already allocated object
0475   decltype(fSize) s;
0476   decltype(fCapacity) cap;
0477   std::ifstream fin(filename, std::ios::binary);
0478   fin.read(reinterpret_cast<char *>(&s), sizeof(s));
0479   if (!fin) return -1;
0480   fin.read(reinterpret_cast<char *>(&cap), sizeof(cap));
0481   if (!fin) return -2;
0482   //  if (cap != fCapacity || s != fSize)
0483   //    std::cerr << " warning: reading from SOA3D with different size\n";
0484
0485   fin.read(reinterpret_cast<char *>(fX), fCapacity * sizeof(T));
0486   if (!fin) return -3;
0487
0488   fin.read(reinterpret_cast<char *>(fY), fCapacity * sizeof(T));
0489   if (!fin) return -4;
0490
0491   fin.read(reinterpret_cast<char *>(fZ), fCapacity * sizeof(T));
0492   if (!fin) return -5;
0493
0494   return fCapacity;
0495 }
0496
0497 #ifdef VECGEOM_CUDA_INTERFACE
0498
0499 template <typename T>
0500 DevicePtr<cuda::SOA3D<T>> SOA3D<T>::CopyToGpu(DevicePtr<T> xGpu, DevicePtr<T> yGpu, DevicePtr<T> zGpu) const
0501 {
0502   xGpu.ToDevice(fX, fSize);
0503   yGpu.ToDevice(fY, fSize);
0504   zGpu.ToDevice(fZ, fSize);
0505
0506   DevicePtr<cuda::SOA3D<T>> gpu_ptr;
0507   gpu_ptr.Allocate();
0508   gpu_ptr.Construct(xGpu, yGpu, zGpu, fSize);
0509 }
0510
0511 template <typename T>
0512 DevicePtr<cuda::SOA3D<T>> SOA3D<T>::CopyToGpu(DevicePtr<T> xGpu, DevicePtr<T> yGpu, DevicePtr<T> zGpu,
0513                                               size_t count) const
0514 {
0515   xGpu.ToDevice(fX, count);
0516   yGpu.ToDevice(fY, count);
0517   zGpu.ToDevice(fZ, count);
0518
0519   DevicePtr<cuda::SOA3D<T>> gpu_ptr;
0520   gpu_ptr.Allocate();
0521   gpu_ptr.Construct(xGpu, yGpu, zGpu, fSize);
0522 }
0523
0524 #endif // VECGEOM_CUDA_INTERFACE
0525 } // namespace VECGEOM_IMPL_NAMESPACE
0526 } // namespace vecgeom
0527
0528 #endif // VECGEOM_BASE_SOA3D_H_