gallo/atrip
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Compare commits

base: gallo:122329eca7bdd3a4ff8c2043c0a27fd1e5176813
Branches Tags
gallo:cuda
gallo:cuda-energy
gallo:openacc
gallo:newtuples
gallo:cuda-test-doubles
gallo:cuda-energy-working
..
compare: gallo:249f1c0b51fc7fff305e0b0ceab5ab353ef78b2c
Branches Tags
gallo:cuda-energy
gallo:cuda
gallo:openacc
gallo:newtuples
gallo:cuda-test-doubles
gallo:cuda-energy-working

No commits in common. "122329eca7bdd3a4ff8c2043c0a27fd1e5176813" and "249f1c0b51fc7fff305e0b0ceab5ab353ef78b2c" have entirely different histories.
122329eca7 ... 249f1c0b51

7 changed files with 219 additions and 511 deletions
Show Stats Expand all files Collapse all files

11
include/atrip/CUDA.hpp
View File

@ -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__; \

15
include/atrip/Equations.hpp
View File

@ -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

171
include/atrip/Operations.hpp
View File

@ -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

6
include/atrip/Slice.hpp
View File

@ -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_)

53
include/atrip/SliceUnion.hpp
View File

@ -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

100
src/atrip/Atrip.cxx
View File

@ -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,31 @@ 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) */
/*
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 defined(HAVE_OCD) || defined(ATRIP_PRINT_TUPLES)
tupleEnergies[abc] = host_tuple_energy;
tupleEnergies[abc] = tupleEnergy;
#endif
energy += host_tuple_energy;
energy += tupleEnergy;
}

370
src/atrip/Equations.cxx
View File

@ -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,162 @@ 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
double 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.};
F energy(0.);
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};
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))
, 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
} // ii
} // jj
} // kk
const double real_part = acc::real(_energy);
acc::sum_in_place(energy, &real_part);
return std::real(energy);
}
template <typename F>
__MAYBE_GLOBAL__
void getEnergySame
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.);
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 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))
, 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
} // ii
} // jj
} // kk
const double real_part = acc::real(_energy);
acc::sum_in_place(energy, &real_part);
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 +360,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)
{ \
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 +433,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 +480,7 @@ void getEnergySame
)
#define MAYBE_CONJ(_conj, _buffer) \
do { \
acc::maybeConjugate<<< \
cuda::maybeConjugate<<< \
\
Atrip::kernelDimensions.ooo.blocks, \
\
@ -523,23 +549,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 +581,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 +597,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 +681,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 +741,7 @@ void getEnergySame
}
}
#endif /* defined(ATRIP_USE_DGEMM) */
#endif
}
@ -749,12 +773,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 +801,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