Use acc::maybeConjugate for cpu code

include/atrip/CUDA.hpp

@@ -11,11 +11,22 @@
 #if defined(HAVE_CUDA) && defined(__CUDACC__)
 #  define __MAYBE_GLOBAL__ __global__
 #  define __MAYBE_DEVICE__ __device__
+#  define __MAYBE_HOST__ __host__
+#  define __INLINE__ __inline__
 #else
 #  define __MAYBE_GLOBAL__
 #  define __MAYBE_DEVICE__
+#  define __MAYBE_HOST__
+#  define __INLINE__ inline
 #endif
 
+#if defined(HAVE_CUDA)
+#define ACC_FUNCALL(fname, i, j, ...) fname<<<(i), (j)>>>(__VA_ARGS__)
+#else
+#define ACC_FUNCALL(fname, i, j, ...) fname(__VA_ARGS__)
+#endif /*  defined(HAVE_CUDA) */
+
+
 #define _CHECK_CUDA_SUCCESS(message, ...)                               \
   do {                                                                  \
     CUresult result = __VA_ARGS__;                                      \

include/atrip/Equations.hpp

@@ -23,6 +23,8 @@
 #include<thrust/device_vector.h>
 #endif
 
+#include<atrip/CUDA.hpp>
+
 
 namespace atrip {
 using ABCTuple = std::array<size_t, 3>;
@@ -32,21 +34,25 @@ using ABCTuples = std::vector<ABCTuple>;
 
 // [[file:~/cuda/atrip/atrip.org::*Energy][Energy:1]]
 template <typename F=double>
-double getEnergyDistinct
+__MAYBE_GLOBAL__
+void getEnergyDistinct
   ( F const epsabc
   , size_t const No
   , F* const epsi
   , F* const Tijk
   , F* const Zijk
+  , double* energy
   );
 
 template <typename F=double>
-double getEnergySame
+__MAYBE_GLOBAL__
+void getEnergySame
   ( F const epsabc
   , size_t const No
   , F* const epsi
   , F* const Tijk
   , F* const Zijk
+  , double* energy
   );
 // Energy:1 ends here
 
@@ -97,6 +103,11 @@ void singlesContribution
     // -- TIJK
     // , DataPtr<F> Tijk
     , DataFieldType<F>* Tijk_
+#if defined(HAVE_CUDA)
+     // -- tmp buffers
+    , DataFieldType<F>* _t_buffer
+    , DataFieldType<F>* _vhhh 
+#endif
     );
 // Doubles contribution:1 ends here
 

include/atrip/Operations.hpp

@@ -0,0 +1,164 @@
+// Copyright 2022 Alejandro Gallo
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef OPERATIONS_HPP_
+#define OPERATIONS_HPP_
+
+#include <atrip/CUDA.hpp>
+#include <atrip/Types.hpp>
+#include <atrip/Complex.hpp>
+
+namespace atrip {
+namespace acc {
+
+  // cuda kernels
+
+  ////
+  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__
+  cuDoubleComplex maybeConjugateScalar(const cuDoubleComplex &a) {
+    return {a.x, -a.y};
+  }
+#endif /*  defined(HAVE_CUDA) */
+
+  template <typename F>
+  __MAYBE_GLOBAL__
+  void maybeConjugate(F* to, F* from, size_t n) {
+    for (size_t i = 0; i < n; ++i) {
+      to[i] = maybeConjugateScalar<F>(from[i]);
+    }
+  }
+
+
+  template <typename F>
+  __MAYBE_DEVICE__ __MAYBE_HOST__
+  void reorder(F* to, F* from, size_t size, size_t I, size_t J, size_t K) {
+    size_t idx = 0;
+    const size_t IDX = I + J*size + K*size*size;
+    for (size_t k = 0; k < size; k++)
+    for (size_t j = 0; j < size; j++)
+    for (size_t i = 0; i < size; i++, idx++)
+      to[idx] += from[IDX];
+  }
+
+  // Multiplication operation
+  //////////////////////////////////////////////////////////////////////////////
+
+  template <typename F>
+  __MAYBE_DEVICE__ __MAYBE_HOST__ __INLINE__
+  F prod(const F &a, const F &b) { return a * b; }
+
+#if defined(HAVE_CUDA)
+  template <>
+  __MAYBE_DEVICE__ __MAYBE_HOST__ __INLINE__
+  cuDoubleComplex prod(const cuDoubleComplex &a, const cuDoubleComplex &b) {
+    return cuCmul(a, b);
+  }
+#endif /*  defined(HAVE_CUDA) */
+
+  // Division operation
+  //////////////////////////////////////////////////////////////////////////////
+
+  template <typename F>
+  __MAYBE_DEVICE__ __MAYBE_HOST__ __INLINE__
+  F div(const F &a, const F &b) { return a / b; }
+
+#if defined(HAVE_CUDA)
+  template <>
+  __MAYBE_DEVICE__ __MAYBE_HOST__ __INLINE__
+  cuDoubleComplex div(const cuDoubleComplex &a, const cuDoubleComplex &b) {
+    return cuCdiv(a, b);
+  }
+#endif /*  defined(HAVE_CUDA) */
+
+  // Real part
+  //////////////////////////////////////////////////////////////////////////////
+
+  template <typename F>
+  __MAYBE_HOST__ __INLINE__
+  double real(F &a) { return std::real(a); }
+
+  template <>
+  __MAYBE_DEVICE__ __MAYBE_HOST__ __INLINE__
+  double real(double &a) {
+    return a;
+  }
+
+#if defined(HAVE_CUDA)
+  template <>
+  __MAYBE_DEVICE__ __MAYBE_HOST__ __INLINE__
+  double real(cuDoubleComplex &a) {
+    return cuCreal(a);
+  }
+#endif /*  defined(HAVE_CUDA) */
+
+  // Substraction operator
+  //////////////////////////////////////////////////////////////////////////////
+
+  template <typename F>
+  __MAYBE_DEVICE__ __MAYBE_HOST__ __INLINE__
+  F sub(const F &a, const F &b) { return a - b; }
+
+#if defined(HAVE_CUDA)
+  template <>
+  __MAYBE_DEVICE__ __MAYBE_HOST__ __INLINE__
+  cuDoubleComplex sub(const cuDoubleComplex &a,
+                               const cuDoubleComplex &b) {
+    return cuCsub(a, b);
+  }
+#endif /*  defined(HAVE_CUDA) */
+
+  // Addition operator
+  //////////////////////////////////////////////////////////////////////////////
+
+  template <typename F>
+  __MAYBE_DEVICE__ __MAYBE_HOST__ __INLINE__
+  F add(const F &a, const F &b) { return a + b; }
+
+#if defined(HAVE_CUDA)
+  template <>
+  __MAYBE_DEVICE__ __MAYBE_HOST__ __INLINE__
+  cuDoubleComplex add(const cuDoubleComplex &a, const cuDoubleComplex &b) {
+    return cuCadd(a, b);
+  }
+#endif /*  defined(HAVE_CUDA) */
+
+  // Sum in place operator
+  //////////////////////////////////////////////////////////////////////////////
+
+  template <typename F>
+  __MAYBE_DEVICE__ __MAYBE_HOST__
+  void sum_in_place(F* to, const F* from) { *to += *from; }
+
+#if defined(HAVE_CUDA)
+  template <>
+  __MAYBE_DEVICE__ __MAYBE_HOST__
+  void sum_in_place(cuDoubleComplex* to, const cuDoubleComplex* from) {
+    to->x += from->x;
+    to->y += from->y;
+  }
+#endif /*  defined(HAVE_CUDA) */
+
+
+}  // namespace acc
+}  // namespace atrip
+
+#endif

include/atrip/Slice.hpp

@@ -352,7 +352,7 @@ Info info;
 
 // [[file:~/cuda/atrip/atrip.org::*Attributes][Attributes:2]]
 DataPtr<F> data;
-#if defined(HAVE_CUDA)
+#if defined(HAVE_CUDA) && !defined (ATRIP_SOURCES_IN_GPU)
     F* mpi_data;
 #endif
 // Attributes:2 ends here
@@ -456,7 +456,7 @@ void unwrapAndMarkReady() {
       if (errorCode != MPI_SUCCESS)
         throw "Atrip: Unexpected error MPI ERROR";
 
-#if defined(HAVE_CUDA)
+#if defined(HAVE_CUDA) && !defined(ATRIP_SOURCES_IN_GPU)
       // copy the retrieved mpi data to the device
       WITH_CHRONO("cuda:memcpy",
                   _CHECK_CUDA_SUCCESS("copying mpi data to device",
@@ -488,7 +488,7 @@ void unwrapAndMarkReady() {
 Slice(size_t size_)
       : info({})
       , data(DataNullPtr)
-#if defined(HAVE_CUDA)
+#if defined(HAVE_CUDA) && !defined(ATRIP_SOURCES_IN_GPU)
       , mpi_data(nullptr)
 #endif
       , size(size_)

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_
@@ -405,6 +420,7 @@ template <typename F=double>
               , sliceSize(std::accumulate(sliceLength.begin(),
                                           sliceLength.end(),
                                           1UL, std::multiplies<size_t>()))
+
 #if defined(ATRIP_SOURCES_IN_GPU)
               , sources(rankMap.nSources())
 #else
@@ -417,6 +433,23 @@ template <typename F=double>
     { // constructor begin
 
       LOG(0,"Atrip") << "INIT SliceUnion: " << name << "\n";
+        printf("sliceSize %d, number of slices %d\n\n\n", sliceSize, sources.size());
+
+#if defined(ATRIP_SOURCES_IN_GPU)
+      for (auto& ptr: sources) {
+        const CUresult sourceError =
+          cuMemAlloc(&ptr, sizeof(F) * sliceSize);
+        if (ptr == 0UL) {
+          throw "UNSUFICCIENT MEMORY ON THE GRAPHIC CARD FOR SOURCES";
+        }
+        if (sourceError != CUDA_SUCCESS) {
+          std::stringstream s;
+          s << "Error allocating memory for sources "
+            << "code " << sourceError << "\n";
+          throw s.str();
+        }
+      }
+#endif
 
       for (auto& ptr: sliceBuffers) {
 #if defined(HAVE_CUDA)
@@ -445,6 +478,34 @@ template <typename F=double>
                      std::inserter(freePointers, freePointers.begin()),
                      [](DataPtr<F> ptr) { return ptr; });
 
+#if defined(HAVE_CUDA)
+      LOG(1,"Atrip") << "warming communication up " << slices.size() << "\n";
+      WITH_CHRONO("cuda:warmup",
+                  int nRanks=Atrip::np, requestCount=0;
+                  int nSends=sliceBuffers.size()*nRanks;
+                  MPI_Request *requests = (MPI_Request*) malloc(nSends*2 * sizeof(MPI_Request));
+                  MPI_Status *statuses = (MPI_Status*) malloc(nSends*2 * sizeof(MPI_Status));
+                  for (int sliceId=0; sliceId<sliceBuffers.size(); sliceId++){
+                    for (int rankId=0; rankId<nRanks; rankId++){
+                      MPI_Isend((void*)SOURCES_DATA(sources[0]),
+                                sliceSize,
+                                traits::mpi::datatypeOf<F>(),
+                                rankId,
+                                100,
+                                universe,
+                                &requests[requestCount++]);
+                      MPI_Irecv((void*)sliceBuffers[sliceId],
+                                sliceSize,
+                                traits::mpi::datatypeOf<F>(),
+                                rankId,
+                                100,
+                                universe,
+                                &requests[requestCount++]);
+                    }
+                  }
+                  MPI_Waitall(nSends*2, requests, statuses);
+                  )
+#endif
 
 
       LOG(1,"Atrip") << "#slices " << slices.size() << "\n";
@@ -527,14 +588,16 @@ template <typename F=double>
       if (slice.info.state == Slice<F>::Fetch) { // if-1
         // TODO: do it through the slice class
         slice.info.state = Slice<F>::Dispatched;
-#if defined(HAVE_CUDA)
-#  if !defined(ATRIP_CUDA_AWARE_MPI) && defined(ATRIP_SOURCES_IN_GPU)
+#if defined(HAVE_CUDA) && defined(ATRIP_SOURCES_IN_GPU)
+#  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

@@ -202,7 +202,7 @@ Atrip::Output Atrip::run(Atrip::Input<F> const& in) {
   _CHECK_CUDA_SUCCESS("Zijk",
                       cuMemAlloc(&Zijk, sizeof(F) * No * No * No));
 #else
-  std::vector<F> &Tai = _Tai, &epsi = _epsi, &epsa = _epsa;
+  DataPtr<F> Tai = _Tai.data(), epsi = _epsi.data(), epsa = _epsa.data();
   Zijk = (DataFieldType<F>*)malloc(No*No*No * sizeof(DataFieldType<F>));
   Tijk = (DataFieldType<F>*)malloc(No*No*No * sizeof(DataFieldType<F>));
 #endif
@@ -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,13 +294,38 @@ 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);
   )
 
   // all tensors
   std::vector< SliceUnion<F>* > unions = {&taphh, &hhha, &abph, &abhh, &tabhh};
 
+#ifdef HAVE_CUDA
+    // TODO: free buffers
+    DataFieldType<F>* _t_buffer;
+    DataFieldType<F>* _vhhh;
+    WITH_CHRONO("double:cuda:alloc",
+    _CHECK_CUDA_SUCCESS("Allocating _t_buffer",
+                        cuMemAlloc((CUdeviceptr*)&_t_buffer,
+                                   No*No*No * sizeof(DataFieldType<F>)));
+    _CHECK_CUDA_SUCCESS("Allocating _vhhh",
+                        cuMemAlloc((CUdeviceptr*)&_vhhh,
+                                   No*No*No * sizeof(DataFieldType<F>)));
+                )
+    //const size_t
+     // bs = Atrip::kernelDimensions.ooo.blocks,
+      //ths = Atrip::kernelDimensions.ooo.threads;
+    //cuda::zeroing<<<bs, ths>>>((DataFieldType<F>*)_t_buffer, NoNoNo);
+    //cuda::zeroing<<<bs, ths>>>((DataFieldType<F>*)_vhhh, NoNoNo);
+#endif
+
   // get tuples for the current rank
   TuplesDistribution *distribution;
 
@@ -639,7 +707,14 @@ Atrip::Output Atrip::run(Atrip::Input<F> const& in) {
                                          tabhh.unwrapSlice(Slice<F>::AC, abc),
                                          tabhh.unwrapSlice(Slice<F>::BC, abc),
                                          // -- TIJK
-                                         (DataFieldType<F>*)Tijk);
+                                         (DataFieldType<F>*)Tijk
+#if defined(HAVE_CUDA)
+                                         // -- tmp buffers
+                                         ,(DataFieldType<F>*)_t_buffer
+                                         ,(DataFieldType<F>*)_vhhh
+#endif
+                                         );
+
                   WITH_RANK << iteration << "-th doubles done\n";
       ))
     }
@@ -667,7 +742,7 @@ Atrip::Output Atrip::run(Atrip::Input<F> const& in) {
                                       (DataFieldType<F>*)Tai,
 #else
       singlesContribution<F>(No, Nv, abc[0], abc[1], abc[2],
-                             Tai.data(),
+                             Tai,
 #endif
                              (DataFieldType<F>*)abhh.unwrapSlice(Slice<F>::AB,
                                                                  abc),
@@ -683,31 +758,71 @@ Atrip::Output Atrip::run(Atrip::Input<F> const& in) {
     // COMPUTE ENERGY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% {{{1
 #if defined(ATRIP_ONLY_DGEMM)
     if (false)
-#endif
+#endif /* defined(ATRIP_ONLY_DGEMM) */
     if (!isFakeTuple(i)) {
-      double tupleEnergy(0.);
+#if defined(HAVE_CUDA)
+      double *tupleEnergy;
+      cuMemAlloc((DataPtr<double>*)&tupleEnergy, sizeof(double));
+#else
+      double _tupleEnergy(0.);
+      double *tupleEnergy = &_tupleEnergy;
+#endif /* defined(HAVE_CUDA) */
 
       int distinct(0);
       if (abc[0] == abc[1]) distinct++;
       if (abc[1] == abc[2]) distinct--;
-      const F epsabc(_epsa[abc[0]] + _epsa[abc[1]] + _epsa[abc[2]]);
+      const double
+        epsabc = std::real(_epsa[abc[0]] + _epsa[abc[1]] + _epsa[abc[2]]);
+
+      DataFieldType<F> _epsabc{epsabc};
 
-      // LOG(0, "AtripCUDA") <<  "doing energy " << i << "distinct " << distinct << "\n";
       WITH_CHRONO("energy",
-/*
-    TODO: think about how to do this on the GPU in the best way possible
-        if ( distinct == 0)
-          tupleEnergy = getEnergyDistinct<F>(epsabc, No, (F*)epsi, (F*)Tijk, (F*)Zijk);
-        else
-          tupleEnergy = getEnergySame<F>(epsabc, No, (F*)epsi, (F*)Tijk, (F*)Zijk);
-*/
-      )
+                  if ( distinct == 0) {
+                    ACC_FUNCALL(getEnergyDistinct<DataFieldType<F>>,
+                                1, 1, // for cuda
+                                _epsabc,
+                                No,
+#if defined(HAVE_CUDA)
+                                (DataFieldType<F>*)epsi,
+                                (DataFieldType<F>*)Tijk,
+                                (DataFieldType<F>*)Zijk,
+#else
+                                epsi,
+                                Tijk,
+                                Zijk,
+#endif
+                                tupleEnergy);
+                  } else {
+                    ACC_FUNCALL(getEnergySame<DataFieldType<F>>,
+                                1, 1, // for cuda
+                                _epsabc,
+                                No,
+#if defined(HAVE_CUDA)
+                                (DataFieldType<F>*)epsi,
+                                (DataFieldType<F>*)Tijk,
+                                (DataFieldType<F>*)Zijk,
+#else
+                                epsi,
+                                Tijk,
+                                Zijk,
+#endif
+                                tupleEnergy);
+                  })
+
+#if defined(HAVE_CUDA)
+        double host_tuple_energy;
+        cuMemcpyDtoH((void*)&host_tuple_energy,
+                     (DataPtr<double>)tupleEnergy,
+                     sizeof(double));
+#else
+        double host_tuple_energy = *tupleEnergy;
+#endif /* defined(HAVE_CUDA) */
 
 #if defined(HAVE_OCD) || defined(ATRIP_PRINT_TUPLES)
-      tupleEnergies[abc] = tupleEnergy;
+      tupleEnergies[abc] = host_tuple_energy;
 #endif
 
-      energy += tupleEnergy;
+      energy += host_tuple_energy;
 
     }
 
@@ -837,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,106 +13,22 @@
 // limitations under the License.
 
 // [[file:~/cuda/atrip/atrip.org::*Prolog][Prolog:2]]
+#include <cstring>
+
 #include<atrip/Equations.hpp>
 
 #include<atrip/CUDA.hpp>
+#include<atrip/Operations.hpp>
 
 namespace atrip {
 // Prolog:2 ends here
 
 
 
-#ifdef HAVE_CUDA
-namespace cuda {
-
-  // cuda kernels
-
-  template <typename F>
-  __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>
-  __device__
-  F maybeConjugateScalar(const F a);
-
-  template <>
-  __device__
-  double maybeConjugateScalar(const double a) { return a; }
-
-  template <>
-  __device__
-  cuDoubleComplex
-  maybeConjugateScalar(const cuDoubleComplex a) {
-    return {a.x, -a.y};
-  }
-
-  template <typename F>
-  __global__
-  void maybeConjugate(F* to, F* from, size_t n) {
-    for (size_t i = 0; i < n; ++i) {
-      to[i] = maybeConjugateScalar<F>(from[i]);
-    }
-  }
-
-
-  template <typename F>
-  __global__
-  void reorder(F* to, F* from, size_t size, size_t I, size_t J, size_t K) {
-    size_t idx = 0;
-    const size_t IDX = I + J*size + K*size*size;
-    for (size_t k = 0; k < size; k++)
-    for (size_t j = 0; j < size; j++)
-    for (size_t i = 0; i < size; i++, idx++)
-      to[idx] += from[IDX];
-  }
-
-  // I mean, really CUDA... really!?
-  template <typename F>
-  __device__
-  F multiply(const F &a, const F &b);
-  template <>
-  __device__
-  double multiply(const double &a, const double &b) { return a * b; }
-
-  template <>
-  __device__
-  cuDoubleComplex multiply(const cuDoubleComplex &a, const cuDoubleComplex &b) {
-    return
-      {a.x * b.x - a.y * b.y,
-       a.x * b.y + a.y * b.x};
-  }
-
-  template <typename F>
-  __device__
-  void sum_in_place(F* to, const F* from);
-
-  template <>
-  __device__
-  void sum_in_place(double* to, const double *from) { *to += *from; }
-
-  template <>
-  __device__
-  void sum_in_place(cuDoubleComplex* to, const cuDoubleComplex* from) {
-    to->x += from->x;
-    to->y += from->y;
-  }
-
-};
-#endif
-
 #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
@@ -133,7 +49,7 @@ namespace cuda {
     _REORDER_BODY_(__VA_ARGS__)                 \
       }
 #if defined(HAVE_CUDA)
-#define GO(__TO, __FROM) cuda::sum_in_place<F>(&__TO, &__FROM);
+#define GO(__TO, __FROM) acc::sum_in_place<F>(&__TO, &__FROM);
 #else
 #define GO(__TO, __FROM) __TO += __FROM;
 #endif
@@ -179,24 +95,182 @@ namespace cuda {
 #undef _IJK_
 #undef GO
 
+#if defined(HAVE_CUDA)
+#  define MIN(a, b) min((a), (b))
+#else
+#  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>
-double getEnergyDistinct
+__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}
+                , eijk(acc::add(acc::add(ei, ej), ek))
+                , denominator(acc::sub(epsabc, eijk))
+                , 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(Tijk[i + No*j + No*No*k]))
+                , B(acc::maybeConjugateScalar(Tijk[i + No*k + No*No*j]))
+                , C(acc::maybeConjugateScalar(Tijk[j + No*i + No*No*k]))
+                , D(acc::maybeConjugateScalar(Tijk[j + No*k + No*No*i]))
+                , E(acc::maybeConjugateScalar(Tijk[k + No*i + No*No*j]))
+                , _F(acc::maybeConjugateScalar(Tijk[k + No*j + No*No*i]))
+                , AU = acc::prod(A, U)
+                , BV = acc::prod(B, V)
+                , CW = acc::prod(C, W)
+                , DX = acc::prod(D, X)
+                , EY = acc::prod(E, Y)
+                , FZ = acc::prod(_F, Z)
+                , UXY = acc::add(U, acc::add(X, Y))
+                , VWZ = acc::add(V, acc::add(W, Z))
+                , ADE = acc::add(A, acc::add(D, E))
+                , BCF = acc::add(B, acc::add(C, _F))
+                // I just might as well write this in CL
+                , _first =  acc::add(AU,
+                                     acc::add(BV,
+                                              acc::add(CW,
+                                                       acc::add(DX,
+                                                                acc::add(EY, FZ)))))
+                , _second = acc::prod(acc::sub(UXY,
+                                               acc::prod(F{-2.0}, VWZ)),
+                                      ADE)
+                , _third = acc::prod(acc::sub(VWZ,
+                                              acc::prod(F{-2.0}, UXY)),
+                                     BCF)
+                , value = acc::add(acc::prod(F{3.0}, _first),
+                                   acc::add(_second,
+                                            _third))
+                , _loop_energy = acc::prod(acc::prod(F{2.0}, value),
+                                           acc::div(acc::prod(facjk, facij),
+                                                    denominator))
+                ;
+              acc::sum_in_place(&_energy, &_loop_energy);
+            } // i
+          } // j
+        } // k
+      } // ii
+    } // jj
+  } // kk
+  const double real_part = acc::real(_energy);
+  acc::sum_in_place(energy, &real_part);
+}
+
+
+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})
+              , eijk(acc::add(acc::add(ei, ej), ek))
+              , denominator(acc::sub(epsabc, eijk))
+              , 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(Tijk[i + No*j + No*No*k]))
+              , B(acc::maybeConjugateScalar(Tijk[j + No*k + No*No*i]))
+              , C(acc::maybeConjugateScalar(Tijk[k + No*i + No*No*j]))
+              , ABC = acc::add(A, acc::add(B, C))
+              , UVW = acc::add(U, acc::add(V, W))
+              , AU = acc::prod(A, U)
+              , BV = acc::prod(B, V)
+              , CW = acc::prod(C, W)
+              , AU_and_BV_and_CW = acc::add(acc::add(AU, BV), CW)
+              , value = acc::sub(acc::prod(F{3.0}, AU_and_BV_and_CW),
+                                 acc::prod(ABC, UVW))
+              , _loop_energy = acc::prod(acc::prod(F{2.0}, value),
+                                         acc::div(acc::prod(facjk, facij),
+                                                  denominator))
+              ;
+
+              acc::sum_in_place(&_energy, &_loop_energy);
+            } // i
+          } // j
+        } // k
+      } // ii
+    } // jj
+  } // kk
+  const double real_part = acc::real(_energy);
+  acc::sum_in_place(energy, &real_part);
+}
+// 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( std::min(No, kk+blockSize) );
+    const size_t kend( MIN(No, kk+blockSize) );
     for (size_t jj(kk); jj<No; jj+=blockSize){
-      const size_t jend( std::min( No, jj+blockSize) );
+      const size_t jend( MIN( No, jj+blockSize) );
       for (size_t ii(jj); ii<No; ii+=blockSize){
-        const size_t iend( std::min( 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;
@@ -215,12 +289,12 @@ double getEnergyDistinct
                 , 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(maybeConjugate<F>(Tijk[i + No*j + No*No*k]))
-                , B(maybeConjugate<F>(Tijk[i + No*k + No*No*j]))
-                , C(maybeConjugate<F>(Tijk[j + No*i + No*No*k]))
-                , D(maybeConjugate<F>(Tijk[j + No*k + No*No*i]))
-                , E(maybeConjugate<F>(Tijk[k + No*i + No*No*j]))
-                , _F(maybeConjugate<F>(Tijk[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
@@ -242,26 +316,28 @@ double getEnergyDistinct
       } // ii
     } // jj
   } // kk
-  return std::real(energy);
+  *_energy = acc::real(energy);
 }
 
 
 template <typename F>
-double getEnergySame
+__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( std::min( kk+blockSize, No) );
+    const size_t kend( MIN( kk+blockSize, No) );
     for (size_t jj(kk); jj<No; jj+=blockSize){
-      const size_t jend( std::min( jj+blockSize, No) );
+      const size_t jend( MIN( jj+blockSize, No) );
       for (size_t ii(jj); ii<No; ii+=blockSize){
-        const size_t iend( std::min( ii+blockSize, No) );
+        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;
@@ -277,9 +353,9 @@ double getEnergySame
               , 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(maybeConjugate<F>(Tijk[i + No*j + No*No*k]))
-              , B(maybeConjugate<F>(Tijk[j + No*k + No*No*i]))
-              , C(maybeConjugate<F>(Tijk[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
@@ -294,49 +370,61 @@ double getEnergySame
       } // ii
     } // jj
   } // kk
-  return std::real(energy);
+  *_energy = acc::real(energy);
 }
 // Energy:2 ends here
+#endif /* defined(ATRIP_NEW_ENERGY) */
 
 // [[file:~/cuda/atrip/atrip.org::*Energy][Energy:3]]
 // instantiate double
 template
-double getEnergyDistinct
-  ( double const epsabc
+__MAYBE_GLOBAL__
+void getEnergyDistinct
+  ( DataFieldType<double> const epsabc
   , size_t const No
-  , double* const epsi
-  , double* const Tijk
-  , double* const Zijk
+  , DataFieldType<double>* const epsi
+  , DataFieldType<double>* const Tijk
+  , DataFieldType<double>* const Zijk
+  , DataFieldType<double>* energy
   );
 
 template
-double getEnergySame
-  ( double const epsabc
+__MAYBE_GLOBAL__
+void getEnergySame
+  ( DataFieldType<double> const epsabc
   , size_t const No
-  , double* const epsi
-  , double* const Tijk
-  , double* const Zijk
+  , DataFieldType<double>* const epsi
+  , DataFieldType<double>* const Tijk
+  , DataFieldType<double>* const Zijk
+  , DataFieldType<double>* energy
   );
 
+// TODO: put this back in
+#if defined(ATRIP_WITH_COMPLEX)
 // instantiate Complex
 template
-double getEnergyDistinct
-  ( Complex const epsabc
+__MAYBE_GLOBAL__
+void getEnergyDistinct
+  ( DataFieldType<Complex> const epsabc
   , size_t const No
-  , Complex* const epsi
-  , Complex* const Tijk
-  , Complex* const Zijk
+  , DataFieldType<Complex>* const epsi
+  , DataFieldType<Complex>* const Tijk
+  , DataFieldType<Complex>* const Zijk
+  , DataFieldType<double>* energy
   );
 
 template
-double getEnergySame
-  ( Complex const epsabc
+__MAYBE_GLOBAL__
+void getEnergySame
+  ( DataFieldType<Complex> const epsabc
   , size_t const No
-  , Complex* const epsi
-  , Complex* const Tijk
-  , Complex* const Zijk
+  , DataFieldType<Complex>* const epsi
+  , DataFieldType<Complex>* const Tijk
+  , DataFieldType<Complex>* const Zijk
+  , DataFieldType<double>* energy
   );
 // Energy:3 ends here
+#endif
 
 // [[file:~/cuda/atrip/atrip.org::*Singles%20contribution][Singles contribution:2]]
   template <typename F> __MAYBE_GLOBAL__
@@ -360,18 +448,26 @@ double getEnergySame
       const size_t ijk = i + j*No + k*NoNo;
 
 #ifdef HAVE_CUDA
-#  define GO(__TPH, __VABIJ)                                          \
-    {                                                                 \
-      const DataFieldType<F> product                                  \
-            = cuda::multiply<DataFieldType<F>>((__TPH), (__VABIJ));   \
-      cuda::sum_in_place<DataFieldType<F>>(&Zijk[ijk], &product);     \
-    }
+
+#define GO(__TPH, __VABIJ)                              \
+  do {                                                  \
+    const DataFieldType<F>                              \
+      product = acc::prod<DataFieldType<F>>((__TPH),    \
+                                            (__VABIJ)); \
+    acc::sum_in_place<DataFieldType<F>>(&Zijk[ijk],     \
+                                        &product);      \
+  } while (0)
+
 #else
-#  define GO(__TPH, __VABIJ) Zijk[ijk] += (__TPH) * (__VABIJ);
+
+#define GO(__TPH, __VABIJ) Zijk[ijk] += (__TPH) * (__VABIJ)
+
 #endif
-      GO(Tph[ a + i * Nv ], VBCij[ j + k * No ])
-      GO(Tph[ b + j * Nv ], VACij[ i + k * No ])
-      GO(Tph[ c + k * Nv ], VABij[ i + j * No ])
+
+      GO(Tph[ a + i * Nv ], VBCij[ j + k * No ]);
+      GO(Tph[ b + j * Nv ], VACij[ i + k * No ]);
+      GO(Tph[ c + k * Nv ], VABij[ i + j * No ]);
+
 #undef GO
     } // for loop j
   }
@@ -433,16 +529,22 @@ double getEnergySame
     // -- TIJK
     // , DataPtr<F> Tijk_
     , DataFieldType<F>* Tijk_
+#if defined(HAVE_CUDA)
+     // -- tmp buffers
+    , DataFieldType<F>* _t_buffer
+    , DataFieldType<F>* _vhhh 
+#endif
      ) {
-
-    const size_t NoNo = No*No;
+    const size_t a = abc[0], b = abc[1], c = abc[2]
+              , NoNo = No*No
+              ;
 
     DataFieldType<F>* Tijk = (DataFieldType<F>*)Tijk_;
 
 #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     \
                        >{},                     \
@@ -480,12 +582,7 @@ double getEnergySame
                   )
 #define MAYBE_CONJ(_conj, _buffer)                                \
   do {                                                            \
-    cuda::maybeConjugate<<<                                       \
-                                                                  \
-                            Atrip::kernelDimensions.ooo.blocks,   \
-                                                                  \
-                              Atrip::kernelDimensions.ooo.threads \
-                                                                  \
+    acc::maybeConjugate<<<1, 1                                    \
                        >>>((DataFieldType<F>*)_conj,              \
                             (DataFieldType<F>*)_buffer,           \
                             NoNoNo);                              \
@@ -538,58 +635,38 @@ double getEnergySame
                   (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)
+  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>)));
-                )
-    const size_t
-      bs = Atrip::kernelDimensions.ooo.blocks,
-      ths = Atrip::kernelDimensions.ooo.threads;
 
 #if !defined(ATRIP_ONLY_DGEMM)
-    cuda::zeroing<<<bs, ths>>>((DataFieldType<F>*)_t_buffer, NoNoNo);
-    cuda::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",
-                cuda::zeroing<<<bs, ths>>>((DataFieldType<F>*)Tijk,
-                                           NoNoNo);
-                )
-#else
-    WITH_CHRONO("double:reorder",
-      for (size_t k = 0; k < NoNoNo; k++) {
-        Tijk[k] = DataFieldType<F>{0.0};
-       })
-#endif
+    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)
@@ -597,7 +674,7 @@ double getEnergySame
 #undef REORDER
 #define MAYBE_CONJ(a, b) do {} while(0)
 #define REORDER(i, j, k) do {} while(0)
-#endif
+#endif /* defined(ATRIP_ONLY_DGEMM) */
 
     // HOLES
     WITH_CHRONO("doubles:holes",
@@ -681,16 +758,16 @@ double getEnergySame
 #ifdef HAVE_CUDA
     // we need to synchronize here since we need
     // the Tijk for next process in the pipeline
-    _CHECK_CUDA_SUCCESS("Synchronizing",
-                        cuCtxSynchronize());
-    _CHECK_CUDA_SUCCESS("Freeing _vhhh",
-                        cuMemFree((CUdeviceptr)_vhhh));
-    _CHECK_CUDA_SUCCESS("Freeing _t_buffer",
-                        cuMemFree((CUdeviceptr)_t_buffer));
+    //_CHECK_CUDA_SUCCESS("Synchronizing",
+    //                    cuCtxSynchronize());
+    //_CHECK_CUDA_SUCCESS("Freeing _vhhh",
+    //                    cuMemFree((CUdeviceptr)_vhhh));
+    //_CHECK_CUDA_SUCCESS("Freeing _t_buffer",
+    //                    cuMemFree((CUdeviceptr)_t_buffer));
 #else
     free(_vhhh);
     free(_t_buffer);
-#endif
+#endif  /* defined(HAVE_CUDA) */
   }
 
   #undef REORDER
@@ -741,7 +818,7 @@ double getEnergySame
       }
 
     }
-#endif
+#endif /* defined(ATRIP_USE_DGEMM) */
   }
 
 
@@ -773,6 +850,12 @@ double getEnergySame
     , DataPtr<double> const TBChh
     // -- TIJK
     , DataFieldType<double>* Tijk
+#if defined(HAVE_CUDA)
+     // -- tmp buffers
+    , DataFieldType<double>* _t_buffer
+    , DataFieldType<double>* _vhhh 
+#endif
+
     );
 
   template
@@ -801,6 +884,12 @@ double getEnergySame
     , DataPtr<Complex> const TBChh
     // -- TIJK
     , DataFieldType<Complex>* Tijk
+#if defined(HAVE_CUDA)
+     // -- tmp buffers
+    , DataFieldType<Complex>* _t_buffer
+    , DataFieldType<Complex>* _vhhh 
+#endif
+
     );
 // Doubles contribution:2 ends here