Add preliminary openacc support, atrip bench not linking

.github/workflows/main.yml

@@ -18,6 +18,8 @@ jobs:
     strategy:
       matrix:
         compiler:
+          - gcc12
+          - gcc11
           - gcc11
           - gcc10
           - gcc9

bench/Makefile.am

@@ -23,6 +23,18 @@ atrip_SOURCES = main.cxx
 atrip_CPPFLAGS = $(AM_CPPFLAGS)
 atrip_LDADD = $(BENCHES_LDADD)
 
+atrip: main.cxx
+	$(NVCXX) -cuda     \
+		-x cu -I../      \
+		$(MPILIBS)       \
+		-I$(srcdir)/     \
+		$(AM_CPPFLAGS)   \
+		$(DEFS)          \
+		$(BENCHES_LDADD) \
+		$(AM_LDFLAGS)    \
+		$< -o $@
+
+endif
 
 if !WITH_CUDA
 ##

configure.ac

@@ -48,7 +48,8 @@ AM_CONDITIONAL([WITH_CLANG_CHECK], [test x${clang_check} = xYES])
 AC_ARG_ENABLE([cuda],
               [AS_HELP_STRING([--enable-cuda],
                               [Build with cuda])],
-                              [WITH_CUDA=yes],
+              [WITH_CUDA=yes
+               WITH_OPENACC=yes],
               [WITH_CUDA=no])
 AC_ARG_VAR([NVCC], [Path to the nvidia cuda compiler.])
 AC_ARG_VAR([CUDA_LDFLAGS], [LDFLAGS to find libraries -lcuda, -lcudart, -lcublas.])
@@ -182,6 +183,8 @@ if test x${WITH_CUDA} = xyes; then
                      -----------------------
   ])
   AC_CHECK_PROGS([NVCC], [nvcc])
+  AC_CHECK_PROGS([NVCXX], [nvc++])
+  MPILIBS=$($MPICXX -show | awk '!($1="")')
   AC_SUBST([CUDA_LDFLAGS])
   AC_DEFINE([HAVE_CUDA],1,[Wether we are using CUDA])
   # TODO: make sure to find cuda and cudart
@@ -227,6 +230,7 @@ AC_MSG_RESULT([
                        ATRIP_LDFLAGS  = $ATRIP_LDFLAGS
                        BLAS = ${BLAS_LIBS}
                        LIBS = ${LIBS}
+                       MPILIBS = $MPILIBS
 ])
 
 AC_OUTPUT

etc/env/raven/cuda-openacc

@@ -0,0 +1,82 @@
+#!/usr/bin/env bash
+mods=(
+  #cuda/11.6
+  nvhpcsdk/22  # for openacc
+  gcc/12
+  openmpi
+  mkl/2020.4
+  autoconf/2.69
+  automake/1.15
+  libtool/2.4.6
+)
+
+
+module purge
+module load ${mods[@]}
+
+LIB_PATH="${NVHPC_CUDA_HOME}/lib64"
+export CUBLAS_LD_PATH="${NVHPC_ROOT}/math_libs/lib64/"
+export CUDA_ROOT=${CUDA_HOME}
+export CUDA_LDFLAGS="-L${LIB_PATH} -lcuda -L${LIB_PATH} -lcudart -L${CUBLAS_LD_PATH} -lcublas"
+export CUDA_CXXFLAGS="-I${CUDA_HOME}/include"
+
+export LD_LIBRARY_PATH="${MKL_HOME}/lib/intel64:${LD_LIBRARY_PATH}"
+
+MPILIBS=$(mpicxx -show | awk '!($1="")')
+export MPILIBS
+export MPINVCXX="nv++ ${MPILIBS}"
+
+ls ${CUBLAS_LD_PATH}/libcublas.so
+ls ${LIB_PATH}/libcudart.so
+
+#export OMPI_CC="nvc"
+#export OMPI_CXX="nvc++"
+
+BLAS_LDFLAGS="-L${PWD}/OpenBLAS-0.3.20/ -lopenblas"
+_openblas_make () {
+
+  [[ -d OpenBLAS-0.3.20/ ]] || {
+    wget https://github.com/xianyi/OpenBLAS/releases/download/v0.3.20/OpenBLAS-0.3.20.tar.gz
+    tar xvzf OpenBLAS-0.3.20.tar.gz
+    cd OpenBLAS-0.3.20/
+    make FC=gfortran CC=gcc USE_OPENMP=1 NUM_THREADS=72 TARGET=SKYLAKEX
+  } && {
+    echo "Openblas built"
+  }
+
+}
+( _openblas_make; )
+
+
+
+cat <<EOF
+
+////////////////////////////////////////////////////////////////////////////////
+                                    info
+////////////////////////////////////////////////////////////////////////////////
+
+
+MKL_HOME         = $MKL_HOME
+BLAS_STATIC_PATH = $BLAS_STATIC_PATH
+
+CUDA_ROOT      = ${CUDA_HOME}
+CUDA_LDFLAGS   = "-L${LIB_PATH} -lcuda -L${LIB_PATH} -lcudart -L${LIB_PATH} -lcublas"
+CUDA_CXXFLAGS  = "-I${CUDA_HOME}/include"
+
+
+
+Consider now runnng the following
+
+../../configure \\
+  --enable-cuda \\
+  --disable-slice \\
+  --with-blas="${BLAS_LDFLAGS}" \\
+  CXX="gcc" \\
+  NVCC="\$MPINVCXX" \\
+  MPICXX="mpicxx"
+
+
+EOF
+	
+
+return

etc/m4/atrip_openacc.m4

@@ -0,0 +1,70 @@
+#   SYNOPSIS
+#
+#   ATRIP_OPENACC([ACTION-SUCCESS], [ACTION-FAILURE])
+#
+#   DESCRIPTION
+#
+#   Check whether the given the -fopenacc flag works with the current language's compiler
+#   or gives an error.
+#
+#   ACTION-SUCCESS/ACTION-FAILURE are shell commands to execute on
+#   success/failure.
+#
+#   LICENSE
+#
+#   Copyright (c) 2023 Alejandro Gallo <aamsgallo@gmail.com>
+#
+#   Copying and distribution of this file, with or without modification, are
+#   permitted in any medium without royalty provided the copyright notice
+#   and this notice are preserved.  This file is offered as-is, without any
+#   warranty.
+
+AC_DEFUN([ATRIP_OPENACC],
+[
+AC_MSG_CHECKING([that the compiler works with the -fopenacc])
+AC_COMPILE_IFELSE([AC_LANG_SOURCE([_ATRIP_OPENACC_SOURCE])],
+                  [
+                   $1
+                   AC_MSG_RESULT([yes])
+                   ],
+                  [
+                   $2
+                   AC_MSG_ERROR([no])
+                   ])
+])dnl DEFUN
+
+m4_define([_ATRIP_OPENACC_SOURCE], [[
+#include <stdio.h>
+#include <stdlib.h>
+#include <openacc.h>
+
+#define SIZE 10
+
+int main(int argc, char **argv) {
+  float matrix[SIZE * SIZE];
+  float result[SIZE * SIZE];
+
+  // Initialize the matrix with random values
+  for (int i = 0; i < SIZE * SIZE; i++) {
+    matrix[i] = rand() / (float)RAND_MAX;
+  }
+
+#pragma acc data                                            \
+        copy(matrix[0:SIZE * SIZE])           \
+        copyout(result[0:SIZE * SIZE])
+  {
+    // Calculate the matrix multiplication
+#pragma acc parallel loop collapse(2)
+    for (int i = 0; i < SIZE; i++) {
+      for (int j = 0; j < SIZE; j++) {
+        float sum = 0.0f;
+        for (int k = 0; k < SIZE; k++) {
+          sum += matrix[i * SIZE + k] * matrix[j * SIZE + k];
+        }
+        result[i * SIZE + j] = sum;
+      }
+    }
+  }
+  return 0;
+}
+]])

etc/m4/ax_cxx_compile_stdcxx.m4

@@ -43,7 +43,7 @@
 #   and this notice are preserved.  This file is offered as-is, without any
 #   warranty.
 
-#serial 14
+#serial 15
 
 dnl  This macro is based on the code from the AX_CXX_COMPILE_STDCXX_11 macro
 dnl  (serial version number 13).
@@ -189,7 +189,11 @@ m4_define([_AX_CXX_COMPILE_STDCXX_testbody_new_in_11], [[
 
 #error "This is not a C++ compiler"
 
-#elif __cplusplus < 201103L
+// MSVC always sets __cplusplus to 199711L in older versions; newer versions
+// only set it correctly if /Zc:__cplusplus is specified as well as a
+// /std:c++NN switch:
+// https://devblogs.microsoft.com/cppblog/msvc-now-correctly-reports-__cplusplus/
+#elif __cplusplus < 201103L && !defined _MSC_VER
 
 #error "This is not a C++11 compiler"
 
@@ -480,7 +484,7 @@ m4_define([_AX_CXX_COMPILE_STDCXX_testbody_new_in_14], [[
 
 #error "This is not a C++ compiler"
 
-#elif __cplusplus < 201402L
+#elif __cplusplus < 201402L && !defined _MSC_VER
 
 #error "This is not a C++14 compiler"
 
@@ -604,7 +608,7 @@ m4_define([_AX_CXX_COMPILE_STDCXX_testbody_new_in_17], [[
 
 #error "This is not a C++ compiler"
 
-#elif __cplusplus < 201703L
+#elif __cplusplus < 201703L && !defined _MSC_VER
 
 #error "This is not a C++17 compiler"
 
@@ -970,7 +974,7 @@ namespace cxx17
 
 }  // namespace cxx17
 
-#endif  // __cplusplus < 201703L
+#endif  // __cplusplus < 201703L && !defined _MSC_VER
 
 ]])
 
@@ -983,7 +987,7 @@ m4_define([_AX_CXX_COMPILE_STDCXX_testbody_new_in_20], [[
 
 #error "This is not a C++ compiler"
 
-#elif __cplusplus < 202002L
+#elif __cplusplus < 202002L && !defined _MSC_VER
 
 #error "This is not a C++20 compiler"
 
@@ -1000,6 +1004,6 @@ namespace cxx20
 
 }  // namespace cxx20
 
-#endif  // __cplusplus < 202002L
+#endif  // __cplusplus < 202002L && !defined _MSC_VER
 
 ]])

include/atrip/CUDA.hpp

@@ -11,22 +11,11 @@
 #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,8 +23,6 @@
 #include<thrust/device_vector.h>
 #endif
 
-#include<atrip/CUDA.hpp>
-
 
 namespace atrip {
 using ABCTuple = std::array<size_t, 3>;
@@ -34,25 +32,21 @@ using ABCTuples = std::vector<ABCTuple>;
 
 // [[file:~/cuda/atrip/atrip.org::*Energy][Energy:1]]
 template <typename F=double>
-__MAYBE_GLOBAL__
-void getEnergyDistinct
+double getEnergyDistinct
   ( F const epsabc
   , size_t const No
   , F* const epsi
   , F* const Tijk
   , F* const Zijk
-  , double* energy
   );
 
 template <typename F=double>
-__MAYBE_GLOBAL__
-void getEnergySame
+double getEnergySame
   ( F const epsabc
   , size_t const No
   , F* const epsi
   , F* const Tijk
   , F* const Zijk
-  , double* energy
   );
 // Energy:1 ends here
 
@@ -103,11 +97,6 @@ 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

@@ -1,171 +0,0 @@
-// 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_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; }
-
-#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) && !defined (ATRIP_SOURCES_IN_GPU)
+#if defined(HAVE_CUDA)
     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) && !defined(ATRIP_SOURCES_IN_GPU)
+#if defined(HAVE_CUDA)
       // 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) && !defined(ATRIP_SOURCES_IN_GPU)
+#if defined(HAVE_CUDA)
       , mpi_data(nullptr)
 #endif
       , size(size_)

include/atrip/SliceUnion.hpp

@@ -405,7 +405,6 @@ 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
@@ -418,23 +417,6 @@ 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)
@@ -463,34 +445,6 @@ 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";
@@ -573,11 +527,12 @@ 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) && defined(ATRIP_SOURCES_IN_GPU)
-#  if !defined(ATRIP_CUDA_AWARE_MPI) 
+#if defined(HAVE_CUDA)
+#  if !defined(ATRIP_CUDA_AWARE_MPI) && defined(ATRIP_SOURCES_IN_GPU)
 #    error "You need CUDA aware MPI to have slices on the GPU"
 #  endif
-        MPI_Irecv((void*)slice.data,
+        slice.mpi_data = (F*)malloc(sizeof(F) * slice.size);
+        MPI_Irecv(slice.mpi_data,
 #else
         MPI_Irecv(slice.data,
 #endif

src/Makefile.am

@@ -7,16 +7,17 @@ AM_CPPFLAGS = $(CTF_CPPFLAGS)
 lib_LIBRARIES = libatrip.a
 
 libatrip_a_CPPFLAGS = -I$(top_srcdir)/include/
-libatrip_a_SOURCES = ./atrip/Blas.cxx ./atrip/Tuples.cxx ./atrip/DatabaseCommunicator.cxx
+libatrip_a_SOURCES =
 NVCC_FILES  = ./atrip/Equations.cxx ./atrip/Complex.cxx ./atrip/Atrip.cxx
-
+NVCC_FILES += ./atrip/Blas.cxx ./atrip/Tuples.cxx ./atrip/DatabaseCommunicator.cxx
 if WITH_CUDA
 NVCC_OBJS = $(patsubst %.cxx,%.nvcc.o,$(NVCC_FILES))
 libatrip_a_CPPFLAGS += $(CUDA_CXXFLAGS)
 libatrip_a_DEPENDENCIES = $(NVCC_OBJS)
 libatrip_a_LIBADD = $(NVCC_OBJS)
 %.nvcc.o: %.cxx
-	$(NVCC) -c -x cu -ccbin="${MPICXX}" -I../ $(CPPFLAGS) $(CTF_CPPFLAGS) $(DEFS) $(libatrip_a_CPPFLAGS) $< -o $@
+	##$(NVCC) -c -x cu -ccbin="${MPICXX}" -I../ $(CPPFLAGS) $(CTF_CPPFLAGS) $(DEFS) $(libatrip_a_CPPFLAGS) $< -o $@
+	$(NVCXX) -cuda $(MPILIBS) -c -x cu -I../ $(CPPFLAGS) $(CTF_CPPFLAGS) $(DEFS) $(libatrip_a_CPPFLAGS) $< -o $@
 
 #./atrip/Equations.o: ./atrip/Equations.cxx
 #	$(NVCC) -c -I../ $(CPPFLAGS) $(libatrip_a_CPPFLAGS) $< -o $@

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
-  DataPtr<F> Tai = _Tai.data(), epsi = _epsi.data(), epsa = _epsa.data();
+  std::vector<F> &Tai = _Tai, &epsi = _epsi, &epsa = _epsa;
   Zijk = (DataFieldType<F>*)malloc(No*No*No * sizeof(DataFieldType<F>));
   Tijk = (DataFieldType<F>*)malloc(No*No*No * sizeof(DataFieldType<F>));
 #endif
@@ -258,25 +258,6 @@ Atrip::Output Atrip::run(Atrip::Input<F> const& in) {
   // 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;
 
@@ -658,14 +639,7 @@ 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
-#if defined(HAVE_CUDA)
-                                         // -- tmp buffers
-                                         ,(DataFieldType<F>*)_t_buffer
-                                         ,(DataFieldType<F>*)_vhhh
-#endif
-                                         );
-
+                                         (DataFieldType<F>*)Tijk);
                   WITH_RANK << iteration << "-th doubles done\n";
       ))
     }
@@ -693,7 +667,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,
+                             Tai.data(),
 #endif
                              (DataFieldType<F>*)abhh.unwrapSlice(Slice<F>::AB,
                                                                  abc),
@@ -709,71 +683,28 @@ Atrip::Output Atrip::run(Atrip::Input<F> const& in) {
     // COMPUTE ENERGY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% {{{1
 #if defined(ATRIP_ONLY_DGEMM)
     if (false)
-#endif /* defined(ATRIP_ONLY_DGEMM) */
+#endif
     if (!isFakeTuple(i)) {
-#if defined(HAVE_CUDA)
-      double *tupleEnergy;
-      cuMemAlloc((DataPtr<double>*)&tupleEnergy, sizeof(double));
-#else
-      double _tupleEnergy(0.);
-      double *tupleEnergy = &_tupleEnergy;
-#endif /* defined(HAVE_CUDA) */
+      double tupleEnergy(0.);
 
       int distinct(0);
       if (abc[0] == abc[1]) distinct++;
       if (abc[1] == abc[2]) distinct--;
-      const double
-        epsabc = std::real(_epsa[abc[0]] + _epsa[abc[1]] + _epsa[abc[2]]);
-
-      DataFieldType<F> _epsabc{epsabc};
+      const F epsabc(_epsa[abc[0]] + _epsa[abc[1]] + _epsa[abc[2]]);
 
+      // LOG(0, "AtripCUDA") <<  "doing energy " << i << "distinct " << distinct << "\n";
       WITH_CHRONO("energy",
-                  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 ( 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 defined(HAVE_OCD) || defined(ATRIP_PRINT_TUPLES)
-      tupleEnergies[abc] = host_tuple_energy;
+      tupleEnergies[abc] = tupleEnergy;
 #endif
 
-      energy += host_tuple_energy;
+      energy += tupleEnergy;
 
     }
 

src/atrip/Equations.cxx

@@ -16,13 +16,96 @@
 #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, \
@@ -50,7 +133,7 @@ namespace atrip {
     _REORDER_BODY_(__VA_ARGS__)                 \
       }
 #if defined(HAVE_CUDA)
-#define GO(__TO, __FROM) acc::sum_in_place<F>(&__TO, &__FROM);
+#define GO(__TO, __FROM) cuda::sum_in_place<F>(&__TO, &__FROM);
 #else
 #define GO(__TO, __FROM) __TO += __FROM;
 #endif
@@ -96,205 +179,181 @@ namespace atrip {
 #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
-
 
 // [[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
-  ) {
+  __MAYBE_DEVICE__
+double getEnergyDistinct
+  (F const epsabc,
+   size_t const No,
+   F* const epsi,
+   F* const Tijk,
+   F* const Zijk) {
   constexpr size_t blockSize=16;
-  F _energy = {0.};
+  F energy(0.);
+#if defined(HAVE_CUDA)
+  #pragma acc kernels
+  for (size_t k(0); k < No; k++) {
+    for (size_t j(k); j < No; j++) {
+      for (size_t i(j); i < No; i++) {
+#else
   for (size_t kk=0; kk<No; kk+=blockSize){
-    const size_t kend( MIN(No, kk+blockSize) );
+    const size_t kend( std::min(No, kk+blockSize) );
     for (size_t jj(kk); jj<No; jj+=blockSize){
-      const size_t jend( MIN( No, jj+blockSize) );
+      const size_t jend( std::min( No, jj+blockSize) );
       for (size_t ii(jj); ii<No; ii+=blockSize){
-        const size_t iend( MIN( No, ii+blockSize) );
+        const size_t iend( std::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++){
+#endif
+              const F ek(epsi[k]);
+              const F ej(epsi[j]);
+              const F facjk = j == k ? F(0.5) : F(1.0);
               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))
+                , 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(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))
+                , 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]))
+                , 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 )
                 ;
-              acc::sum_in_place(&_energy, &_loop_energy);
+              energy += 2.0 * value / denominator * facjk * facij;
             } // i
           } // j
         } // k
+#if !defined(HAVE_CUDA)
       } // ii
     } // jj
   } // kk
-  const double real_part = acc::real(_energy);
-  acc::sum_in_place(energy, &real_part);
+#endif
+  return std::real(energy);
 }
 
 
 template <typename F>
-__MAYBE_GLOBAL__
-void getEnergySame
+  __MAYBE_DEVICE__
+double 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.};
+  F energy = F(0.);
+#if defined(HAVE_CUDA)
+  #pragma acc kernels
+  for (size_t k(0); k < No; k++) {
+    for (size_t j(k); j < No; j++) {
+      for (size_t i(j); i < No; i++) {
+#else
   for (size_t kk=0; kk<No; kk+=blockSize){
-    const size_t kend( MIN( kk+blockSize, No) );
+    const size_t kend( std::min( kk+blockSize, No) );
     for (size_t jj(kk); jj<No; jj+=blockSize){
-      const size_t jend( MIN( jj+blockSize, No) );
+      const size_t jend( std::min( jj+blockSize, No) );
       for (size_t ii(jj); ii<No; ii+=blockSize){
-        const size_t iend( MIN( ii+blockSize, No) );
+        const size_t iend( std::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++){
+#endif
+              const F facjk( j == k ? F(0.5) : F(1.0));
+              const F ek(epsi[k]);
+              const F ej(epsi[j]);
               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))
+              , 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(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))
+              , 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]))
+              , value
+                = F(3.0) * ( A * U
+                           + B * V
+                           + C * W
+                           )
+                - ( A + B + C ) * ( U + V + W )
               ;
-
-              acc::sum_in_place(&_energy, &_loop_energy);
+              energy += F(2.0) * value / denominator * facjk * facij;
             } // i
           } // j
         } // k
+#if !defined(HAVE_CUDA)
       } // ii
     } // jj
   } // kk
-  const double real_part = acc::real(_energy);
-  acc::sum_in_place(energy, &real_part);
+#endif
+  return std::real(energy);
 }
 // Energy:2 ends here
 
 // [[file:~/cuda/atrip/atrip.org::*Energy][Energy:3]]
 // instantiate double
 template
-__MAYBE_GLOBAL__
-void getEnergyDistinct
-  ( DataFieldType<double> const epsabc
+double getEnergyDistinct
+  ( double const epsabc
   , size_t const No
-  , DataFieldType<double>* const epsi
-  , DataFieldType<double>* const Tijk
-  , DataFieldType<double>* const Zijk
-  , DataFieldType<double>* energy
+  , double* const epsi
+  , double* const Tijk
+  , double* const Zijk
   );
 
 template
-__MAYBE_GLOBAL__
-void getEnergySame
-  ( DataFieldType<double> const epsabc
+double getEnergySame
+  ( double const epsabc
   , size_t const No
-  , DataFieldType<double>* const epsi
-  , DataFieldType<double>* const Tijk
-  , DataFieldType<double>* const Zijk
-  , DataFieldType<double>* energy
+  , double* const epsi
+  , double* const Tijk
+  , double* const Zijk
   );
 
 // instantiate Complex
 template
-__MAYBE_GLOBAL__
-void getEnergyDistinct
-  ( DataFieldType<Complex> const epsabc
+double getEnergyDistinct
+  ( Complex const epsabc
   , size_t const No
-  , DataFieldType<Complex>* const epsi
-  , DataFieldType<Complex>* const Tijk
-  , DataFieldType<Complex>* const Zijk
-  , DataFieldType<double>* energy
+  , Complex* const epsi
+  , Complex* const Tijk
+  , Complex* const Zijk
   );
 
 template
-__MAYBE_GLOBAL__
-void getEnergySame
-  ( DataFieldType<Complex> const epsabc
+double getEnergySame
+  ( Complex const epsabc
   , size_t const No
-  , DataFieldType<Complex>* const epsi
-  , DataFieldType<Complex>* const Tijk
-  , DataFieldType<Complex>* const Zijk
-  , DataFieldType<double>* energy
+  , Complex* const epsi
+  , Complex* const Tijk
+  , Complex* const Zijk
   );
 // Energy:3 ends here
 
@@ -320,26 +379,18 @@ void getEnergySame
       const size_t ijk = i + j*No + k*NoNo;
 
 #ifdef HAVE_CUDA
-
-#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)
-
+#  define GO(__TPH, __VABIJ)                                          \
+    {                                                                 \
+      const DataFieldType<F> product                                  \
+            = cuda::multiply<DataFieldType<F>>((__TPH), (__VABIJ));   \
+      cuda::sum_in_place<DataFieldType<F>>(&Zijk[ijk], &product);     \
+    }
 #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
   }
@@ -401,15 +452,9 @@ void 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 a = abc[0], b = abc[1], c = abc[2]
-              , NoNo = No*No
-              ;
+
+    const size_t NoNo = No*No;
 
     DataFieldType<F>* Tijk = (DataFieldType<F>*)Tijk_;
 
@@ -454,7 +499,7 @@ void getEnergySame
                   )
 #define MAYBE_CONJ(_conj, _buffer)                                \
   do {                                                            \
-    acc::maybeConjugate<<<                                        \
+    cuda::maybeConjugate<<<                                       \
                                                                   \
                             Atrip::kernelDimensions.ooo.blocks,   \
                                                                   \
@@ -523,23 +568,23 @@ void getEnergySame
     F one{1.0}, m_one{-1.0}, zero{0.0};
     const size_t NoNoNo = No*NoNo;
 #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
+    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;
-    acc::zeroing<<<bs, ths>>>((DataFieldType<F>*)_t_buffer, NoNoNo);
-    acc::zeroing<<<bs, ths>>>((DataFieldType<F>*)_vhhh, NoNoNo);
+
+#if !defined(ATRIP_ONLY_DGEMM)
+    cuda::zeroing<<<bs, ths>>>((DataFieldType<F>*)_t_buffer, NoNoNo);
+    cuda::zeroing<<<bs, ths>>>((DataFieldType<F>*)_vhhh, NoNoNo);
 #endif
 
 #else
@@ -555,17 +600,15 @@ void getEnergySame
     // Set Tijk to zero
 #if defined(HAVE_CUDA) && !defined(ATRIP_ONLY_DGEMM)
     WITH_CHRONO("double:reorder",
-                acc::zeroing<<<bs, ths>>>((DataFieldType<F>*)Tijk,
+                cuda::zeroing<<<bs, ths>>>((DataFieldType<F>*)Tijk,
                                            NoNoNo);
                 )
-#endif
-
-#if !defined(HAVE_CUDA)
+#else
     WITH_CHRONO("double:reorder",
       for (size_t k = 0; k < NoNoNo; k++) {
         Tijk[k] = DataFieldType<F>{0.0};
        })
-#endif /* !defined(HAVE_CUDA) */
+#endif
 
 
 #if defined(ATRIP_ONLY_DGEMM)
@@ -573,7 +616,7 @@ void getEnergySame
 #undef REORDER
 #define MAYBE_CONJ(a, b) do {} while(0)
 #define REORDER(i, j, k) do {} while(0)
-#endif /* defined(ATRIP_ONLY_DGEMM) */
+#endif
 
     // HOLES
     WITH_CHRONO("doubles:holes",
@@ -657,16 +700,16 @@ void 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  /* defined(HAVE_CUDA) */
+#endif
   }
 
   #undef REORDER
@@ -717,7 +760,7 @@ void getEnergySame
       }
 
     }
-#endif /* defined(ATRIP_USE_DGEMM) */
+#endif
   }
 
 
@@ -749,12 +792,6 @@ void getEnergySame
     , DataPtr<double> const TBChh
     // -- TIJK
     , DataFieldType<double>* Tijk
-#if defined(HAVE_CUDA)
-     // -- tmp buffers
-    , DataFieldType<double>* _t_buffer
-    , DataFieldType<double>* _vhhh 
-#endif
-
     );
 
   template
@@ -783,12 +820,6 @@ void 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