Use acc::maybeConjugate for cpu code

include/atrip/Operations.hpp

@@ -24,20 +24,13 @@ namespace acc {
 
   // cuda kernels
 
-  template <typename F>
-  __MAYBE_GLOBAL__
-  void zeroing(F* a, size_t n) {
-    F zero = {0};
-    for (size_t i = 0; i < n; i++) {
-      a[i] = zero;
-    }
-  }
-
   ////
   template <typename F>
   __MAYBE_DEVICE__ __MAYBE_HOST__ __INLINE__
   F maybeConjugateScalar(const F &a) { return a; }
 
+  // TODO: instantiate for std::complex<double>
+
 #if defined(HAVE_CUDA)
   template <>
   __MAYBE_DEVICE__ __MAYBE_HOST__ __INLINE__

include/atrip/SliceUnion.hpp

@@ -200,7 +200,7 @@ template <typename F=double>
                            : Slice<F>::Fetch
                            ;
           if (blank.info.state == Slice<F>::SelfSufficient) {
-#if defined(HAVE_CUDA)
+#if defined(HAVE_CUDA) && !defined(ATRIP_SOURCES_IN_GPU)
             const size_t _size = sizeof(F) * sliceSize;
             // TODO: this is code duplication with downstairs
             if (freePointers.size() == 0) {
@@ -221,7 +221,6 @@ template <typename F=double>
                                                          (void*)SOURCES_DATA(sources[from.source]),
                                                          sizeof(F) * sliceSize));
                         ))
-
 #else
             blank.data = SOURCES_DATA(sources[from.source]);
 #endif
@@ -388,6 +387,22 @@ template <typename F=double>
       }
     }
 
+    static size_t
+    getSize(const std::vector<size_t> sliceLength,
+            const std::vector<size_t> paramLength,
+            const size_t np,
+            const MPI_Comm global_world) {
+        const RankMap<F> rankMap(paramLength, np, global_world);
+        const size_t
+          nSources = rankMap.nSources(),
+          sliceSize = std::accumulate(sliceLength.begin(),
+                                      sliceLength.end(),
+                                      1UL,
+                                      std::multiplies<size_t>());
+        return nSources * sliceSize;
+      }
+
+
     // CONSTRUCTOR
     SliceUnion( std::vector<typename Slice<F>::Type> sliceTypes_
               , std::vector<size_t> sliceLength_
@@ -574,12 +589,15 @@ template <typename F=double>
         // TODO: do it through the slice class
         slice.info.state = Slice<F>::Dispatched;
 #if defined(HAVE_CUDA) && defined(ATRIP_SOURCES_IN_GPU)
-#  if !defined(ATRIP_CUDA_AWARE_MPI) 
+#  if !defined(ATRIP_CUDA_AWARE_MPI)
 #    error "You need CUDA aware MPI to have slices on the GPU"
 #  endif
         MPI_Irecv((void*)slice.data,
+#elif defined(HAVE_CUDA) && !defined(ATRIP_SOURCES_IN_GPU)
+        slice.mpi_data = (F*)malloc(sizeof(F) * slice.size);
+        MPI_Irecv(slice.mpi_data,
 #else
-        MPI_Irecv(slice.data,
+        MPI_Irecv((void*)slice.data,
 #endif
                   slice.size,
                   traits::mpi::datatypeOf<F>(),

src/atrip/Atrip.cxx

@@ -160,9 +160,9 @@ Atrip::Output Atrip::run(Atrip::Input<F> const& in) {
 
 
   LOG(0,"Atrip") << "ooo blocks: "
-		 << Atrip::kernelDimensions.ooo.blocks << "\n";
+    << Atrip::kernelDimensions.ooo.blocks << "\n";
   LOG(0,"Atrip") << "ooo threads per block: "
-		 << Atrip::kernelDimensions.ooo.threads << "\n";
+    << Atrip::kernelDimensions.ooo.threads << "\n";
 #endif
 
   // allocate the three scratches, see piecuch
@@ -235,11 +235,54 @@ Atrip::Output Atrip::run(Atrip::Input<F> const& in) {
     MPI_Comm_size(child_comm, &child_size);
   }
 
+  // a, b, c, d, e, f and P => Nv
+  // H                      => No
+  // total_source_sizes contains a list of the number of elements
+  // in all sources of every tensor union, therefore nSlices * sliceSize
+  const std::vector<size_t> total_source_sizes = {
+    // ABPH
+    SliceUnion<F>::getSize({Nv, No}, {Nv, Nv}, (size_t)np, universe),
+    // ABHH
+    SliceUnion<F>::getSize({No, No}, {Nv, Nv}, (size_t)np, universe),
+    // TABHH
+    SliceUnion<F>::getSize({No, No}, {Nv, Nv}, (size_t)np, universe),
+    // TAPHH
+    SliceUnion<F>::getSize({Nv, No, No}, {Nv}, (size_t)np, universe),
+    // HHHA
+    SliceUnion<F>::getSize({No, No, No}, {Nv}, (size_t)np, universe),
+  };
+
+  const size_t
+    total_source_size = sizeof(DataFieldType<F>)
+                      * std::accumulate(total_source_sizes.begin(),
+                                        total_source_sizes.end(),
+                                        0UL);
+
+#if defined(HAVE_CUDA)
+    DataPtr<F> all_sources_pointer;
+    cuMemAlloc(&all_sources_pointer, total_source_size);
+#else
+    DataPtr<F>
+      all_sources_pointer = (DataPtr<F>)malloc(total_source_size);
+#endif
+  size_t _source_pointer_idx = 0;
+
   // BUILD SLICES PARAMETRIZED BY NV x NV =============================={{{1
   WITH_CHRONO("nv-nv-slices",
     LOG(0,"Atrip") << "building NV x NV slices\n";
+    // TODO
+    // DataPtr<F> offseted_pointer = all_sources_pointer
+    //                             * total_source_sizes[_source_pointer_idx++];
     ABPH<F> abph(*in.Vppph, (size_t)No, (size_t)Nv, (size_t)np, child_comm, universe);
+
+    // TODO
+    // DataPtr<F> offseted_pointer = all_sources_pointer
+    //                             * total_source_sizes[_source_pointer_idx++];
     ABHH<F> abhh(*in.Vpphh, (size_t)No, (size_t)Nv, (size_t)np, child_comm, universe);
+
+    // TODO
+    // DataPtr<F> offseted_pointer = all_sources_pointer
+    //                             * total_source_sizes[_source_pointer_idx++];
     TABHH<F> tabhh(*in.Tpphh, (size_t)No, (size_t)Nv, (size_t)np, child_comm, universe);
   )
 
@@ -251,7 +294,13 @@ Atrip::Output Atrip::run(Atrip::Input<F> const& in) {
   // BUILD SLICES PARAMETRIZED BY NV ==================================={{{1
   WITH_CHRONO("nv-slices",
     LOG(0,"Atrip") << "building NV slices\n";
+    // TODO
+    // DataPtr<F> offseted_pointer = all_sources_pointer
+    //                             * total_source_sizes[_source_pointer_idx++];
     TAPHH<F> taphh(*in.Tpphh, (size_t)No, (size_t)Nv, (size_t)np, child_comm, universe);
+    // TODO
+    // DataPtr<F> offseted_pointer = all_sources_pointer
+    //                             * total_source_sizes[_source_pointer_idx++];
     HHHA<F>  hhha(*in.Vhhhp, (size_t)No, (size_t)Nv, (size_t)np, child_comm, universe);
   )
 
@@ -903,5 +952,5 @@ Atrip::Output Atrip::run(Atrip::Input<F> const& in) {
 }
 // instantiate
 template Atrip::Output Atrip::run(Atrip::Input<double> const& in);
-template Atrip::Output Atrip::run(Atrip::Input<Complex> const& in);
+// template Atrip::Output Atrip::run(Atrip::Input<Complex> const& in);
 // Main:1 ends here

src/atrip/Complex.cxx

@@ -21,11 +21,6 @@ namespace atrip {
   template <> double maybeConjugate(const double a) { return a; }
   template <> Complex maybeConjugate(const Complex a) { return std::conj(a); }
 
-#if defined(HAVE_CUDA)
-
-#endif
-
-
   namespace traits {
     template <typename F> bool isComplex() { return false; }
     template <> bool isComplex<double>() { return false; }

src/atrip/Equations.cxx

@@ -13,6 +13,8 @@
 // limitations under the License.
 
 // [[file:~/cuda/atrip/atrip.org::*Prolog][Prolog:2]]
+#include <cstring>
+
 #include<atrip/Equations.hpp>
 
 #include<atrip/CUDA.hpp>
@@ -25,11 +27,8 @@ namespace atrip {
 
 #if defined(HAVE_CUDA)
 #define FOR_K()                                             \
-  for (size_t kmin = blockIdx.x * blockDim.x + threadIdx.x, \
-         k = kmin,                                          \
-         idx = kmin * size * size * size;                   \
-       k < (kmin < size) ? kmin + 1 : size;                 \
-       k++)
+  const size_t k = blockIdx.x * blockDim.x + threadIdx.x; \
+  size_t idx = k*size*size;
 #else
 #define FOR_K() for (size_t k=0, idx=0; k < size; k++)
 #endif
@@ -102,6 +101,7 @@ namespace atrip {
 #  define MIN(a, b) std::min((a), (b))
 #endif
 
+#if defined(ATRIP_NEW_ENERGY)
 
 // [[file:~/cuda/atrip/atrip.org::*Energy][Energy:2]]
 template <typename F>
@@ -250,6 +250,131 @@ void getEnergySame
 }
 // Energy:2 ends here
 
+#else
+
+// [[file:~/cuda/atrip/atrip.org::*Energy][Energy:2]]
+template <typename F>
+__MAYBE_GLOBAL__
+void getEnergyDistinct
+  ( F const epsabc
+  , size_t const No
+  , F* const epsi
+  , F* const Tijk
+  , F* const Zijk
+  , double* _energy
+  ) {
+  constexpr size_t blockSize=16;
+  F energy(0.);
+  for (size_t kk=0; kk<No; kk+=blockSize){
+    const size_t kend( MIN(No, kk+blockSize) );
+    for (size_t jj(kk); jj<No; jj+=blockSize){
+      const size_t jend( MIN( No, jj+blockSize) );
+      for (size_t ii(jj); ii<No; ii+=blockSize){
+        const size_t iend( MIN( No, ii+blockSize) );
+        for (size_t k(kk); k < kend; k++){
+          const F ek(epsi[k]);
+          const size_t jstart = jj > k ? jj : k;
+          for (size_t j(jstart); j < jend; j++){
+            F const ej(epsi[j]);
+            F const facjk = j == k ? F(0.5) : F(1.0);
+            size_t istart = ii > j ? ii : j;
+            for (size_t i(istart); i < iend; i++){
+              const F
+                  ei(epsi[i])
+                , facij = i == j ? F(0.5) : F(1.0)
+                , denominator(epsabc - ei - ej - ek)
+                , U(Zijk[i + No*j + No*No*k])
+                , V(Zijk[i + No*k + No*No*j])
+                , W(Zijk[j + No*i + No*No*k])
+                , X(Zijk[j + No*k + No*No*i])
+                , Y(Zijk[k + No*i + No*No*j])
+                , Z(Zijk[k + No*j + No*No*i])
+                , A(acc::maybeConjugateScalar<F>(Tijk[i + No*j + No*No*k]))
+                , B(acc::maybeConjugateScalar<F>(Tijk[i + No*k + No*No*j]))
+                , C(acc::maybeConjugateScalar<F>(Tijk[j + No*i + No*No*k]))
+                , D(acc::maybeConjugateScalar<F>(Tijk[j + No*k + No*No*i]))
+                , E(acc::maybeConjugateScalar<F>(Tijk[k + No*i + No*No*j]))
+                , _F(acc::maybeConjugateScalar<F>(Tijk[k + No*j + No*No*i]))
+                , value
+                  = 3.0 * ( A * U
+                            + B * V
+                            + C * W
+                            + D * X
+                            + E * Y
+                            + _F * Z )
+                 + ( ( U + X + Y )
+                   - 2.0 * ( V + W + Z )
+                   ) * ( A + D + E )
+                 + ( ( V + W + Z )
+                   - 2.0 * ( U + X + Y )
+                   ) * ( B + C + _F )
+                ;
+              energy += 2.0 * value / denominator * facjk * facij;
+            } // i
+          } // j
+        } // k
+      } // ii
+    } // jj
+  } // kk
+  *_energy = acc::real(energy);
+}
+
+
+template <typename F>
+__MAYBE_GLOBAL__
+void getEnergySame
+  ( F const epsabc
+  , size_t const No
+  , F* const epsi
+  , F* const Tijk
+  , F* const Zijk
+  , double* _energy
+  ) {
+  constexpr size_t blockSize = 16;
+  F energy = F(0.);
+  for (size_t kk=0; kk<No; kk+=blockSize){
+    const size_t kend( MIN( kk+blockSize, No) );
+    for (size_t jj(kk); jj<No; jj+=blockSize){
+      const size_t jend( MIN( jj+blockSize, No) );
+      for (size_t ii(jj); ii<No; ii+=blockSize){
+        const size_t iend( MIN( ii+blockSize, No) );
+        for (size_t k(kk); k < kend; k++){
+          const F ek(epsi[k]);
+          const size_t jstart = jj > k ? jj : k;
+          for(size_t j(jstart); j < jend; j++){
+            const F facjk( j == k ? F(0.5) : F(1.0));
+            const F ej(epsi[j]);
+            const size_t istart = ii > j ? ii : j;
+            for(size_t i(istart); i < iend; i++){
+              const F
+                ei(epsi[i])
+              , facij ( i==j ? F(0.5) : F(1.0))
+              , denominator(epsabc - ei - ej - ek)
+              , U(Zijk[i + No*j + No*No*k])
+              , V(Zijk[j + No*k + No*No*i])
+              , W(Zijk[k + No*i + No*No*j])
+              , A(acc::maybeConjugateScalar<F>(Tijk[i + No*j + No*No*k]))
+              , B(acc::maybeConjugateScalar<F>(Tijk[j + No*k + No*No*i]))
+              , C(acc::maybeConjugateScalar<F>(Tijk[k + No*i + No*No*j]))
+              , value
+                = F(3.0) * ( A * U
+                           + B * V
+                           + C * W
+                           )
+                - ( A + B + C ) * ( U + V + W )
+              ;
+              energy += F(2.0) * value / denominator * facjk * facij;
+            } // i
+          } // j
+        } // k
+      } // ii
+    } // jj
+  } // kk
+  *_energy = acc::real(energy);
+}
+// Energy:2 ends here
+#endif /* defined(ATRIP_NEW_ENERGY) */
+
 // [[file:~/cuda/atrip/atrip.org::*Energy][Energy:3]]
 // instantiate double
 template
@@ -274,6 +399,8 @@ void getEnergySame
   , DataFieldType<double>* energy
   );
 
+// TODO: put this back in
+#if defined(ATRIP_WITH_COMPLEX)
 // instantiate Complex
 template
 __MAYBE_GLOBAL__
@@ -297,6 +424,7 @@ void getEnergySame
   , DataFieldType<double>* energy
   );
 // Energy:3 ends here
+#endif
 
 // [[file:~/cuda/atrip/atrip.org::*Singles%20contribution][Singles contribution:2]]
   template <typename F> __MAYBE_GLOBAL__
@@ -416,7 +544,7 @@ void getEnergySame
 #if defined(ATRIP_USE_DGEMM)
 #if defined(HAVE_CUDA)
 #define REORDER(__II, __JJ, __KK)               \
-  reorder<<<bs, ths>>>(reorder_proxy<           \
+  reorder<<<1, No>>>(reorder_proxy<           \
                        DataFieldType<F>,        \
                        __II ## __JJ ## __KK     \
                        >{},                     \
@@ -454,13 +582,8 @@ void getEnergySame
                   )
 #define MAYBE_CONJ(_conj, _buffer)                                \
   do {                                                            \
-    acc::maybeConjugate<<<                                        \
-                                                                  \
-                            Atrip::kernelDimensions.ooo.blocks,   \
-                                                                  \
-                              Atrip::kernelDimensions.ooo.threads \
-                                                                  \
-                        >>>((DataFieldType<F>*)_conj,             \
+    acc::maybeConjugate<<<1, 1                                    \
+                       >>>((DataFieldType<F>*)_conj,              \
                             (DataFieldType<F>*)_buffer,           \
                             NoNoNo);                              \
   } while (0)
@@ -511,61 +634,39 @@ void getEnergySame
                   _t_buffer,                    \
                   (int const*)&NoNo             \
                   )
-#define MAYBE_CONJ(_conj, _buffer)              \
-  do {                                          \
-    for (size_t __i = 0; __i < NoNoNo; ++__i) { \
-      _conj[__i]                                \
-        = maybeConjugate<F>(_buffer[__i]);      \
-    }                                           \
-  } while (0)
+#define MAYBE_CONJ(_conj, _buffer)               \
+  acc::maybeConjugate((DataFieldType<F>*)_conj,  \
+                      (DataFieldType<F>*)_buffer,\
+                      NoNoNo);
 #endif
 
     F one{1.0}, m_one{-1.0}, zero{0.0};
     const size_t NoNoNo = No*NoNo;
+
+// Zeroing vectors
 #ifdef HAVE_CUDA
-//    DataFieldType<F>* _t_buffer;
-//    DataFieldType<F>* _vhhh;
-//    WITH_CHRONO("double:cuda:alloc",
-//    _CHECK_CUDA_SUCCESS("Allocating _t_buffer",
-//                        cuMemAlloc((CUdeviceptr*)&_t_buffer,
-//                                   NoNoNo * sizeof(DataFieldType<F>)));
-//    _CHECK_CUDA_SUCCESS("Allocating _vhhh",
-//                        cuMemAlloc((CUdeviceptr*)&_vhhh,
-//                                   NoNoNo * sizeof(DataFieldType<F>)));
-//                )
+
 #if !defined(ATRIP_ONLY_DGEMM)
-    // we still have to zero this
-    const size_t
-      bs = Atrip::kernelDimensions.ooo.blocks,
-      ths = Atrip::kernelDimensions.ooo.threads;
-    acc::zeroing<<<bs, ths>>>((DataFieldType<F>*)_t_buffer, NoNoNo);
-    acc::zeroing<<<bs, ths>>>((DataFieldType<F>*)_vhhh, NoNoNo);
+    {
+    const size_t elements = NoNoNo * sizeof(DataFieldType<F>)/4;
+    WITH_CHRONO("double:zeroing",
+                _CHECK_CUDA_SUCCESS("Zeroing Tijk",
+                  cuMemsetD32_v2((CUdeviceptr)Tijk, 0x00, elements));
+                _CHECK_CUDA_SUCCESS("Zeroing t buffer",
+                  cuMemsetD32_v2((CUdeviceptr)_t_buffer, 0x00, elements));
+                _CHECK_CUDA_SUCCESS("Zeroing vhhh buffer",
+                  cuMemsetD32_v2((CUdeviceptr)_vhhh, 0x00, elements));
+                )
+    }
 #endif
 
 #else
     DataFieldType<F>* _t_buffer = (DataFieldType<F>*)malloc(NoNoNo * sizeof(F));
     DataFieldType<F>* _vhhh = (DataFieldType<F>*)malloc(NoNoNo * sizeof(F));
-    DataFieldType<F> zero_h{0.0};
-    for (size_t i=0; i < NoNoNo; i++) {
-      _t_buffer[i] = zero_h;
-      _vhhh[i] = zero_h;
-    }
-#endif
-
-    // Set Tijk to zero
-#if defined(HAVE_CUDA) && !defined(ATRIP_ONLY_DGEMM)
-    WITH_CHRONO("double:reorder",
-                acc::zeroing<<<bs, ths>>>((DataFieldType<F>*)Tijk,
-                                           NoNoNo);
-                )
-#endif
-
-#if !defined(HAVE_CUDA)
-    WITH_CHRONO("double:reorder",
-      for (size_t k = 0; k < NoNoNo; k++) {
-        Tijk[k] = DataFieldType<F>{0.0};
-       })
-#endif /* !defined(HAVE_CUDA) */
+    std::memset((void*)_t_buffer, 0x00, NoNoNo * sizeof(DataFieldType<F>));
+    std::memset((void*)_vhhh,     0x00, NoNoNo * sizeof(DataFieldType<F>));
+    std::memset((void*)Tijk,      0x00, NoNoNo * sizeof(DataFieldType<F>));
+#endif /* HAVE_CUDA */
 
 
 #if defined(ATRIP_ONLY_DGEMM)