include/atrip/Operations.hpp

@@ -38,8 +38,6 @@ namespace acc {
   __MAYBE_DEVICE__ __MAYBE_HOST__ __INLINE__
   F maybeConjugateScalar(const F &a) { return a; }
 
-  // TODO: instantiate for std::complex<double>
-
 #if defined(HAVE_CUDA)
   template <>
   __MAYBE_DEVICE__ __MAYBE_HOST__ __INLINE__

include/atrip/SliceUnion.hpp

@@ -200,7 +200,7 @@ template <typename F=double>
                            : Slice<F>::Fetch
                            ;
           if (blank.info.state == Slice<F>::SelfSufficient) {
-#if defined(HAVE_CUDA) && !defined(ATRIP_SOURCES_IN_GPU)
+#if defined(HAVE_CUDA)
             const size_t _size = sizeof(F) * sliceSize;
             // TODO: this is code duplication with downstairs
             if (freePointers.size() == 0) {
@@ -221,6 +221,7 @@ template <typename F=double>
                                                          (void*)SOURCES_DATA(sources[from.source]),
                                                          sizeof(F) * sliceSize));
                         ))
+
 #else
             blank.data = SOURCES_DATA(sources[from.source]);
 #endif
@@ -464,30 +465,30 @@ template <typename F=double>
 
 #if defined(HAVE_CUDA)
       LOG(1,"Atrip") << "warming communication up " << slices.size() << "\n";
-      WITH_CHRONO("cuda:warmup",
+      WITH_CHRONO("cuda:warmup",
-                  int nRanks=Atrip::np, requestCount=0;
+                   int nRanks=Atrip::np, requestCount=0;
-                  int nSends=sliceBuffers.size()*nRanks;
+                   int nSends=sliceBuffers.size()*nRanks;
-                  MPI_Request *requests = (MPI_Request*) malloc(nSends*2 * sizeof(MPI_Request));
+                   MPI_Request *requests = (MPI_Request*) malloc(nSends*2 * sizeof(MPI_Request));
-                  MPI_Status *statuses = (MPI_Status*) malloc(nSends*2 * sizeof(MPI_Status));
+                   MPI_Status *statuses = (MPI_Status*) malloc(nSends*2 * sizeof(MPI_Status));
-                  for (int sliceId=0; sliceId<sliceBuffers.size(); sliceId++){
+                   for (int sliceId=0; sliceId<sliceBuffers.size(); sliceId++){
-                    for (int rankId=0; rankId<nRanks; rankId++){
+                     for (int rankId=0; rankId<nRanks; rankId++){
-                      MPI_Isend((void*)SOURCES_DATA(sources[0]),
+                       MPI_Isend((void*)SOURCES_DATA(sources[0]),
-                                sliceSize,
+                                 sliceSize,
-                                traits::mpi::datatypeOf<F>(),
+                                 traits::mpi::datatypeOf<F>(),
-                                rankId,
+                                 rankId,
-                                100,
+                                 100,
-                                universe,
+                                 universe,
-                                &requests[requestCount++]);
+                                 &requests[requestCount++]);
-                      MPI_Irecv((void*)sliceBuffers[sliceId],
+                       MPI_Irecv((void*)sliceBuffers[sliceId],
-                                sliceSize,
+                                 sliceSize,
-                                traits::mpi::datatypeOf<F>(),
+                                 traits::mpi::datatypeOf<F>(),
-                                rankId,
+                                 rankId,
-                                100,
+                                 100,
-                                universe,
+                                 universe,
-                                &requests[requestCount++]);
+                                 &requests[requestCount++]);
-                    }
+                     }
-                  }
+                   }
-                  MPI_Waitall(nSends*2, requests, statuses);
+                   MPI_Waitall(nSends*2, requests, statuses);
                   )
 #endif
 
@@ -577,11 +578,8 @@ template <typename F=double>
 #    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((void*)slice.data,
+        MPI_Irecv(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
-  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
@@ -693,7 +693,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),
@@ -733,30 +733,18 @@ Atrip::Output Atrip::run(Atrip::Input<F> const& in) {
                                 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);
                   })
 
@@ -903,5 +891,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/Equations.cxx

@@ -25,8 +25,11 @@ namespace atrip {
 
 #if defined(HAVE_CUDA)
 #define FOR_K()                                             \
-  const size_t k = blockIdx.x * blockDim.x + threadIdx.x; \
+  for (size_t kmin = blockIdx.x * blockDim.x + threadIdx.x, \
-  size_t idx = 0;
+         k = kmin,                                          \
+         idx = kmin * size * size * size;                   \
+       k < (kmin < size) ? kmin + 1 : size;                 \
+       k++)
 #else
 #define FOR_K() for (size_t k=0, idx=0; k < size; k++)
 #endif
@@ -99,7 +102,6 @@ namespace atrip {
 #  define MIN(a, b) std::min((a), (b))
 #endif
 
-#if defined(ATRIP_NEW_ENERGY)
 
 // [[file:~/cuda/atrip/atrip.org::*Energy][Energy:2]]
 template <typename F>
@@ -248,131 +250,6 @@ void getEnergySame
 }
 // Energy:2 ends here
 
-#else
-
-// [[file:~/cuda/atrip/atrip.org::*Energy][Energy:2]]
-template <typename F>
-__MAYBE_GLOBAL__
-void getEnergyDistinct
-  ( F const epsabc
-  , size_t const No
-  , F* const epsi
-  , F* const Tijk
-  , F* const Zijk
-  , double* _energy
-  ) {
-  constexpr size_t blockSize=16;
-  F energy(0.);
-  for (size_t kk=0; kk<No; kk+=blockSize){
-    const size_t kend( MIN(No, kk+blockSize) );
-    for (size_t jj(kk); jj<No; jj+=blockSize){
-      const size_t jend( MIN( No, jj+blockSize) );
-      for (size_t ii(jj); ii<No; ii+=blockSize){
-        const size_t iend( MIN( No, ii+blockSize) );
-        for (size_t k(kk); k < kend; k++){
-          const F ek(epsi[k]);
-          const size_t jstart = jj > k ? jj : k;
-          for (size_t j(jstart); j < jend; j++){
-            F const ej(epsi[j]);
-            F const facjk = j == k ? F(0.5) : F(1.0);
-            size_t istart = ii > j ? ii : j;
-            for (size_t i(istart); i < iend; i++){
-              const F
-                  ei(epsi[i])
-                , facij = i == j ? F(0.5) : F(1.0)
-                , denominator(epsabc - ei - ej - ek)
-                , U(Zijk[i + No*j + No*No*k])
-                , V(Zijk[i + No*k + No*No*j])
-                , W(Zijk[j + No*i + No*No*k])
-                , X(Zijk[j + No*k + No*No*i])
-                , Y(Zijk[k + No*i + No*No*j])
-                , Z(Zijk[k + No*j + No*No*i])
-                , A(acc::maybeConjugateScalar<F>(Tijk[i + No*j + No*No*k]))
-                , B(acc::maybeConjugateScalar<F>(Tijk[i + No*k + No*No*j]))
-                , C(acc::maybeConjugateScalar<F>(Tijk[j + No*i + No*No*k]))
-                , D(acc::maybeConjugateScalar<F>(Tijk[j + No*k + No*No*i]))
-                , E(acc::maybeConjugateScalar<F>(Tijk[k + No*i + No*No*j]))
-                , _F(acc::maybeConjugateScalar<F>(Tijk[k + No*j + No*No*i]))
-                , value
-                  = 3.0 * ( A * U
-                            + B * V
-                            + C * W
-                            + D * X
-                            + E * Y
-                            + _F * Z )
-                 + ( ( U + X + Y )
-                   - 2.0 * ( V + W + Z )
-                   ) * ( A + D + E )
-                 + ( ( V + W + Z )
-                   - 2.0 * ( U + X + Y )
-                   ) * ( B + C + _F )
-                ;
-              energy += 2.0 * value / denominator * facjk * facij;
-            } // i
-          } // j
-        } // k
-      } // ii
-    } // jj
-  } // kk
-  *_energy = acc::real(energy);
-}
-
-
-template <typename F>
-__MAYBE_GLOBAL__
-void getEnergySame
-  ( F const epsabc
-  , size_t const No
-  , F* const epsi
-  , F* const Tijk
-  , F* const Zijk
-  , double* _energy
-  ) {
-  constexpr size_t blockSize = 16;
-  F energy = F(0.);
-  for (size_t kk=0; kk<No; kk+=blockSize){
-    const size_t kend( MIN( kk+blockSize, No) );
-    for (size_t jj(kk); jj<No; jj+=blockSize){
-      const size_t jend( MIN( jj+blockSize, No) );
-      for (size_t ii(jj); ii<No; ii+=blockSize){
-        const size_t iend( MIN( ii+blockSize, No) );
-        for (size_t k(kk); k < kend; k++){
-          const F ek(epsi[k]);
-          const size_t jstart = jj > k ? jj : k;
-          for(size_t j(jstart); j < jend; j++){
-            const F facjk( j == k ? F(0.5) : F(1.0));
-            const F ej(epsi[j]);
-            const size_t istart = ii > j ? ii : j;
-            for(size_t i(istart); i < iend; i++){
-              const F
-                ei(epsi[i])
-              , facij ( i==j ? F(0.5) : F(1.0))
-              , denominator(epsabc - ei - ej - ek)
-              , U(Zijk[i + No*j + No*No*k])
-              , V(Zijk[j + No*k + No*No*i])
-              , W(Zijk[k + No*i + No*No*j])
-              , A(acc::maybeConjugateScalar<F>(Tijk[i + No*j + No*No*k]))
-              , B(acc::maybeConjugateScalar<F>(Tijk[j + No*k + No*No*i]))
-              , C(acc::maybeConjugateScalar<F>(Tijk[k + No*i + No*No*j]))
-              , value
-                = F(3.0) * ( A * U
-                           + B * V
-                           + C * W
-                           )
-                - ( A + B + C ) * ( U + V + W )
-              ;
-              energy += F(2.0) * value / denominator * facjk * facij;
-            } // i
-          } // j
-        } // k
-      } // ii
-    } // jj
-  } // kk
-  *_energy = acc::real(energy);
-}
-// Energy:2 ends here
-#endif /* defined(ATRIP_NEW_ENERGY) */
-
 // [[file:~/cuda/atrip/atrip.org::*Energy][Energy:3]]
 // instantiate double
 template
@@ -397,8 +274,6 @@ void getEnergySame
   , DataFieldType<double>* energy
   );
 
-// TODO: put this back in
-#if defined(ATRIP_WITH_COMPLEX)
 // instantiate Complex
 template
 __MAYBE_GLOBAL__
@@ -422,7 +297,6 @@ void getEnergySame
   , DataFieldType<double>* energy
   );
 // Energy:3 ends here
-#endif
 
 // [[file:~/cuda/atrip/atrip.org::*Singles%20contribution][Singles contribution:2]]
   template <typename F> __MAYBE_GLOBAL__
@@ -542,7 +416,7 @@ void getEnergySame
 #if defined(ATRIP_USE_DGEMM)
 #if defined(HAVE_CUDA)
 #define REORDER(__II, __JJ, __KK)               \
-  reorder<<<1, No>>>(reorder_proxy<           \
+  reorder<<<bs, ths>>>(reorder_proxy<           \
                        DataFieldType<F>,        \
                        __II ## __JJ ## __KK     \
                        >{},                     \
@@ -659,15 +533,11 @@ void getEnergySame
 //                        cuMemAlloc((CUdeviceptr*)&_vhhh,
 //                                   NoNoNo * sizeof(DataFieldType<F>)));
 //                )
-#if !defined(ATRIP_ONLY_DGEMM)
+//    const size_t
-    // we still have to zero this
+//      bs = Atrip::kernelDimensions.ooo.blocks,
-    const size_t
+//      ths = Atrip::kernelDimensions.ooo.threads;
-      bs = Atrip::kernelDimensions.ooo.blocks,
+    //cuda::zeroing<<<bs, ths>>>((DataFieldType<F>*)_t_buffer, NoNoNo);
-      ths = Atrip::kernelDimensions.ooo.threads;
+    //cuda::zeroing<<<bs, ths>>>((DataFieldType<F>*)_vhhh, NoNoNo);
-    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));
@@ -684,14 +554,12 @@ void getEnergySame
                 acc::zeroing<<<bs, ths>>>((DataFieldType<F>*)Tijk,
                                            NoNoNo);
                 )
-#endif
+#else
-
-#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) */
+#endif /* defined(HAVE_CUDA) && !defined(ATRIP_ONLY_DGEMM) */
 
 
 #if defined(ATRIP_ONLY_DGEMM)