include/atrip/Operations.hpp

@@ -24,6 +24,15 @@ 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__

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

src/atrip/Complex.cxx

@@ -21,6 +21,11 @@ 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,8 +13,6 @@
 // limitations under the License.
 
 // [[file:~/cuda/atrip/atrip.org::*Prolog][Prolog:2]]
-#include <cstring>
-
 #include<atrip/Equations.hpp>
 
 #include<atrip/CUDA.hpp>
@@ -28,7 +26,7 @@ namespace atrip {
 #if defined(HAVE_CUDA)
 #define FOR_K()                                             \
   const size_t k = blockIdx.x * blockDim.x + threadIdx.x; \
-  size_t idx = k*size*size;
+  size_t idx = 0;
 #else
 #define FOR_K() for (size_t k=0, idx=0; k < size; k++)
 #endif
@@ -582,8 +580,13 @@ void getEnergySame
                   )
 #define MAYBE_CONJ(_conj, _buffer)                                \
   do {                                                            \
-    acc::maybeConjugate<<<1, 1                                    \
-                       >>>((DataFieldType<F>*)_conj,              \
+    acc::maybeConjugate<<<                                        \
+                                                                  \
+                            Atrip::kernelDimensions.ooo.blocks,   \
+                                                                  \
+                              Atrip::kernelDimensions.ooo.threads \
+                                                                  \
+                        >>>((DataFieldType<F>*)_conj,             \
                             (DataFieldType<F>*)_buffer,           \
                             NoNoNo);                              \
   } while (0)
@@ -634,39 +637,61 @@ void getEnergySame
                   _t_buffer,                    \
                   (int const*)&NoNo             \
                   )
-#define MAYBE_CONJ(_conj, _buffer)               \
-  acc::maybeConjugate((DataFieldType<F>*)_conj,  \
-                      (DataFieldType<F>*)_buffer,\
-                      NoNoNo);
+#define MAYBE_CONJ(_conj, _buffer)              \
+  do {                                          \
+    for (size_t __i = 0; __i < NoNoNo; ++__i) { \
+      _conj[__i]                                \
+        = maybeConjugate<F>(_buffer[__i]);      \
+    }                                           \
+  } while (0)
 #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)
-    {
-    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));
-                )
-    }
+    // 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);
 #endif
 
 #else
     DataFieldType<F>* _t_buffer = (DataFieldType<F>*)malloc(NoNoNo * sizeof(F));
     DataFieldType<F>* _vhhh = (DataFieldType<F>*)malloc(NoNoNo * sizeof(F));
-    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 */
+    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) */
 
 
 #if defined(ATRIP_ONLY_DGEMM)