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   SOA3D(SOA3D<T> const &other);
0053
0054   SOA3D() = default;
0055
0056   VECCORE_ATT_HOST_DEVICE
0057   SOA3D &operator=(SOA3D<T> const &other);
0058
0059   VECCORE_ATT_HOST_DEVICE
0060   ~SOA3D();
0061
0062   VECCORE_ATT_HOST_DEVICE
0063   VECGEOM_FORCE_INLINE
0064   size_t size() const;
0065
0066   VECCORE_ATT_HOST_DEVICE
0067   VECGEOM_FORCE_INLINE
0068   size_t capacity() const;
0069
0070   VECCORE_ATT_HOST_DEVICE
0071   VECGEOM_FORCE_INLINE
0072   void resize(size_t newSize);
0073
0074   VECCORE_ATT_HOST_DEVICE
0075   VECGEOM_FORCE_INLINE
0076   void reserve(size_t newCapacity);
0077
0078   VECGEOM_FORCE_INLINE
0079   void clear();
0080
0081   // Element access methods
0082
0083   VECCORE_ATT_HOST_DEVICE
0084   VECGEOM_FORCE_INLINE
0085   Vector3D<T> operator[](size_t index) const;
0086
0087   VECCORE_ATT_HOST_DEVICE
0088   VECGEOM_FORCE_INLINE
0089   T x(size_t index) const;
0090
0091   VECCORE_ATT_HOST_DEVICE
0092   VECGEOM_FORCE_INLINE
0093   T &x(size_t index);
0094
0095   VECCORE_ATT_HOST_DEVICE
0096   VECGEOM_FORCE_INLINE
0097   T *x();
0098
0099   VECCORE_ATT_HOST_DEVICE
0100   VECGEOM_FORCE_INLINE
0101   T const *x() const;
0102
0103   VECCORE_ATT_HOST_DEVICE
0104   VECGEOM_FORCE_INLINE
0105   T y(size_t index) const;
0106
0107   VECCORE_ATT_HOST_DEVICE
0108   VECGEOM_FORCE_INLINE
0109   T &y(size_t index);
0110
0111   VECCORE_ATT_HOST_DEVICE
0112   VECGEOM_FORCE_INLINE
0113   T *y();
0114
0115   VECCORE_ATT_HOST_DEVICE
0116   VECGEOM_FORCE_INLINE
0117   T const *y() const;
0118
0119   VECCORE_ATT_HOST_DEVICE
0120   VECGEOM_FORCE_INLINE
0121   T z(size_t index) const;
0122
0123   VECCORE_ATT_HOST_DEVICE
0124   VECGEOM_FORCE_INLINE
0125   T &z(size_t index);
0126
0127   VECCORE_ATT_HOST_DEVICE
0128   VECGEOM_FORCE_INLINE
0129   T *z();
0130
0131   VECCORE_ATT_HOST_DEVICE
0132   VECGEOM_FORCE_INLINE
0133   T const *z() const;
0134
0135   // Element manipulation methods
0136
0137   VECCORE_ATT_HOST_DEVICE
0138   VECGEOM_FORCE_INLINE
0139   void set(size_t index, T x, T y, T z);
0140
0141   VECCORE_ATT_HOST_DEVICE
0142   VECGEOM_FORCE_INLINE
0143   void set(size_t index, Vector3D<T> const &vec);
0144
0145   VECCORE_ATT_HOST_DEVICE
0146   VECGEOM_FORCE_INLINE
0147   void push_back(T x, T y, T z);
0148
0149   VECCORE_ATT_HOST_DEVICE
0150   VECGEOM_FORCE_INLINE
0151   void push_back(Vector3D<T> const &vec);
0152
0153 #ifdef VECGEOM_CUDA_INTERFACE
0154   DevicePtr<cuda::SOA3D<T>> CopyToGpu(DevicePtr<T> xGpu, DevicePtr<T> yGpu, DevicePtr<T> zGpu) const;
0155   DevicePtr<cuda::SOA3D<T>> CopyToGpu(DevicePtr<T> xGpu, DevicePtr<T> yGpu, DevicePtr<T> zGpu, size_t size) const;
0156 #endif // VECGEOM_CUDA_INTERFACE
0157
0158   void ToFile(std::string /*filename*/) const;
0159   int FromFile(std::string /*filename*/);
0160
0161 private:
0162   VECCORE_ATT_HOST_DEVICE
0163   void Allocate();
0164
0165   VECCORE_ATT_HOST_DEVICE
0166   void Deallocate();
0167 };
0168
0169 #if defined(GCC_DIAG_POP_NEEDED)
0170
0171 #pragma GCC diagnostic pop
0172 #undef GCC_DIAG_POP_NEEDED
0173
0174 #endif
0175
0176 template <typename T>
0177 VECCORE_ATT_HOST_DEVICE
0178 SOA3D<T>::SOA3D(T *xval, T *yval, T *zval, size_t sz)
0179     : fAllocated(false), fSize(sz), fCapacity(fSize), fX(xval), fY(yval), fZ(zval)
0180 {
0181 }
0182
0183 template <typename T>
0184 VECCORE_ATT_HOST_DEVICE
0185 SOA3D<T>::SOA3D(size_t sz) : fSize(sz), fCapacity(sz)
0186 {
0187   Allocate();
0188 }
0189
0190 template <typename T>
0191 SOA3D<T>::SOA3D(SOA3D<T> const &rhs)
0192     : fAllocated(false), fSize(rhs.fSize), fCapacity(rhs.fCapacity)
0193 {
0194   if (rhs.fAllocated) {
0195     Allocate();
0196     copy(rhs.fX, rhs.fX + rhs.fSize, fX);
0197     copy(rhs.fY, rhs.fY + rhs.fSize, fY);
0198     copy(rhs.fZ, rhs.fZ + rhs.fSize, fZ);
0199   } else {
0200     fX = rhs.fX;
0201     fY = rhs.fY;
0202     fZ = rhs.fZ;
0203   }
0204 }
0205
0206 template <typename T>
0207 VECCORE_ATT_HOST_DEVICE
0208 SOA3D<T> &SOA3D<T>::operator=(SOA3D<T> const &rhs)
0209 {
0210 #ifndef VECCORE_CUDA_DEVICE_COMPILATION
0211   fSize     = rhs.fSize;
0212   fCapacity = rhs.fCapacity;
0213   Deallocate();
0214   if (rhs.fAllocated && rhs.fSize > 0) {
0215     Allocate();
0216     copy(rhs.fX, rhs.fX + rhs.fSize, fX);
0217     copy(rhs.fY, rhs.fY + rhs.fSize, fY);
0218     copy(rhs.fZ, rhs.fZ + rhs.fSize, fZ);
0219   } else {
0220     fX = rhs.fX;
0221     fY = rhs.fY;
0222     fZ = rhs.fZ;
0223   }
0224 #else
0225   fAllocated = false;
0226   fSize      = rhs.fSize;
0227   fCapacity  = rhs.fCapacity;
0228   fX         = rhs.fX;
0229   fY         = rhs.fY;
0230   fZ         = rhs.fZ;
0231 #endif
0232   return *this;
0233 }
0234
0235 template <typename T>
0236 VECCORE_ATT_HOST_DEVICE
0237 SOA3D<T>::~SOA3D()
0238 {
0239 #ifndef VECCORE_CUDA_DEVICE_COMPILATION
0240   Deallocate();
0241 #endif
0242 }
0243
0244 template <typename T>
0245 VECCORE_ATT_HOST_DEVICE
0246 size_t SOA3D<T>::size() const
0247 {
0248   return fSize;
0249 }
0250
0251 template <typename T>
0252 VECCORE_ATT_HOST_DEVICE
0253 size_t SOA3D<T>::capacity() const
0254 {
0255   return fCapacity;
0256 }
0257
0258 template <typename T>
0259 VECCORE_ATT_HOST_DEVICE
0260 void SOA3D<T>::resize(size_t newSize)
0261 {
0262   assert(newSize <= fCapacity);
0263   fSize = newSize;
0264 }
0265
0266 template <typename T>
0267 VECCORE_ATT_HOST_DEVICE
0268 void SOA3D<T>::reserve(size_t newCapacity)
0269 {
0270   fCapacity = newCapacity;
0271   T *xNew, *yNew, *zNew;
0272   xNew  = AlignedAllocate<T>(fCapacity);
0273   yNew  = AlignedAllocate<T>(fCapacity);
0274   zNew  = AlignedAllocate<T>(fCapacity);
0275   fSize = (fSize > fCapacity) ? fCapacity : fSize;
0276   if (fX && fY && fZ) {
0277     copy(fX, fX + fSize, xNew);
0278     copy(fY, fY + fSize, yNew);
0279     copy(fZ, fZ + fSize, zNew);
0280   }
0281   Deallocate();
0282   fX         = xNew;
0283   fY         = yNew;
0284   fZ         = zNew;
0285   fAllocated = true;
0286 }
0287
0288 template <typename T>
0289 void SOA3D<T>::clear()
0290 {
0291   Deallocate();
0292   fSize     = 0;
0293   fCapacity = 0;
0294 }
0295
0296 template <typename T>
0297 VECCORE_ATT_HOST_DEVICE
0298 void SOA3D<T>::Allocate()
0299 {
0300   if (fCapacity == 0) return;
0301
0302   fX         = AlignedAllocate<T>(fCapacity);
0303   fY         = AlignedAllocate<T>(fCapacity);
0304   fZ         = AlignedAllocate<T>(fCapacity);
0305   fAllocated = true;
0306 }
0307
0308 template <typename T>
0309 VECCORE_ATT_HOST_DEVICE
0310 void SOA3D<T>::Deallocate()
0311 {
0312   if (fAllocated) {
0313     AlignedFree(fX);
0314     AlignedFree(fY);
0315     AlignedFree(fZ);
0316   }
0317   fAllocated = false;
0318 }
0319
0320 template <typename T>
0321 VECCORE_ATT_HOST_DEVICE
0322 Vector3D<T> SOA3D<T>::operator[](size_t index) const
0323 {
0324   return Vector3D<T>(fX[index], fY[index], fZ[index]);
0325 }
0326
0327 template <typename T>
0328 VECCORE_ATT_HOST_DEVICE
0329 T SOA3D<T>::x(size_t index) const
0330 {
0331   return fX[index];
0332 }
0333
0334 template <typename T>
0335 VECCORE_ATT_HOST_DEVICE
0336 T &SOA3D<T>::x(size_t index)
0337 {
0338   return fX[index];
0339 }
0340
0341 template <typename T>
0342 VECCORE_ATT_HOST_DEVICE
0343 T *SOA3D<T>::x()
0344 {
0345   return fX;
0346 }
0347
0348 template <typename T>
0349 VECCORE_ATT_HOST_DEVICE
0350 T const *SOA3D<T>::x() const
0351 {
0352   return fX;
0353 }
0354
0355 template <typename T>
0356 VECCORE_ATT_HOST_DEVICE
0357 T SOA3D<T>::y(size_t index) const
0358 {
0359   return fY[index];
0360 }
0361
0362 template <typename T>
0363 VECCORE_ATT_HOST_DEVICE
0364 T &SOA3D<T>::y(size_t index)
0365 {
0366   return fY[index];
0367 }
0368
0369 template <typename T>
0370 VECCORE_ATT_HOST_DEVICE
0371 T *SOA3D<T>::y()
0372 {
0373   return fY;
0374 }
0375
0376 template <typename T>
0377 VECCORE_ATT_HOST_DEVICE
0378 T const *SOA3D<T>::y() const
0379 {
0380   return fY;
0381 }
0382
0383 template <typename T>
0384 VECCORE_ATT_HOST_DEVICE
0385 T SOA3D<T>::z(size_t index) const
0386 {
0387   return fZ[index];
0388 }
0389
0390 template <typename T>
0391 VECCORE_ATT_HOST_DEVICE
0392 T &SOA3D<T>::z(size_t index)
0393 {
0394   return fZ[index];
0395 }
0396
0397 template <typename T>
0398 VECCORE_ATT_HOST_DEVICE
0399 T *SOA3D<T>::z()
0400 {
0401   return fZ;
0402 }
0403
0404 template <typename T>
0405 VECCORE_ATT_HOST_DEVICE
0406 T const *SOA3D<T>::z() const
0407 {
0408   return fZ;
0409 }
0410
0411 template <typename T>
0412 VECGEOM_FORCE_INLINE
0413 VECCORE_ATT_HOST_DEVICE
0414 void SOA3D<T>::set(size_t index, T xval, T yval, T zval)
0415 {
0416 // not asserting in case of NVCC -- still getting annoying
0417 // errors on CUDA < 8.0
0418 #ifndef VECCORE_CUDA
0419   assert(index < fCapacity);
0420 #endif
0421   fX[index] = xval;
0422   fY[index] = yval;
0423   fZ[index] = zval;
0424 }
0425
0426 template <typename T>
0427 VECGEOM_FORCE_INLINE
0428 VECCORE_ATT_HOST_DEVICE
0429 void SOA3D<T>::set(size_t index, Vector3D<T> const &vec)
0430 {
0431 // not asserting in case of NVCC -- still getting annoying
0432 // errors on CUDA < 8.0
0433 #ifndef VECCORE_CUDA
0434   assert(index < fCapacity);
0435 #endif
0436   fX[index] = vec[0];
0437   fY[index] = vec[1];
0438   fZ[index] = vec[2];
0439 }
0440
0441 template <typename T>
0442 VECCORE_ATT_HOST_DEVICE
0443 void SOA3D<T>::push_back(T xval, T yval, T zval)
0444 {
0445   fX[fSize] = xval;
0446   fY[fSize] = yval;
0447   fZ[fSize] = zval;
0448   ++fSize;
0449 }
0450
0451 template <typename T>
0452 VECCORE_ATT_HOST_DEVICE
0453 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 }
0526 } // End global namespace
0527
0528 #endif // VECGEOM_BASE_SOA3D_H_