diff --git a/Compute.cpp b/Compute.cpp index 2a90187d..04282c4a 100644 --- a/Compute.cpp +++ b/Compute.cpp @@ -9,6 +9,10 @@ #include "Space.h" #include "gravity.h" +#ifdef CUDA +#include "PEList.h" +#endif + int decodeReqID(int reqID); void Compute::setOpt(Opt *_opt){ @@ -306,6 +310,7 @@ void GravityCompute::recvdParticles(ExternalGravityParticle *part,int num,int ch #ifdef BENCHMARK_TIME_COMPUTE computeTimePart += CmiWallTimer() - startTime; #endif + CkAssert(tp->particleInterRemote != NULL); tp->particleInterRemote[chunk] += computed * num; tp->finishBucket(reqIDlist); CkAssert(state->counterArrays[1][chunk] >= 0); @@ -341,6 +346,7 @@ void GravityCompute::nodeRecvdEvent(TreePiece *owner, int chunk, State *state, i state->counterArrays[1][chunk] --; CkAssert(state->counterArrays[1][chunk] >= 0); if (state->counterArrays[1][chunk] == 0) { + CkAssert(owner->particleInterRemote != NULL); cacheGravPart[CkMyPe()].finishedChunk(chunk, owner->particleInterRemote[chunk]); #ifdef CHECK_WALK_COMPLETIONS CkPrintf("[%d] finishedChunk %d GravityCompute::nodeRecvdEvent\n", owner->getIndex(), chunk); @@ -377,20 +383,18 @@ void ListCompute::nodeRecvdEvent(TreePiece *owner, int chunk, State *state, int // flush the interactions remaining in the state DoubleWalkState *ds = (DoubleWalkState *)state; - if(ds->nodeLists.totalNumInteractions > 0){ - sendNodeInteractionsToGpu(ds, owner); - resetCudaNodeState(ds); - } - if(ds->particleLists.totalNumInteractions > 0){ - sendPartInteractionsToGpu(ds, owner); - resetCudaPartState(ds); - } + treeProxy[owner->getIndex()].flushInteractionsToGpu((intptr_t) this, + (intptr_t) ds); } #endif #if COSMO_PRINT_BK > 1 CkPrintf("[%d] FINISHED CHUNK %d from nodeRecvdEvent\n", owner->getIndex(), chunk); #endif + if (owner->particleInterRemote == NULL) { + CkPrintf("ERROR [%d] TP %d particleInterRemote NULL at chunk %d (ListCompute::nodeRecvdEvent)\n", CkMyPe(), owner->getIndex(), chunk); + CkAbort("particleInterRemote NULL"); + } cacheGravPart[CkMyPe()].finishedChunk(chunk, owner->particleInterRemote[chunk]); #ifdef CHECK_WALK_COMPLETIONS CkPrintf("[%d] finishedChunk %d ListCompute::nodeRecvdEvent\n", owner->getIndex(), chunk); @@ -948,19 +952,17 @@ void ListCompute::recvdParticles(ExternalGravityParticle *part,int num,int chunk if (remainingChunk == 0) { #ifdef CUDA if (!bUseCpu) { - if(state->nodeLists.totalNumInteractions > 0){ - sendNodeInteractionsToGpu(state, tp); - resetCudaNodeState(state); - } - if(state->particleLists.totalNumInteractions > 0){ - sendPartInteractionsToGpu(state, tp); - resetCudaPartState(state); - } + treeProxy[tp->getIndex()].flushInteractionsToGpu((intptr_t) this, + (intptr_t) state); } #endif #if COSMO_PRINT_BK > 1 CkPrintf("[%d] FINISHED CHUNK %d from recvdParticles\n", tp->getIndex(), chunk); #endif + if (tp->particleInterRemote == NULL) { + CkPrintf("ERROR [%d] TP %d particleInterRemote NULL at chunk %d (ListCompute::recvdParticles)\n", CkMyPe(), tp->getIndex(), chunk); + CkAbort("particleInterRemote NULL"); + } cacheGravPart[CkMyPe()].finishedChunk(chunk, tp->particleInterRemote[chunk]); #ifdef CHECK_WALK_COMPLETIONS CkPrintf("[%d] finishedChunk %d ListCompute::recvdParticles\n", tp->getIndex(), chunk); @@ -1057,19 +1059,13 @@ CudaRequest *GenericList::serialize(TreePiece *tp){ int *affectedBuckets = NULL; if(totalNumInteractions > 0){ -#ifdef PINNED_HOST_MEMORY - allocatePinnedHostMemory((void **)&flatlists, totalNumInteractions*sizeof(T)); - allocatePinnedHostMemory((void **)&markers, (numFilledBuckets+1)*sizeof(int)); - allocatePinnedHostMemory((void **)&starts, (numFilledBuckets)*sizeof(int)); - allocatePinnedHostMemory((void **)&sizes, (numFilledBuckets)*sizeof(int)); -#else flatlists = (T *) malloc(totalNumInteractions*sizeof(T)); markers = (int *) malloc((numFilledBuckets+1)*sizeof(int)); starts = (int *) malloc(numFilledBuckets*sizeof(int)); sizes = (int *) malloc(numFilledBuckets*sizeof(int)); -#endif affectedBuckets = new int[numFilledBuckets]; + CkAssert(tp->bGPUBufferFilled == 1); // Needed for getBucketParameters() // populate flat lists int listslen = lists.length(); for(int i = 0; i < listslen; i++){ @@ -1077,12 +1073,7 @@ CudaRequest *GenericList::serialize(TreePiece *tp){ if(listilen > 0){ memcpy(&flatlists[listpos], lists[i].getVec(), listilen*sizeof(T)); markers[curbucket] = listpos; - if(tp->largePhase()){ - getBucketParameters(tp, i, starts[curbucket], sizes[curbucket]); - } - else{ - getActiveBucketParameters(tp, i, starts[curbucket], sizes[curbucket]); - } + getBucketParameters(tp, i, starts[curbucket], sizes[curbucket]); affectedBuckets[curbucket] = i; listpos += listilen; curbucket++; @@ -1090,6 +1081,7 @@ CudaRequest *GenericList::serialize(TreePiece *tp){ } markers[numFilledBuckets] = listpos; CkAssert(listpos == totalNumInteractions); + CkAssert(starts[0] >= 0); // Ensure tree serialization has already finished } CudaRequest *request = new CudaRequest; @@ -1135,8 +1127,6 @@ CudaRequest *GenericList::serialize(TreePiece *tp){ } } - //CkPrintf("[%d] Offloading %d particle interactions\n", tp->getIndex(), numParticleInteractions); - // create flat lists and associated data structures // allocate memory and flatten lists only if there // are interactions to transfer. @@ -1149,17 +1139,10 @@ CudaRequest *GenericList::serialize(TreePiece *tp){ int *affectedBuckets = NULL; if(totalNumInteractions > 0){ -#ifdef PINNED_HOST_MEMORY - allocatePinnedHostMemory((void **)&flatlists, numParticleInteractions*sizeof(ILCell)); - allocatePinnedHostMemory((void **)&markers, (numFilledBuckets+1)*sizeof(int)); - allocatePinnedHostMemory((void **)&starts, (numFilledBuckets)*sizeof(int)); - allocatePinnedHostMemory((void **)&sizes, (numFilledBuckets)*sizeof(int)); -#else flatlists = (ILCell *) malloc(numParticleInteractions*sizeof(ILCell)); markers = (int *) malloc((numFilledBuckets+1)*sizeof(int)); starts = (int *) malloc(numFilledBuckets*sizeof(int)); sizes = (int *) malloc(numFilledBuckets*sizeof(int)); -#endif affectedBuckets = new int[numFilledBuckets]; // populate flat lists @@ -1185,18 +1168,14 @@ CudaRequest *GenericList::serialize(TreePiece *tp){ listpos++; } } - if(tp->largePhase()){ - getBucketParameters(tp, i, starts[curbucket], sizes[curbucket]); - } - else{ - getActiveBucketParameters(tp, i, starts[curbucket], sizes[curbucket]); - } + getBucketParameters(tp, i, starts[curbucket], sizes[curbucket]); affectedBuckets[curbucket] = i; curbucket++; } } markers[numFilledBuckets] = listpos; CkAssert(listpos == numParticleInteractions); + CkAssert(starts[0] >= 0); // Ensure tree serialization has already finished } CudaRequest *request = new CudaRequest; @@ -1391,146 +1370,6 @@ template int calcParticleForces(TreePiece *tp, int b, int activeRung return computed; } -#ifdef GPU_LOCAL_TREE_WALK - -// XXX This appears to be identical to cudaCallback(), I think it can -// be deleted --trq. -void cudaCallbackForAllBuckets(void *param, void *msg) { - CudaRequest *data = (CudaRequest *)param; - int *affectedBuckets = data->affectedBuckets; - TreePiece *tp = (TreePiece*)data->tp; - DoubleWalkState *state = (DoubleWalkState *)data->state; - int bucket; - - int numBucketsDone = data->numBucketsPlusOne-1; - - // bucket computations finished - // - for(int i = 0; i < numBucketsDone; i++){ - bucket = affectedBuckets[i]; - state->counterArrays[0][bucket]--; - tp->finishBucket(bucket); - } - - // free data structures - if(numBucketsDone > 0){ - delete [] data->affectedBuckets; - } - freePinnedHostMemory(data->list); - freePinnedHostMemory(data->bucketMarkers); - freePinnedHostMemory(data->bucketStarts); - freePinnedHostMemory(data->bucketSizes); - -#ifdef CUDA_PRINT_ERRORS - printf("cudaCallbackForAllBuckets: %s\n", cudaGetErrorString( cudaGetLastError() ) ); -#endif - - delete ((CkCallback *)data->cb); - delete data; -} - -/** - * This function is designed to send an ignition signal to the GPU manager. - * To make minor change to existing ChaNGa code, we mimic a nodeGravityCompute - * request. This request will eventually call our GPU local tree walk kernel. - */ -void ListCompute::sendLocalTreeWalkTriggerToGpu(State *state, TreePiece *tp, - int activeRung, int startBucket, int endBucket) { - int numFilledBuckets = 0; - for (int i = startBucket; i < endBucket; ++i) { - if (tp->bucketList[i]->rungs >= activeRung) { - ++numFilledBuckets; - } - } - - // No necessary to call GPU kernel if there is no active bucket in - // current tree piece - if (numFilledBuckets == 0) { - return; - } - - int *affectedBuckets = new int[numFilledBuckets]; - - // Set up a series of dummy parameters to match existing function interfaces - int dummyTotalNumInteractions = 1; - int dummyCurBucket = 0; - ILCell *dummyFlatlists = NULL; - int *dummyNodeMarkers = NULL; - int *dummyStarts = NULL; - int *dummySizes = NULL; - - // XXX I think this can be deleted --trq. -#ifdef PINNED_HOST_MEMORY - allocatePinnedHostMemory((void **)&dummyFlatlists, dummyTotalNumInteractions * - sizeof(ILCell)); - allocatePinnedHostMemory((void **)&dummyNodeMarkers, (numFilledBuckets+1) * - sizeof(int)); - allocatePinnedHostMemory((void **)&dummyStarts, numFilledBuckets * - sizeof(int)); - allocatePinnedHostMemory((void **)&dummySizes, numFilledBuckets * sizeof(int)); -#else - dummyFlatlists = (ILCell *) malloc(dummy_totalNumInteractions*sizeof(ILCell)); - dummyNodeMarkers = (int *) malloc((numFilledBuckets+1)*sizeof(int)); - dummyStarts = (int *) malloc(numFilledBuckets*sizeof(int)); - dummySizes = (int *) malloc(numFilledBuckets*sizeof(int)); -#endif - -// No need to memset the interaction list array since we're not using it at all -// And, we can't directly memset it. - ILCell temp_ilc; - memcpy(&dummyFlatlists[0], &temp_ilc, dummyTotalNumInteractions * - sizeof(ILCell)); - for (int i = startBucket; i < endBucket; ++i) { - if (tp->bucketList[i]->rungs >= activeRung) { - ((DoubleWalkState *)state)->counterArrays[0][i] ++; - dummyNodeMarkers[dummyCurBucket] = tp->bucketList[i]->nodeArrayIndex; - int tempNum = 0; - memcpy(&dummyStarts[dummyCurBucket], &tempNum, sizeof(int)); - memcpy(&dummySizes[dummyCurBucket], &tempNum, sizeof(int)); - affectedBuckets[dummyCurBucket] = i; - dummyCurBucket++; - } - } - - CudaRequest *request = new CudaRequest; - - request->d_localMoments = tp->d_localMoments; - request->d_localParts = tp->d_localParts; - request->d_localVars = tp->d_localVars; - request->sMoments = tp->sMoments; - request->sCompactParts = tp->sCompactParts; - request->sVarParts = tp->sVarParts; - request->stream = tp->stream; - - request->numBucketsPlusOne = numFilledBuckets+1; - request->affectedBuckets = affectedBuckets; - request->tp = (void *)tp; - request->state = (void *)state; - request->node = true; - request->remote = false; - request->firstParticle = tp->FirstGPUParticleIndex; - request->lastParticle = tp->LastGPUParticleIndex; - // In DataManager serializes the local tree so that the root of the local tree - // will always be the 0th element in the moments array. - request->rootIdx = 0; - request->theta = theta; - request->thetaMono = thetaMono; - request->nReplicas = tp->nReplicas; - request->fperiod = tp->fPeriod.x; - request->fperiodY = tp->fPeriod.y; - request->fperiodZ = tp->fPeriod.z; - request->cb = new CkCallback(cudaCallbackForAllBuckets, request); - - request->list = (void *)dummyFlatlists; - request->bucketMarkers = dummyNodeMarkers; - request->bucketStarts = dummyStarts; - request->bucketSizes = dummySizes; - request->numInteractions = dummyTotalNumInteractions; - - TreePieceCellListDataTransferLocal(request); -} -#endif //GPU_LOCAL_TREE_WALK - /// @brief Check for computation /// Computation can be done on buckets indexed from start to end /// @param state_ State to be checked @@ -1544,9 +1383,6 @@ void ListCompute::stateReady(State *state_, TreePiece *tp, int chunk, int start, DoubleWalkState *state = (DoubleWalkState *)state_; #ifdef CUDA - bool resume = state->resume; - int numActiveBuckets = tp->numActiveBuckets; - // for local particles //std::map &lpref = tp->dm->getLocalPartsOnGpuTable(); // for cached particles @@ -1666,7 +1502,6 @@ void ListCompute::stateReady(State *state_, TreePiece *tp, int chunk, int start, continue; } #endif - int computed = 0; // check whether it is already on the gpu int index = node->nodeArrayIndex; @@ -1821,82 +1656,33 @@ void ListCompute::stateReady(State *state_, TreePiece *tp, int chunk, int start, // local particles if(hasLocalLists){ - if(tp->largePhase()){ - for(int level = 0; level <= maxlevel; level++){ - CkVec &lpilist = state->lplists[level]; - for(int i = 0; i < lpilist.length(); i++){ - LocalPartInfo &lpi = lpilist[i]; + for(int level = 0; level <= maxlevel; level++){ + CkVec &lpilist = state->lplists[level]; + for(int i = 0; i < lpilist.length(); i++){ + LocalPartInfo &lpi = lpilist[i]; #if defined CHANGA_REFACTOR_WALKCHECK_INTERLIST || defined CHANGA_REFACTOR_PRINT_INTERACTIONS - NodeKey key = lpi.key; + NodeKey key = lpi.key; #endif #if defined CHANGA_REFACTOR_WALKCHECK_INTERLIST - tp->addToBucketChecklist(b, key); - tp->combineKeys(key, b); + tp->addToBucketChecklist(b, key); + tp->combineKeys(key, b); #endif - int gpuIndex = lpi.nd->bucketArrayIndex; + int gpuIndex = lpi.nd->bucketArrayIndex; #ifdef CHANGA_REFACTOR_PRINT_INTERACTIONS - { - int thisIndex = tp->getIndex(); - if(b == TEST_BUCKET && tp->getIndex() == TEST_TP){ - Vector3D &vec = lpi.offset; - CkPrintf("[%d]: remote: %d resume: %d bucket %d with local part %ld (%d), (%1.0f,%1.0f,%1.0f)\n", thisIndex, getOptType() == Remote, state->resume, b, key, gpuIndex, vec.x, vec.y, vec.z); - } - } -#endif - - if(gpuIndex >= 0){ - - // put bucket in interaction list - Vector3D &off = lpi.offset; - ILPart tilp(gpuIndex, encodeOffset(0, off.x, off.y, off.z), lpi.numParticles); - state->particleLists.push_back(b, tilp, state, tp); - if(state->partOffloadReady()){ - // enough nodes to offload - sendPartInteractionsToGpu(state, tp); - resetCudaPartState(state); - } + { + int thisIndex = tp->getIndex(); + if(b == TEST_BUCKET && tp->getIndex() == TEST_TP){ + Vector3D &vec = lpi.offset; + CkPrintf("[%d]: remote: %d resume: %d bucket %d with local part %ld (%d), (%1.0f,%1.0f,%1.0f)\n", thisIndex, getOptType() == Remote, state->resume, b, key, gpuIndex, vec.x, vec.y, vec.z); } } - } - } - else{ // small phase, need to attach particle data to state->particles - for(int level = 0; level <= maxlevel; level++){ - CkVec &lpilist = state->lplists[level]; - for(int i = 0; i < lpilist.length(); i++){ - LocalPartInfo &lpi = lpilist[i]; -#if defined CHANGA_REFACTOR_WALKCHECK_INTERLIST || defined CHANGA_REFACTOR_PRINT_INTERACTIONS - NodeKey key = lpi.key; #endif -#if defined CHANGA_REFACTOR_WALKCHECK_INTERLIST - tp->addToBucketChecklist(b, key); - tp->combineKeys(key, b); -#endif - - int gpuIndex = lpi.nd->bucketArrayIndex; - if(gpuIndex < 0){ - CkAssert(state->particles != NULL); - gpuIndex = state->particles->length(); - lpi.nd->bucketArrayIndex = gpuIndex; - state->markedBuckets.push_back(lpi.nd); - for(int j = 0; j < lpi.numParticles; j++){ - state->particles->push_back(CompactPartData(lpi.particles[j])); - } - } -#ifdef CHANGA_REFACTOR_PRINT_INTERACTIONS - { - int thisIndex = tp->getIndex(); - if(b == TEST_BUCKET && tp->getIndex() == TEST_TP){ - Vector3D &vec = lpi.offset; - CkPrintf("[%d]: remote: %d resume: %d bucket %d with local part %ld (%d), (%1.0f,%1.0f,%1.0f)\n", thisIndex, getOptType() == Remote, state->resume, b, key, gpuIndex, vec.x, vec.y, vec.z); - } - } -#endif - CkAssert(gpuIndex >= 0); + if(gpuIndex >= 0){ - // put bucket in interaction list + // put bucket in interaction list Vector3D &off = lpi.offset; ILPart tilp(gpuIndex, encodeOffset(0, off.x, off.y, off.z), lpi.numParticles); state->particleLists.push_back(b, tilp, state, tp); @@ -1978,7 +1764,6 @@ void ListCompute::initCudaState(DoubleWalkState *state, int numBuckets, int node /// @brief Reset node array after interactions have been sent to the GPU. void ListCompute::resetCudaNodeState(DoubleWalkState *state){ - GenericTreeNode *tmp; state->nodeLists.reset(); if(state->nodes){ state->nodes->length() = 0; @@ -1995,70 +1780,37 @@ void ListCompute::resetCudaPartState(DoubleWalkState *state){ } } -/// A group of node or particle interactions has been completed by the GPU. -void cudaCallback(void *param, void *msg){ - CudaRequest *data = (CudaRequest *)param; - // Paranoid about data corruption here. - CkAssert(((CkCallback *)data->cb)->type == CkCallback::callbackType::callCFn); - CkAssert(((CkCallback *)data->cb)->d.cfn.fn == cudaCallback); - CkAssert(data->numBucketsPlusOne > 0); - CkAssert(data->numInteractions >= 0); - - int *affectedBuckets = data->affectedBuckets; - TreePiece *tp = (TreePiece*)data->tp; - DoubleWalkState *state = (DoubleWalkState *)data->state; - int bucket; +/// @brief Entry method to send interaction lists gathered by the +/// treewalk to the GPU for calculation. This has a busy wait on +/// bGPUBufferFilled because the placement of buckets on the GPU needs +/// to be known before the interactions are sent. +/// @param ptr_lc a pointer to the ListCompute object containing the +/// interactions to be calculated. This is declared as an opaque +/// pointer so it can be passed via a Charm++ entry method. +/// @param ptr_ds a pointer to the DoubleWalkState object associated +/// with the walk. This is declared as an opaque pointer so it can +/// be passed via a Charm++ entry method. + +void TreePiece::flushInteractionsToGpu(intptr_t ptr_lc, + intptr_t ptr_ds) { + if(bGPUBufferFilled == 0) { // Need to wait for the buffer + // to be filled before sending + // node interaction lists. + thisProxy[thisIndex].flushInteractionsToGpu(ptr_lc, ptr_ds); + return; + } - int numBucketsDone = data->numBucketsPlusOne-1; - -#if COSMO_PRINT_BK > 1 - CkPrintf("[%d] CUDACALLBACK(node/part: %d, remote: %d, resume: %d), numBucketsDone: %d\n", - tp->getIndex(), - data->node, data->remote, state->resume, - numBucketsDone); -#endif - // bucket computations finished - // -#ifdef CHANGA_REFACTOR_MEMCHECK - CkPrintf("memcheck before cudaCallback\n"); - CmiMemoryCheck(); -#endif - for(int i = 0; i < numBucketsDone; i++){ - bucket = affectedBuckets[i]; - state->counterArrays[0][bucket]--; -#if COSMO_PRINT_BK > 1 - CkPrintf("[%d] bucket %d numAddReq: %d,%d\n", tp->getIndex(), bucket, tp->getSRemoteGravityState()->counterArrays[0][bucket], tp->getSLocalGravityState()->counterArrays[0][bucket]); -#endif - //CkPrintf("[%d] bucket %d numAddReq: %d\n", tp->getIndex(), bucket, state->counterArrays[0][bucket]); - tp->finishBucket(bucket); - } + ListCompute *lc = (ListCompute *)ptr_lc; + DoubleWalkState *ds = (DoubleWalkState *)ptr_ds; -#ifdef CHANGA_REFACTOR_MEMCHECK - CkPrintf("memcheck before cudaCallback freeing\n"); - CmiMemoryCheck(); -#endif - - // free data structures - if(numBucketsDone > 0){ - delete [] data->affectedBuckets; - } - freePinnedHostMemory(data->list); - freePinnedHostMemory(data->bucketMarkers); - freePinnedHostMemory(data->bucketStarts); - freePinnedHostMemory(data->bucketSizes); - if(data->missedNodes) { - freePinnedHostMemory(data->missedNodes); - } - if(data->missedParts) { - freePinnedHostMemory(data->missedParts); - } - - delete ((CkCallback *)data->cb); - delete data; -#ifdef CHANGA_REFACTOR_MEMCHECK - CkPrintf("memcheck after cudaCallback freeing\n"); - CmiMemoryCheck(); -#endif + if(ds->nodeLists.totalNumInteractions > 0){ + lc->sendNodeInteractionsToGpu(ds, this); + lc->resetCudaNodeState(ds); + } + if(ds->particleLists.totalNumInteractions > 0){ + lc->sendPartInteractionsToGpu(ds, this); + lc->resetCudaPartState(ds); + } } void ListCompute::sendNodeInteractionsToGpu(DoubleWalkState *state, @@ -2078,12 +1830,6 @@ void ListCompute::sendNodeInteractionsToGpu(DoubleWalkState *state, data->remote = (getOptType() == Remote); data->missedNodes = NULL; data->missedParts = NULL; - - data->d_localMoments = tp->d_localMoments; - data->d_localParts = tp->d_localParts; - data->d_localVars = tp->d_localVars; - data->d_remoteMoments = tp->d_remoteMoments; - data->stream = tp->stream; #ifdef CUDA_PRINT_TRANSFERRED_INTERACTIONS CkPrintf("*************\n"); @@ -2119,42 +1865,36 @@ void ListCompute::sendNodeInteractionsToGpu(DoubleWalkState *state, #endif OptType type = getOptType(); - data->cb = new CkCallback(cudaCallback, data); - - if(data->numBucketsPlusOne == 1) { // Nothing for the GPU to do - cudaCallback(data, NULL); // Clean up - return; - } if(type == Local){ #ifdef HAPI_TRACE tp->localNodeInteractions += state->nodeLists.totalNumInteractions; #endif - TreePieceCellListDataTransferLocal(data); + peNodeLocalListProxy.ckLocalBranch()->sendList(tp, data); } else if(type == Remote && !state->resume){ #ifdef HAPI_TRACE tp->remoteNodeInteractions += state->nodeLists.totalNumInteractions; #endif - TreePieceCellListDataTransferRemote(data); + peNodeRemoteListProxy.ckLocalBranch()->sendList(tp, data); } else if(type == Remote && state->resume){ CudaMultipoleMoments *missedNodes = state->nodes->getVec(); size_t len = sizeof(CudaMultipoleMoments)*state->nodes->length(); CkAssert(missedNodes); - allocatePinnedHostMemory(&data->missedNodes, len); - memcpy(data->missedNodes, missedNodes, len); + data->missedNodes = state->nodes->getVec(); data->sMissed = len; #ifdef HAPI_TRACE tp->remoteResumeNodeInteractions += state->nodeLists.totalNumInteractions; #endif - TreePieceCellListDataTransferRemoteResume(data); + peNodeRemoteResumeListProxy.ckLocalBranch()->sendList(tp, data); } else {CkAssert(0);} #ifdef CHANGA_REFACTOR_MEMCHECK CkPrintf("memcheck after sendNodeInteractionsToGpu\n"); CmiMemoryCheck(); #endif + delete data; } void ListCompute::sendPartInteractionsToGpu(DoubleWalkState *state, @@ -2175,11 +1915,6 @@ void ListCompute::sendPartInteractionsToGpu(DoubleWalkState *state, data->missedNodes = NULL; data->missedParts = NULL; - data->d_localParts = tp->d_localParts; - data->d_localVars = tp->d_localVars; - data->d_remoteParts = tp->d_remoteParts; - data->stream = tp->stream; - #ifdef CUDA_PRINT_TRANSFERRED_INTERACTIONS CkPrintf("*************\n"); CkPrintf("[%d] partLists: ", thisIndex); @@ -2213,50 +1948,36 @@ void ListCompute::sendPartInteractionsToGpu(DoubleWalkState *state, #endif OptType type = getOptType(); - data->cb = new CkCallback(cudaCallback, data); - - if(data->numBucketsPlusOne == 1) { // Nothing for the GPU to do - cudaCallback(data, NULL); // Clean up - return; - } if(type == Local){ #ifdef HAPI_TRACE tp->localPartInteractions += state->particleLists.totalNumInteractions; #endif - if(tp->largePhase()){ - TreePiecePartListDataTransferLocal(data); - } - else{ - CompactPartData *parts = state->particles->getVec(); - int leng = state->particles->length(); - TreePiecePartListDataTransferLocalSmallPhase(data, parts, leng); - tp->clearMarkedBuckets(state->markedBuckets); - } + pePartLocalListProxy.ckLocalBranch()->sendList(tp, data); } else if(type == Remote && !state->resume){ #ifdef HAPI_TRACE tp->remotePartInteractions += state->particleLists.totalNumInteractions; #endif - TreePiecePartListDataTransferRemote(data); + pePartRemoteListProxy.ckLocalBranch()->sendList(tp, data); } else if(type == Remote && state->resume){ CompactPartData *missedParts = state->particles->getVec(); size_t len = sizeof(CompactPartData)*state->particles->length(); CkAssert(missedParts); - allocatePinnedHostMemory(&data->missedParts, len); - memcpy(data->missedParts, missedParts, len); + data->missedParts = state->particles->getVec(); data->sMissed = len; #ifdef HAPI_TRACE tp->remoteResumePartInteractions += state->particleLists.totalNumInteractions; #endif - TreePiecePartListDataTransferRemoteResume(data); + pePartRemoteResumeListProxy.ckLocalBranch()->sendList(tp, data); } else {CkAssert(0);} #ifdef CHANGA_REFACTOR_MEMCHECK CkPrintf("memcheck after sendPartInteractionsToGpu\n"); CmiMemoryCheck(); #endif + delete data; } #endif diff --git a/Compute.h b/Compute.h index 295be5d4..a11d4499 100644 --- a/Compute.h +++ b/Compute.h @@ -185,9 +185,6 @@ class ListCompute : public Compute{ void enableCpu() {bUseCpu = 1;} #ifdef CUDA -#ifdef GPU_LOCAL_TREE_WALK - void sendLocalTreeWalkTriggerToGpu(State *state, TreePiece *tp, int activeRung, int startBucket, int endBucket); -#endif //GPU_LOCAL_TREE_WALK void sendNodeInteractionsToGpu(DoubleWalkState *state, TreePiece *tp); void sendPartInteractionsToGpu(DoubleWalkState *state, TreePiece *tp); #endif diff --git a/CudaFunctions.h b/CudaFunctions.h index f9d461a5..c3ad5f4f 100644 --- a/CudaFunctions.h +++ b/CudaFunctions.h @@ -4,8 +4,8 @@ #ifdef CUDA #include "hapi.h" #include "HostCUDA.h" -void TreePieceDataTransferBasic(CudaRequest *data, CudaDevPtr *ptr); -void TreePieceDataTransferBasicCleanup(CudaDevPtr *ptr); +void DataTransferBasic(CudaRequest *data, CudaDevPtr *ptr, const char* functionTag); +void DataTransferBasicCleanup(CudaDevPtr *ptr, cudaStream_t stream, const char* functionTag); #ifdef GPU_LOCAL_TREE_WALK __global__ void gpuLocalTreeWalk( @@ -36,7 +36,6 @@ __global__ void nodeGravityComputation( __global__ void ZeroVars(VariablePartData *particleVars, int nVars); -#ifdef CUDA_2D_TB_KERNEL __global__ void particleGravityComputation( CompactPartData *targetCores, VariablePartData *targetVars, @@ -46,18 +45,6 @@ __global__ void particleGravityComputation( int *bucketStarts, int *bucketSizes, cudatype fperiod); -#else -__global__ void particleGravityComputation( - CompactPartData *targetCores, - VariablePartData *targetVars, - CompactPartData *sourceCores, - ILPart *ils, - int *ilmarks, - int *bucketStarts, - int *bucketSizes, - cudatype fperiod); -#endif - #endif #endif diff --git a/DataManager.cpp b/DataManager.cpp index 52110130..adf322a4 100644 --- a/DataManager.cpp +++ b/DataManager.cpp @@ -12,6 +12,8 @@ #include "hapi.h" #include "cuda_typedef.h" #include "SFC.h" +#include "MemoryPool.h" +#include "PEList.h" #endif #include "Compute.h" @@ -33,14 +35,36 @@ void DataManager::init() { oldNumChunks = 0; chunkRoots = NULL; #ifdef CUDA + // Initialize host memory pool (per OS process) + hostPoolInit(); + + // Initialize GPU memory pool (node-level synchronization inside gpuPoolInit) + gpuPoolInit(); + treePiecesDone = 0; treePiecesDonePrefetch = 0; - treePiecesDoneLocalComputation = 0; - treePiecesDoneRemoteChunkComputation = 0; - treePiecesWantParticlesBack = 0; + bEwaldDone = false; + PEsWantParticlesBack = 0; treePiecesParticlesUpdated = 0; - gpuFree = true; - + memLog = new MemLog(); + lockMemLog = CmiCreateLock(); + bGpuMemLogger = 0; // Default disabled + cpuMemLog = new MemLog(); + lockCpuMemLog = CmiCreateLock(); + bCpuMemLogger = 0; // Default disabled + dHostPoolTargetCapacityGB = 5.0; // Default target capacity + nHostPoolMinCapacityPerBucketMB = 50; // Default minimum per bucket + bHostPoolDebug = 0; // Default disabled + + d_localMoments = nullptr; + d_localParts = nullptr; + d_localVars = nullptr; + d_remoteMoments = nullptr; + d_remoteParts = nullptr; + ewtGPU = nullptr; + cachedData = nullptr; + bLocalDataTransferred.store(false); + bRemoteDataTransferred.store(false); #endif Cool = CoolInit(); LWData = LymanWernerTableInit(); @@ -51,14 +75,9 @@ void DataManager::init() { } #ifdef CUDA -/// @brief Initialize CUDA streams -/// @param _numStreams Total number of streams to create -void DataManager::createStreams(int _numStreams, const CkCallback& cb) { - numStreams = _numStreams; - streams = new cudaStream_t[numStreams]; - for (int i = 0; i < numStreams; i++) { - hapiCheck(cudaStreamCreate(&streams[i])); - } +/// @brief Initialize CUDA stream for this DataManager +void DataManager::createStream(const CkCallback& cb) { + cudaStreamCreate(&stream); contribute(cb); } #endif @@ -152,8 +171,7 @@ void DataManager::notifyPresence(Tree::GenericTreeNode *root, TreePiece *tp) { CmiLock(__nodelock); registeredTreePieces.push_back(TreePieceDescriptor(tp, root)); #ifdef CUDA - //gpuFree = true; - //registeredTreePieceIndices.push_back(index); + registeredPEs.insert(tp->getParentPE()); #if COSMO_PRINT_BK > 1 CkPrintf("(%d) notifyPresence called by %d, length: %d\n", CkMyPe(), tp->getIndex(), registeredTreePieces.length()); #endif @@ -164,10 +182,12 @@ void DataManager::notifyPresence(Tree::GenericTreeNode *root, TreePiece *tp) { /// \brief Clear registeredTreePieces on this node. void DataManager::clearRegisteredPieces(const CkCallback& cb) { registeredTreePieces.removeAll(); +#ifdef CUDA + registeredPEs.clear(); +#endif contribute(cb); } - /// \brief Build a local tree inside the node. /// /// This will be an exact superset of all the trees in this @@ -249,10 +269,18 @@ void DataManager::combineLocalTrees(CkReductionMsg *msg) { } } + TreePiece *tp = registeredTreePieces[0].treePiece; + // Get periodic parameters into DataManager + // and build Ewald Table for use by all TPs on this node + if(tp->bEwald) { + fPeriod = tp->fPeriod; + dEwhCut = tp->dEwhCut; + fEwCut = tp->fEwCut; + nReplicas = tp->nReplicas; + EwaldInit(); + } + } -#ifdef CUDA - gpuFree = true; -#endif contribute(*(CkCallback*)msg->getData()); delete msg; } @@ -450,6 +478,13 @@ void DataManager::resetReadOnly(Parameters param, const CkCallback &cb) bUseCkLoopPar = param.bUseCkLoopPar; #else bUseCkLoopPar = 0; +#endif +#ifdef CUDA + bGpuMemLogger = param.bGpuMemLogger; + bCpuMemLogger = param.bCpuMemLogger; + dHostPoolTargetCapacityGB = param.dHostPoolTargetCapacityGB; + nHostPoolMinCapacityPerBucketMB = param.nHostPoolMinCapacityPerBucketMB; + bHostPoolDebug = param.bHostPoolDebug; #endif contribute(cb); // parameter structure requires some cleanup @@ -486,33 +521,246 @@ void DataManager::serializeLocalTree(){ CmiUnlock(__nodelock); } +/// @brief Get the data produced by DataManager::EwaldInit and launch +/// the Ewald kernel on the GPU +void DataManager::startEwaldGPU() { + // Early exit when no Ewald work to do: no particles on this node, or GPU + // buffers were never allocated (e.g. no-local-parts run). + // + // We must still call cudaFinishAllBuckets(1) here. Normally that happens in + // finishEwaldGPU() (the callback passed to DataManagerEwald). When we skip + // the kernel, finishEwaldGPU never runs, so we replicate its behavior here + // to avoid deadlock: TPs are waiting for their buckets to be marked finished. + if (savedNumTotalParticles <= 0 || d_localParts == nullptr || d_localVars == nullptr) { + for (int i = 0; i < registeredTreePieces.length(); i++) { + if(registeredTreePieces[i].treePiece->getNumActiveParticles() > 0) { + int in = registeredTreePieces[i].treePiece->getIndex(); + treePieces[in].cudaFinishAllBuckets(1); + } + } + return; + } + + int nEwhLoop = ewt.size(); + if (nEwhLoop > NEWH) { + CkAbort("DataManager::startEwaldGPU: nEwhLoop (%d) exceeds NEWH (%d); increase NEWH in EwaldCUDA.h", + nEwhLoop, NEWH); + } + + if(ewtGPU == nullptr) { +#ifdef PINNED_HOST_MEMORY + const char* funcTag = "DataManager::startEwaldGPU"; + hostMalloc(&ewtGPU, sizeof(EwtData)*NEWH, funcTag); + hostMalloc(&cachedData, sizeof(EwaldReadOnlyData), funcTag); +#else + ewtGPU = (EwtData *) malloc(sizeof(EwtData)*NEWH); + cachedData = (EwaldReadOnlyData *) malloc(sizeof(EwaldReadOnlyData)); +#endif + } + + MultipoleMoments *mm = &root->moments; + for (int i=0; imomcRoot.m = (cudatype) momcRoot.m; + cachedData->momcRoot.xx = (cudatype) momcRoot.xx; + cachedData->momcRoot.yy = (cudatype) momcRoot.yy; + cachedData->momcRoot.xy = (cudatype) momcRoot.xy; + cachedData->momcRoot.xz = (cudatype) momcRoot.xz; + cachedData->momcRoot.yz = (cudatype) momcRoot.yz; + cachedData->momcRoot.xxx = (cudatype) momcRoot.xxx; + cachedData->momcRoot.xyy = (cudatype) momcRoot.xyy; + cachedData->momcRoot.xxy = (cudatype) momcRoot.xxy; + cachedData->momcRoot.yyy = (cudatype) momcRoot.yyy; + cachedData->momcRoot.xxz = (cudatype) momcRoot.xxz; + cachedData->momcRoot.yyz = (cudatype) momcRoot.yyz; + cachedData->momcRoot.xyz = (cudatype) momcRoot.xyz; + cachedData->momcRoot.xxxx = (cudatype) momcRoot.xxxx; + cachedData->momcRoot.xyyy = (cudatype) momcRoot.xyyy; + cachedData->momcRoot.xxxy = (cudatype) momcRoot.xxxy; + cachedData->momcRoot.yyyy = (cudatype) momcRoot.yyyy; + cachedData->momcRoot.xxxz = (cudatype) momcRoot.xxxz; + cachedData->momcRoot.yyyz = (cudatype) momcRoot.yyyz; + cachedData->momcRoot.xxyy = (cudatype) momcRoot.xxyy; + cachedData->momcRoot.xxyz = (cudatype) momcRoot.xxyz; + cachedData->momcRoot.xyyz = (cudatype) momcRoot.xyyz; + cachedData->momcRoot.zz = (cudatype) momcRoot.zz; + cachedData->momcRoot.xzz = (cudatype) momcRoot.xzz; + cachedData->momcRoot.yzz = (cudatype) momcRoot.yzz; + cachedData->momcRoot.zzz = (cudatype) momcRoot.zzz; + cachedData->momcRoot.xxzz = (cudatype) momcRoot.xxzz; + cachedData->momcRoot.xyzz = (cudatype) momcRoot.xyzz; + cachedData->momcRoot.xzzz = (cudatype) momcRoot.xzzz; + cachedData->momcRoot.yyzz = (cudatype) momcRoot.yyzz; + cachedData->momcRoot.yzzz = (cudatype) momcRoot.yzzz; + cachedData->momcRoot.zzzz = (cudatype) momcRoot.zzzz; +#else + cachedData->mm.xx = (cudatype) mm->xx; + cachedData->mm.xy = (cudatype) mm->xy; + cachedData->mm.xz = (cudatype) mm->xz; + cachedData->mm.yy = (cudatype) mm->yy; + cachedData->mm.yz = (cudatype) mm->yz; + cachedData->mm.zz = (cudatype) mm->zz; +#endif + cudatype L = fPeriod.x; + cudatype alpha = 2.0f/L; + + cachedData->mm.totalMass = (cudatype) mm->totalMass; + cachedData->mm.cmx = (cudatype) mm->cm.x; + cachedData->mm.cmy = (cudatype) mm->cm.y; + cachedData->mm.cmz = (cudatype) mm->cm.z; + cachedData->fEwCut = (cudatype) fEwCut; + cachedData->nReps = nReplicas; + cachedData->nEwReps = (int) ceil(fEwCut); + cachedData->nEwhLoop = nEwhLoop; + cachedData->L = L; + cachedData->alpha = alpha; + cachedData->alpha2 = alpha*alpha; + cachedData->k1 = (cudatype) M_PI/(alpha*alpha*L*L*L); + cachedData->ka = (cudatype) 2.0*alpha/sqrt(M_PI); + cachedData->fEwCut2 = (cudatype) fEwCut*fEwCut*L*L; + + /* + Break between Taylor expansion for small r and multipole expansion. + This value is for double precision. + roData->fInner2 = (cudatype) 1.2e-3*L*L; + The following is for single precision. The CUDA version currently uses + erff() and erfcf(). If these ever get changed then the following line + needs to be changed accordingly. + */ + cachedData->fInner2 = (cudatype) 1.1e-2*L*L; + + ewaldCallback + = new CkCallback(CkIndex_DataManager::finishEwaldGPU(), CkMyNode(), dMProxy); + + DataManagerEwald(d_localParts, d_localVars, ewtGPU, cachedData, + savedNumTotalParticles, stream, ewaldCallback); +} + +/// @brief Callback from Ewald kernel launch on GPU +void DataManager::finishEwaldGPU() { + delete ewaldCallback; + + for(int i = 0; i < registeredTreePieces.length(); i++){ + if(registeredTreePieces[i].treePiece->getNumActiveParticles() > 0) { + int in = registeredTreePieces[i].treePiece->getIndex(); + treePieces[in].cudaFinishAllBuckets(1); + } + } + transferParticleVarsBack(true); +} + + +/// @brief Callback from local tree walk on GPU +/// Call finishBucket for all buckets and TreePieces +/// Start Ewald calculation if enabled +void DataManager::finishLocalWalk() { + delete localWalkCallback; + delete lwReq; + +#ifdef GPU_LOCAL_TREE_WALK +#ifdef PINNED_HOST_MEMORY + // Direct free for large local tree buffers (bypasses pool but logs analytics) + const char* funcTag = "DataManager::finishLocalWalk"; + hostFree(bufLocalMoments, funcTag); + hostFree(bufLocalParts, funcTag); +#else + free(bufLocalMoments); + free(bufLocalParts); +#endif + + for(int i = 0; i < registeredTreePieces.length(); i++){ + if(registeredTreePieces[i].treePiece->getNumActiveParticles() > 0) { + int in = registeredTreePieces[i].treePiece->getIndex(); + treePieces[in].cudaFinishAllBuckets(0); + } + } +#endif + + if (registeredTreePieces[0].treePiece->bEwald) { + startEwaldGPU(); + } + else { + transferParticleVarsBack(true); + } +} + /// @brief Callback from local data transfer to GPU -/// Indicate the transfer is done, and start the local gravity walks -/// on the treepieces on this node. +/// Indicate the transfer is done, and start the local gravity walk +/// in one big kernel launch void DataManager::startLocalWalk() { delete localTransferCallback; + bLocalDataTransferred.store(true); + + // Delayed local kernels can now launch (remote may have completed first). + int pe; + int firstPE = CkNodeFirst(CkMyNode()); + int nPEs = CkNodeSize(CkMyNode()); + for (int i = 0; i < numPEListProxies; i++) { + for (int j = 0; j < nPEs; j++) { + pe = firstPE + j; + (*(PEListProxies[i]))[pe].tryLaunchDelayedKernel(); + } + } + // We arent calculating local gravity on the CPU, but bookkeeping + // still needs to be handled for(int i = 0; i < registeredTreePieces.length(); i++){ - if(verbosity > 1) CkPrintf("[%d] GravityLocal %d\n", CkMyPe(), i); - int in = registeredTreePieces[i].treePiece->getIndex(); - treePieces[in].commenceCalculateGravityLocal((intptr_t)d_localMoments, - (intptr_t)d_localParts, - (intptr_t)d_localVars, - (intptr_t)streams, numStreams, - sMoments, sCompactParts, sVarParts); - if(registeredTreePieces[0].treePiece->bEwald) { - EwaldMsg *msg = new (8*sizeof(int)) EwaldMsg; - msg->fromInit = false; - // Make priority lower than gravity or smooth. - *((int *)CkPriorityPtr(msg)) = 3*numTreePieces + in + 1; - CkSetQueueing(msg,CK_QUEUEING_IFIFO); - treePieces[in].calculateEwald(msg); - } + /// Pieces with no active particles have been short-cut. + if(registeredTreePieces[i].treePiece->getNumActiveParticles() > 0) { + int in = registeredTreePieces[i].treePiece->getIndex(); + treePieces[in].commenceCalculateGravityLocal(); + } } - freePinnedHostMemory(bufLocalMoments); - freePinnedHostMemory(bufLocalParts); - freePinnedHostMemory(bufLocalVars); +#ifdef GPU_LOCAL_TREE_WALK + localWalkCallback + = new CkCallback(CkIndex_DataManager::finishLocalWalk(), CkMyNode(), dMProxy); + + lwReq = new CudaRequest; + + lwReq->d_localMoments = d_localMoments; + lwReq->d_localParts = d_localParts; + lwReq->d_localVars = d_localVars; + lwReq->sMoments = sMoments; + lwReq->sCompactParts = sCompactParts; + lwReq->sVarParts = sVarParts; + lwReq->stream = stream; + + lwReq->numBucketsPlusOne = 0; + lwReq->affectedBuckets = 0; + lwReq->tp = this; + lwReq->state = NULL; + lwReq->node = true; + lwReq->remote = false; + lwReq->firstParticle = 0; + lwReq->lastParticle = savedNumTotalParticles-1; + + lwReq->rootIdx = 0; + lwReq->theta = theta; + lwReq->thetaMono = thetaMono; + lwReq->nReplicas = registeredTreePieces[0].treePiece->nReplicas; + lwReq->fperiod = registeredTreePieces[0].treePiece->fPeriod.x; + lwReq->fperiodY = registeredTreePieces[0].treePiece->fPeriod.y; + lwReq->fperiodZ = registeredTreePieces[0].treePiece->fPeriod.z; + lwReq->cb = localWalkCallback; + + lwReq->list = NULL; + lwReq->bucketMarkers = NULL; + lwReq->bucketStarts = NULL; + lwReq->bucketSizes = NULL; + lwReq->numInteractions = 0; + + DataManagerLocalTreeWalk(lwReq); +#else + finishLocalWalk(); +#endif } /// @brief Callback from remote data transfer to GPU. @@ -520,28 +768,23 @@ void DataManager::startLocalWalk() { /// remote walk. void DataManager::resumeRemoteChunk() { if(verbosity > 1) CkPrintf("[%d] resumeRemoteChunk registered: %lu\n", CkMyPe(), registeredTreePieces.length()); - int chunk = 0; - chunk = currentChunkBuffers->chunk; delete currentChunkBuffers->moments; delete currentChunkBuffers->particles; delete currentChunkBuffers->cb; delete currentChunkBuffers; - - if(bufRemoteMoments != NULL) - freePinnedHostMemory(bufRemoteMoments); - if(bufRemoteParts != NULL) - freePinnedHostMemory(bufRemoteParts); - - // resume each treepiece's startRemoteChunk, now that the nodes - // are properly labeled and the particles accounted for - for(int i = 0; i < registeredTreePieces.length(); i++){ - if(verbosity > 1) CkPrintf("[%d] resumeRemoteChunk %d\n", CkMyPe(), i); - int in = registeredTreePieces[i].treePiece->getIndex(); -#if COSMO_PRINT_BK > 1 - CkPrintf("(%d) dm->%d\n", CkMyPe(), in); -#endif - treePieces[in].continueStartRemoteChunk(chunk, (intptr_t)d_remoteMoments, (intptr_t)d_remoteParts); - } + bRemoteDataTransferred.store(true); + + // Check and see if the remote walks already finished and are waiting + // to launch their GPU kernels + int pe; + int firstPE = CkNodeFirst(CkMyNode()); + int nPEs = CkNodeSize(CkMyNode()); + for (int i = 0; i < numPEListProxies; i++) { + for (int j = 0; j < nPEs; j++) { + pe = firstPE + j; + (*(PEListProxies[i]))[pe].tryLaunchDelayedKernel(); + } + } } /// @brief record when all TreePieces have finished their prefetch. @@ -553,6 +796,7 @@ void DataManager::donePrefetch(int chunk){ treePiecesDonePrefetch++; if(treePiecesDonePrefetch == registeredTreePieces.length()){ treePiecesDonePrefetch = 0; + #ifdef HAPI_TRACE double starttime = CmiWallTimer(); #endif @@ -560,45 +804,60 @@ void DataManager::donePrefetch(int chunk){ #ifdef HAPI_TRACE traceUserBracketEvent(CUDA_SER_TREE, starttime, CmiWallTimer()); #endif - PendingBuffers *buffers = currentChunkBuffers; - if(gpuFree){ - gpuFree = false; - lastChunkMoments = buffers->moments->length(); - lastChunkParticles = buffers->particles->length(); - - CkCallback *remoteChunkTransferCallback - = new CkCallback(CkIndex_DataManager::resumeRemoteChunk(), CkMyNode(), - dMProxy); - buffers->cb = remoteChunkTransferCallback; - // XXX copies can be saved here. - size_t sRemMoments = lastChunkMoments*sizeof(CudaMultipoleMoments); - if(sRemMoments > 0) { - allocatePinnedHostMemory((void **)&bufRemoteMoments, sRemMoments); - memcpy(bufRemoteMoments, buffers->moments->getVec(), sRemMoments); - } - else - bufRemoteMoments = NULL; - size_t sRemParts = lastChunkParticles*sizeof(CompactPartData); - if(sRemParts > 0) { - allocatePinnedHostMemory((void **)&bufRemoteParts, sRemParts); - memcpy(bufRemoteParts, buffers->particles->getVec(), sRemParts); - } - else - bufRemoteParts = NULL; - - // Transfer moments and particle cores to gpu - DataManagerTransferRemoteChunk(bufRemoteMoments, sRemMoments, - bufRemoteParts, sRemParts, - (void **)&d_remoteMoments, (void **)&d_remoteParts, - streams[0], - remoteChunkTransferCallback); + // Commence the remote walk + for(int i = 0; i < registeredTreePieces.length(); i++){ + if(verbosity > 1) CkPrintf("[%d] continueStartRemoteChunk %d\n", CkMyNode(), i); + if(registeredTreePieces[i].treePiece->getNumActiveParticles() > 0) { + int in = registeredTreePieces[i].treePiece->getIndex(); +#if COSMO_PRINT_BK > 1 + CkPrintf("(%d) dm->%d\n", CkMyNode(), in); +#endif + treePieces[in].continueStartRemoteChunk(0); // 0 is the chunk number; + // assuming only one remote chunk + } } - else{ - // enqueue pendingbuffers - pendingChunkTransferQ.enq(buffers); - } - + + PendingBuffers *buffers = currentChunkBuffers; + lastChunkMoments = buffers->moments->length(); + lastChunkParticles = buffers->particles->length(); + + CkCallback *remoteChunkTransferCallback + = new CkCallback(CkIndex_DataManager::resumeRemoteChunk(), CkMyNode(), + dMProxy); + buffers->cb = remoteChunkTransferCallback; + // XXX copies can be saved here. + size_t sRemMoments = lastChunkMoments*sizeof(CudaMultipoleMoments); + if(sRemMoments > 0) { +#ifdef PINNED_HOST_MEMORY + const char* funcTag = "DataManager::transferPrefetch"; + hostMalloc(&bufRemoteMoments, sRemMoments, funcTag); +#else + bufRemoteMoments = (CudaMultipoleMoments *) malloc(sRemMoments); +#endif + memcpy(bufRemoteMoments, buffers->moments->getVec(), sRemMoments); + } + else + bufRemoteMoments = NULL; + size_t sRemParts = lastChunkParticles*sizeof(CompactPartData); + if(sRemParts > 0) { +#ifdef PINNED_HOST_MEMORY + const char* funcTag = "DataManager::transferPrefetch"; + hostMalloc(&bufRemoteParts, sRemParts, funcTag); +#else + bufRemoteParts = (CompactPartData *) malloc(sRemParts); +#endif + memcpy(bufRemoteParts, buffers->particles->getVec(), sRemParts); + } + else + bufRemoteParts = NULL; + + // Transfer moments and particle cores to gpu + DataManagerTransferRemoteChunk(bufRemoteMoments, sRemMoments, + bufRemoteParts, sRemParts, + (void **)&d_remoteMoments, (void **)&d_remoteParts, + stream, + remoteChunkTransferCallback); } CmiUnlock(__nodelock); } @@ -633,6 +892,8 @@ static inline void addTreeNodeToList(GenericTreeNode *nd, CudaMultipoleMoments cmm(nd->moments); cmm.lesser_corner = nd->boundingBox.lesser_corner; cmm.greater_corner = nd->boundingBox.greater_corner; + cmm.children[0] = -1; + cmm.children[0] = -1; list.push_back(cmm); index++; } @@ -648,7 +909,6 @@ PendingBuffers *DataManager::serializeRemoteChunk(GenericTreeNode *node){ int numTreePieces = registeredTreePieces.length(); int numNodes = 0; int numParticles = 0; - int totalNumBuckets = 0; cacheType *wholeNodeCache = cacheNode.ckLocalBranch()->getCache(); cacheType *ctNode = &wholeNodeCache[chunk]; @@ -844,25 +1104,54 @@ void DataManager::serializeLocal(GenericTreeNode *nodeRoot){ #endif size_t sLocalParts = numParticles*sizeof(CompactPartData); size_t sLocalMoments = localMoments.length()*sizeof(CudaMultipoleMoments); - allocatePinnedHostMemory((void **)&bufLocalParts, sLocalParts); - allocatePinnedHostMemory((void **)&bufLocalMoments, sLocalMoments); +#ifdef PINNED_HOST_MEMORY + // Bypass pool for large local tree buffers (but log analytics) + const char* funcTag = "DataManager::sendLocalData"; + hostMalloc(&bufLocalParts, sLocalParts, funcTag); + hostMalloc(&bufLocalMoments, sLocalMoments, funcTag); +#else + bufLocalParts = (CompactPartData *) malloc(sLocalParts); + bufLocalMoments = (CudaMultipoleMoments *) malloc(sLocalMoments); +#endif int pTPindex = 0; treePiecesBufferFilled = 0; for(int i = 0; i < numTreePieces; i++){ - treePieces[registeredTreePieces[i].treePiece->getIndex()].fillGPUBuffer((intptr_t) bufLocalParts, - (intptr_t) bufLocalMoments, (intptr_t) localMoments.getVec(), pTPindex, - numParticles, (intptr_t) nodeRoot); + // Give fillGPUBuffer priority over SPH + // We want the GPU kernels to start ASAP to minimize idle time on the host + fillGPUMsg *msg = new (8*sizeof(int)) fillGPUMsg; + msg->partIndex = pTPindex; + msg->nParts = numParticles; + *((int *)CkPriorityPtr(msg)) = -100000000 - i; + CkSetQueueing(msg,CK_QUEUEING_IFIFO); + treePieces[registeredTreePieces[i].treePiece->getIndex()].fillGPUBuffer(msg); pTPindex += registeredTreePieces[i].treePiece->getDMNumParticles(); } } +#ifdef GPU_LOCAL_TREE_WALK_DEBUG +static std::unordered_set myLoopCheck; +int maxDepth(int iNode, CkVec& localMoments, int inDepth = 0) { + if(iNode == -1) return 0; + if(inDepth > 64) + CkPrintf("%d: indepth %d\n", CkMyNode(), inDepth); + if(myLoopCheck.find(iNode) != myLoopCheck.end()) + CkPrintf("%d: duplicate node %d depth %d type %d\n", CkMyNode(), iNode, inDepth, + localMoments[iNode].type); + CkAssert(myLoopCheck.find(iNode) == myLoopCheck.end()); + myLoopCheck.insert(iNode); + CkAssert((iNode >= -1) && (iNode < (int) localMoments.length())); + int depth = 1 + max(maxDepth(localMoments[iNode].children[0], localMoments, inDepth+1), + maxDepth(localMoments[iNode].children[1], localMoments, inDepth+1)); + return depth; +} +#endif + /// /// @brief After all pieces have filled the buffer, initiate the transfer. /// @param numParticles total number of particles on this node -/// @param node root of tree /// -void DataManager::transferLocalToGPU(int numParticles, GenericTreeNode *node) +void DataManager::transferLocalToGPU(int numParticles) { CmiLock(__nodelock); treePiecesBufferFilled++; @@ -877,7 +1166,39 @@ void DataManager::transferLocalToGPU(int numParticles, GenericTreeNode *node) double starttime = CmiWallTimer(); #ifdef GPU_LOCAL_TREE_WALK - transformLocalTreeRecursive(node, localMoments); + transformLocalTreeRecursive(root, localMoments); + // Sanitiy check on tree +#ifdef GPU_LOCAL_TREE_WALK_DEBUG + myLoopCheck.clear(); + for(int i = 0; i < localMoments.length(); i++) { + CkAssert(localMoments[i].type != NonLocal); + CkAssert(localMoments[i].type != Empty); + CkAssert(localMoments[i].type != CachedEmpty); + for (int j = 0; j < 2; j ++) { + int child = localMoments[i].children[j]; + if((child < -1) || (child >= (int) localMoments.length())) + CkPrintf("%d %d: Bad moment %d of %ld: type %d child %d %d %d\n", + CkMyNode(), CkMyPe(), + i, localMoments.length(), + localMoments[i].type, + child, + localMoments[i].children[0], + localMoments[i].children[1]); + if(child >= 0) { + if(myLoopCheck.find(child) != myLoopCheck.end()) + CkPrintf("%d: duplicate child node %d type %d child %d index %d\n", + CkMyNode(), i, localMoments[i].type, j, child); + CkAssert(myLoopCheck.find(child) == myLoopCheck.end()); + myLoopCheck.insert(child); + } + CkAssert((child >= -1) && (child < (int) localMoments.length())); + } + } + myLoopCheck.clear(); + CkPrintf("%d %d: maxDepth: %d\n", CkMyNode(), CkMyPe(), + maxDepth(0, localMoments)); +#endif + #endif //GPU_LOCAL_TREE_WALK #ifdef HAPI_TRACE traceUserBracketEvent(SER_LOCAL_TRANSFORM, starttime, CmiWallTimer()); @@ -897,13 +1218,11 @@ void DataManager::transferLocalToGPU(int numParticles, GenericTreeNode *node) traceUserBracketEvent(SER_LOCAL_MEMCPY, starttime, CmiWallTimer()); #endif - allocatePinnedHostMemory((void **)&bufLocalVars, sLocalVars); - // Transfer moments and particle cores to gpu DataManagerTransferLocalTree(bufLocalMoments, sLocalMoments, bufLocalParts, - sLocalParts, bufLocalVars, sLocalVars, + sLocalParts, sLocalVars, (void **)&d_localMoments, (void **)&d_localParts, (void **)&d_localVars, - streams[0], numParticles, + stream, numParticles, localTransferCallback); } @@ -913,6 +1232,8 @@ void DataManager::transformLocalTreeRecursive(GenericTreeNode *node, CkVecgetType(); int node_index = node->nodeArrayIndex; + CkAssert(type != NonLocal); + if(type == Empty || type == CachedEmpty){ // skip return; } else if(type == Bucket || type == NonLocalBucket) { @@ -938,14 +1259,17 @@ void DataManager::transformLocalTreeRecursive(GenericTreeNode *node, CkVecnumChildren(); i++){ GenericTreeNode *child = node->getChildren(i); - int child_index = child->nodeArrayIndex; - localMoments[node_index].children[i] = child_index; - transformLocalTreeRecursive(child, localMoments); - - // child_index == -1 can indicate an empty node or a non-local node. - if (child_index != -1 && localMoments[child_index].bucketSize > 0) { - localMoments[node_index].bucketStart = std::min(localMoments[node_index].bucketStart, localMoments[child_index].bucketStart); - localMoments[node_index].bucketSize += localMoments[child_index].bucketSize; + if(child->getType() != NonLocal) { // NonLocal nodes are not in + // the localMoments vector. + int child_index = child->nodeArrayIndex; + localMoments[node_index].children[i] = child_index; + transformLocalTreeRecursive(child, localMoments); + + // child_index == -1 can indicate an empty node or a non-local node. + if (child_index != -1 && localMoments[child_index].bucketSize > 0) { + localMoments[node_index].bucketStart = std::min(localMoments[node_index].bucketStart, localMoments[child_index].bucketStart); + localMoments[node_index].bucketSize += localMoments[child_index].bucketSize; + } } } } @@ -957,17 +1281,24 @@ void updateParticlesCallback(void *, void *); /// @brief Copy particle accelerations back from GPU to host memory and /// deallocate the device memory /// This is triggered when all TreePieces call finishBucket -void DataManager::transferParticleVarsBack(){ +/// @param bFromEwald the checkin is from the finishEwaldGPU OR Ewald +/// was skipped (non-periodic) +void DataManager::transferParticleVarsBack(bool bFromEwald){ UpdateParticlesStruct *data; CmiLock(__nodelock); - treePiecesWantParticlesBack++; - if(treePiecesWantParticlesBack == registeredTreePieces.length()){ - treePiecesWantParticlesBack = 0; + if(bFromEwald) + bEwaldDone = true; + else + PEsWantParticlesBack++; + if(PEsWantParticlesBack == registeredPEs.size()*numPEListProxies && bEwaldDone){ + bEwaldDone = false; + PEsWantParticlesBack = 0; VariablePartData *buf; if(savedNumTotalParticles > 0){ #ifdef PINNED_HOST_MEMORY - allocatePinnedHostMemory((void **)&buf, savedNumTotalParticles*sizeof(VariablePartData)); + const char* funcTag = "DataManager::transferParticleVarsBack"; + hostMalloc(&buf, savedNumTotalParticles*sizeof(VariablePartData), funcTag); #else buf = (VariablePartData *) malloc(savedNumTotalParticles*sizeof(VariablePartData)); #endif @@ -982,22 +1313,12 @@ void DataManager::transferParticleVarsBack(){ data->buf = buf; data->size = savedNumTotalParticles; - if(verbosity > 1) CkPrintf("[%d] transferParticleVarsBack\n", CkMyPe()); + if(verbosity > 1) CkPrintf("[%d] transferParticleVarsBack\n", CkMyNode()); TransferParticleVarsBack(buf, savedNumTotalParticles*sizeof(VariablePartData), d_localVars, - streams[0], + stream, data->cb); - - cudaFree(d_localMoments); - cudaFree(d_localParts); - cudaFree(d_localVars); - cudaFree(d_remoteMoments); - cudaFree(d_remoteParts); - -#ifdef CUDA_PRINT_ERRORS - printf("transferParticleVarsBack: %s\n", cudaGetErrorString( cudaGetLastError() ) ); -#endif } CmiUnlock(__nodelock); } @@ -1005,10 +1326,8 @@ void DataManager::transferParticleVarsBack(){ void DataManager::updateParticles(UpdateParticlesStruct *data){ int partIndex = 0; - VariablePartData *deviceParticles = data->buf; - #ifdef CUDA_PRINT_TRANSFER_BACK_PARTICLES - CkPrintf("(%d) In DM::updateParticles %d tps\n", CkMyPe(), registeredTreePieces.length()); + CkPrintf("(%d) In DM::updateParticles %d tps\n", CkMyNode(), (int) registeredTreePieces.length()); #endif for(int i = 0; i < registeredTreePieces.length(); i++){ @@ -1045,17 +1364,64 @@ void DataManager::updateParticlesFreeMemory(UpdateParticlesStruct *data) if(treePiecesParticlesUpdated == registeredTreePieces.length()){ treePiecesParticlesUpdated = 0; + // Free host buffers for remote chunk data +#ifdef PINNED_HOST_MEMORY + const char* funcTagHost = "DataManager::updateParticlesFreeMemory"; + if(bufRemoteMoments != NULL) + hostFree(bufRemoteMoments, funcTagHost); + if(bufRemoteParts != NULL) + hostFree(bufRemoteParts, funcTagHost); +#else + if(bufRemoteMoments != NULL) + free(bufRemoteMoments); + if(bufRemoteParts != NULL) + free(bufRemoteParts); +#endif + + // Free device memory + const char* funcTag = "DataManager::updateParticlesFreeMemory"; + gpuFree(d_localMoments, stream, funcTag); + gpuFree(d_localParts, stream, funcTag); + gpuFree(d_localVars, stream, funcTag); + gpuPoolFree(d_remoteMoments, stream, funcTag); + gpuPoolFree(d_remoteParts, stream, funcTag); + + bLocalDataTransferred.store(false); + bRemoteDataTransferred.store(false); + // Set device pointers to nullptr + d_localMoments = nullptr; + d_localParts = nullptr; + d_localVars = nullptr; + d_remoteMoments = nullptr; + d_remoteParts = nullptr; + if(data->size > 0){ #ifdef PINNED_HOST_MEMORY - freePinnedHostMemory(data->buf); + const char* funcTag = "DataManager::updateParticlesFreeMemory"; + hostFree(data->buf, funcTag); #else free(data->buf); #endif } delete (data->cb); delete data; + cudaDeviceSynchronize(); + } CmiUnlock(__nodelock); } +// After after each BigStep, trims the host pool to reclaim memory +void DataManager::trimHostPool(double targetCapacityGB, const CkCallback& cb){ + hostPoolTrim(targetCapacityGB, nHostPoolMinCapacityPerBucketMB); + contribute(cb); +} + +// After trimming, refills hot buckets to prepare for next timestep +void DataManager::refillHostPool(const CkCallback& cb){ + hostPoolAdaptiveRefill(); + contribute(cb); +} + + #endif // CUDA diff --git a/DataManager.h b/DataManager.h index 4c2c5fe6..80ddf52c 100644 --- a/DataManager.h +++ b/DataManager.h @@ -9,11 +9,17 @@ #include #include +#include #include #include "GenericTreeNode.h" #include "ParallelGravity.decl.h" #include "lymanwerner.h" +#ifdef CUDA +#include "memlog.h" +#include "MemoryPool.h" +#endif + #if CHARM_VERSION > 60401 && CMK_BALANCED_INJECTION_API #include "ckBIconfig.h" #endif @@ -30,6 +36,14 @@ struct TreePieceDescriptor{ } }; +/* IBM brain damage */ +#undef hz +/// @brief Coefficients for the Fourier space part of the Ewald sum. +typedef struct ewaldTable { + double hx,hy,hz; + double hCfac,hSfac; +} EWT; + #ifdef CUDA struct UpdateParticlesStruct{ @@ -79,18 +93,26 @@ class DataManager : public CBase_DataManager { // holds chare array indices of registered treepieces CkVec registeredTreePieces; #ifdef CUDA - //CkVec registeredTreePieceIndices; /// @brief counter for the number of tree nodes that are /// replicated by TreePieces that share the same address space. int cumNumReplicatedNodes; + /// Counter for TreePieces ready for serializeLocalTree() int treePiecesDone; int savedChunk; + /// Counter for TreePieces which have completed their remote Prefetch. int treePiecesDonePrefetch; - int treePiecesDoneLocalComputation; - // XXX - assumes that only one chunk can be on the gpu - // at a given time - int treePiecesDoneRemoteChunkComputation; - int treePiecesWantParticlesBack; +public: + /// Indicates that the local Tree data is on the GPU + std::atomic bLocalDataTransferred; + /// Indicates that the remote Tree data is on the GPU + std::atomic bRemoteDataTransferred; +protected: + /// Indicates that Ewald is done or is not needed. + bool bEwaldDone; + /// Counter for PEs that are ready to get their acclerations updated. + int PEsWantParticlesBack; + /// Keep track of which PEs have TreePieces + std::set registeredPEs; /// Reference count for Pieces that have finished updating /// their acclerations. int treePiecesParticlesUpdated; @@ -111,11 +133,13 @@ class DataManager : public CBase_DataManager { // TreePiece counter for multi-threaded GPU host buffer copy int treePiecesBufferFilled; - // can the gpu accept a chunk of remote particles/nodes? - bool gpuFree; + // Request for local GPU tree walk + CudaRequest *lwReq; - /// Callback pointer to pass to HAPI. + /// Callback pointers to pass to HAPI. CkCallback *localTransferCallback; + CkCallback *localWalkCallback; + CkCallback *ewaldCallback; PendingBuffers *currentChunkBuffers; // queue that stores all pending chunk transfers @@ -135,21 +159,11 @@ class DataManager : public CBase_DataManager { CudaMultipoleMoments *bufLocalMoments; /// host buffer to transfer local particles to GPU CompactPartData *bufLocalParts; - /// host buffer to transfer initial accelerations to GPU - VariablePartData *bufLocalVars; - // Pointers to particle and tree data on GPU - CudaMultipoleMoments *d_localMoments; - CudaMultipoleMoments *d_remoteMoments; - CompactPartData *d_localParts; - CompactPartData *d_remoteParts; - VariablePartData *d_localVars; - size_t sMoments; - size_t sCompactParts; - size_t sVarParts; + EwtData *ewtGPU; + EwaldReadOnlyData *cachedData; - int numStreams; - cudaStream_t *streams; + cudaStream_t stream; // For data transfers and local tree walk #endif /// The root of the combined trees @@ -166,6 +180,8 @@ class DataManager : public CBase_DataManager { Tree::NodeLookupType chunkRootTable; public: + CProxy_TreePiece getTreePieces() { return treePieces; } + CkVec getRegisteredTreePieces() { return registeredTreePieces; } /* ** Cooling @@ -186,13 +202,53 @@ class DataManager : public CBase_DataManager { /// @brief Lock for accessing hmstarlog from TreePieces CmiNodeLock lockHMStarLog; +#ifdef CUDA + // Pointers to particle and tree data on GPU + CudaMultipoleMoments *d_localMoments; + CudaMultipoleMoments *d_remoteMoments; + CompactPartData *d_localParts; + CompactPartData *d_remoteParts; + VariablePartData *d_localVars; + size_t sMoments; + size_t sCompactParts; + size_t sVarParts; + + /// @brief log of CUDA memory events. + MemLog *memLog; + /// @brief Lock for accessing memlog from CUDA wrappers + CmiNodeLock lockMemLog; + /// @brief Flag to enable GPU memory logging + int bGpuMemLogger; + + /// @brief log of CPU (host) memory events. + MemLog *cpuMemLog; + /// @brief Lock for accessing CPU memlog + CmiNodeLock lockCpuMemLog; + /// @brief Flag to enable CPU memory logging + int bCpuMemLogger; + /// @brief Target capacity for host pool trim in GB + double dHostPoolTargetCapacityGB; + /// @brief Minimum capacity per bucket during trim in MB + int nHostPoolMinCapacityPerBucketMB; + /// @brief Flag to enable host pool debug output + int bHostPoolDebug; + + /// Allow TreePiece::fillGPUBuffer to access node-wide buffers and data + CudaMultipoleMoments* getLocalMoments() { return localMoments.getVec(); } + CudaMultipoleMoments* getBufLocalMoments() { return bufLocalMoments; } + CompactPartData* getBufLocalParts() { return bufLocalParts; } +#endif + DataManager(const CkArrayID& treePieceID); DataManager(CkMigrateMessage *); void startLocalWalk(); + void finishLocalWalk(); void resumeRemoteChunk(); #ifdef CUDA - void createStreams(int _numStreams, const CkCallback& cb); + void startEwaldGPU(); + void finishEwaldGPU(); + void createStream(const CkCallback& cb); void donePrefetch(int chunk); // serialize remote chunk wrapper void serializeLocalTree(); @@ -203,10 +259,10 @@ class DataManager : public CBase_DataManager { // actual serialization methods PendingBuffers *serializeRemoteChunk(GenericTreeNode *); void serializeLocal(GenericTreeNode *); - void transferLocalToGPU(int nParts, GenericTreeNode *node); + void transferLocalToGPU(int nParts); void freeLocalTreeMemory(); void freeRemoteChunkMemory(int chunk); - void transferParticleVarsBack(); + void transferParticleVarsBack(bool bFromEwald); void updateParticles(UpdateParticlesStruct *data); void updateParticlesFreeMemory(UpdateParticlesStruct *data); void initiateNextChunkTransfer(); @@ -225,6 +281,10 @@ class DataManager : public CBase_DataManager { } nodeTable.clear(); +#ifdef CUDA + cudaStreamDestroy(stream); +#endif + CoolFinalize(Cool); LymanWernerTableFinalize(LWData); delete starLog; @@ -232,10 +292,19 @@ class DataManager : public CBase_DataManager { CmiDestroyLock(lockStarLog); CmiDestroyLock(lockHMStarLog); #ifdef CUDA - for (int i = 0; i < numStreams; i++) { - cudaStreamDestroy(streams[i]); - } - delete[] streams; + if(ewtGPU != nullptr) { +#ifdef PINNED_HOST_MEMORY + const char* funcTag = "DataManager::finishEwaldGPU"; + hostFree(ewtGPU, funcTag); + hostFree(cachedData, funcTag); +#else + free(ewtGPU); + free(cachedData); +#endif + } + + delete memLog; + CmiDestroyLock(lockMemLog); #endif } @@ -285,10 +354,32 @@ class DataManager : public CBase_DataManager { else return NULL; } inline Tree::GenericTreeNode *getRoot() { return root; } + + /// table for the Ewald h-loop + std::vector ewt; +#ifdef HEXADECAPOLE + MOMC momcRoot; /* complete moments of root */ +#endif + /// Periodic size of box + Vector3D fPeriod; + /// Limit of Ewald h-loop + double dEwhCut; + double fEwCut; + int nReplicas; + void EwaldInit(); + void initCooling(double dGmPerCcUnit, double dComovingGmPerCcUnit, double dErgPerGmUnit, double dSecUnit, double dKpcUnit, COOLPARAM inParam, const CkCallback& cb); void initStarLog(std::string _fileName, const CkCallback &cb); + void initMemLog(std::string _fileName, int bGpuMemLoggerFlag, const CkCallback &cb); + void initCpuMemLog(std::string _fileName, int bCpuMemLoggerFlag, const CkCallback &cb); +#ifdef CUDA + void flushMemLog(const CkCallback& cb); + void flushCpuMemLog(const CkCallback& cb); + void trimHostPool(double targetCapacityGB, const CkCallback& cb); + void refillHostPool(const CkCallback& cb); +#endif void initLWData(const CkCallback& cb); void initHMStarLog(std::string _fileName, const CkCallback &cb); void dmCoolTableRead(double *dTableData, int nData, const CkCallback& cb); @@ -329,8 +420,8 @@ class ProjectionsControl : public CBase_ProjectionsControl { #ifdef CUDA // GPUs are assigned to nodes in a round-robin fashion. This allows the user to define // one virtual node per device and utilize multiple GPUs on a single node - // Beacuse devices are assigned per-PE, this is a convenient place to call setDevice - // Note that this code has nothing to do with initalizing projections + // Because devices are assigned per-PE, this is a convenient place to call setDevice + // Note that this code has nothing to do with initializing projections int numGpus; cudaGetDeviceCount(&numGpus); cudaSetDevice(CmiMyNode() % numGpus); @@ -344,6 +435,15 @@ class ProjectionsControl : public CBase_ProjectionsControl { #endif } ProjectionsControl(CkMigrateMessage *m) : CBase_ProjectionsControl(m) { +#ifdef CUDA + // GPUs are assigned to nodes in a round-robin fashion. This allows the user to define + // one virtual node per device and utilize multiple GPUs on a single node + // Because devices are assigned per-PE, this is a convenient place to call setDevice + // Note that this code has nothing to do with initializing projections + int numGpus; + cudaGetDeviceCount(&numGpus); + cudaSetDevice(CmiMyNode() % numGpus); +#endif setBIconfig(); LBTurnCommOff(); #ifndef LB_MANAGER_VERSION diff --git a/Ewald.cpp b/Ewald.cpp index 94ac2b61..ccd96720 100644 --- a/Ewald.cpp +++ b/Ewald.cpp @@ -74,7 +74,7 @@ void TreePiece::BucketEwald(GenericTreeNode *req, int nReps,double fEwCut) #ifndef BENCHMARK_NO_WORK GravityParticle *p; #ifdef HEXADECAPOLE - MOMC mom = momcRoot; + MOMC mom = dm->momcRoot; MultipoleMoments momQuad = root->moments; double xx,xxx,xxy,xxz,yy,yyy,yyz,xyy,zz,zzz,xzz,yzz,xy,xyz,xz,yz; double Q4mirx,Q4miry,Q4mirz,Q4mir,Q4x,Q4y,Q4z; @@ -262,7 +262,8 @@ void TreePiece::BucketEwald(GenericTreeNode *req, int nReps,double fEwCut) ** Total = (22,36) ** = 58 */ - for (i=0;i &ewt = dm->ewt; + for (i=0;imoments.mom),1.0f,root->moments.getRadius()); @@ -305,24 +305,13 @@ void TreePiece::EwaldInit() L = fPeriod.x; alpha = 2.0/L; k4 = M_PI*M_PI/(alpha*alpha*L*L); - i = 0; + ewt.clear(); for (hx=-hReps;hx<=hReps;++hx) { for (hy=-hReps;hy<=hReps;++hy) { for (hz=-hReps;hz<=hReps;++hz) { h2 = hx*hx + hy*hy + hz*hz; if (h2 == 0) continue; if (h2 > dEwhCut*dEwhCut) continue; - if (i == nMaxEwhLoop) { - nMaxEwhLoop *= 2; - /* avoid realloc() */ - EWT *ewtTmp = new EWT[nMaxEwhLoop]; - assert(ewtTmp != NULL); - for(int j = 0; j < i; j++) { - ewtTmp[j] = ewt[j]; - } - delete[] ewt; - ewt = ewtTmp; - } gam[0] = exp(-k4*h2)/(M_PI*h2*L); gam[1] = 2*M_PI/L*gam[0]; gam[2] = -2*M_PI/L*gam[1]; @@ -363,191 +352,15 @@ void TreePiece::EwaldInit() QEVAL(root->moments, gam,hx,hy,hz, ax,ay,az,mfacs); #endif - ewt[i].hx = 2*M_PI/L*hx; - ewt[i].hy = 2*M_PI/L*hy; - ewt[i].hz = 2*M_PI/L*hz; - ewt[i].hCfac = mfacc; - ewt[i].hSfac = mfacs; - ++i; + EWT ewtTmp; + ewtTmp.hx = 2*M_PI/L*hx; + ewtTmp.hy = 2*M_PI/L*hy; + ewtTmp.hz = 2*M_PI/L*hz; + ewtTmp.hCfac = mfacc; + ewtTmp.hSfac = mfacs; + ewt.push_back(ewtTmp); } } } - nEwhLoop = i; - EwaldMsg *msg = new (8*sizeof(int)) EwaldMsg; - msg->fromInit = true; - // Make priority lower than gravity or smooth. - *((int *)CkPriorityPtr(msg)) = 3*numTreePieces * numChunks + thisIndex + 1; - CkSetQueueing(msg,CK_QUEUEING_IFIFO); - thisProxy[thisIndex].calculateEwald(msg); } - - -/// @brief Transfer Ewald data to GPU memory and launch EwaldHost kernel -void TreePiece::EwaldGPU(){ - /* when not using CUDA, definition is required because - EwaldGPU is an entry method - */ - -#ifdef SPCUDA - if(NumberOfGPUParticles == 0 || myNumActiveParticles == 0){ - for (int i=0; imoments; - - cudatype L = fPeriod.x; - cudatype alpha = 2.0f/L; - - ewtTable = (EwtData*) h_idata->ewt; - roData = (EwaldReadOnlyData*) h_idata->cachedData; - - int nActive = 0; - if(largephase){ - int *markers = (int *) h_idata->EwaldMarkers; - int IDX = 0; - int IDXend = 0; - int IDXstart = 0; - for(int b = 0; b < numBuckets; b++){ - GenericTreeNode *bucketNode = bucketList[b]; - GravityParticle *buckparts = bucketNode->particlePointer; - IDXstart = bucketNode->firstParticle; - IDXend = bucketNode->lastParticle; - IDX = bucketNode->bucketArrayIndex; /*First particle index on GPU*/ - for(int j = IDXstart; j <= IDXend; j++){ /*Go through all particles in bucket*/ - if(buckparts[j - IDXstart].rung < activeRung){IDX++; continue;} - markers[nActive] = IDX; - IDX++; - nActive++; - } - } - h_idata->EwaldRange[0] = 0; - h_idata->EwaldRange[1] = nActive - 1; - assert(nActive <= NumberOfGPUParticles); - assert(nActive == myNumActiveParticles); - }else{ - nActive = NumberOfGPUParticles; - h_idata->EwaldRange[0] = FirstGPUParticleIndex; - h_idata->EwaldRange[1] = LastGPUParticleIndex; - } - - for (int i=0; imomcRoot.m = (cudatype) momcRoot.m; - roData->momcRoot.xx = (cudatype) momcRoot.xx; - roData->momcRoot.yy = (cudatype) momcRoot.yy; - roData->momcRoot.xy = (cudatype) momcRoot.xy; - roData->momcRoot.xz = (cudatype) momcRoot.xz; - roData->momcRoot.yz = (cudatype) momcRoot.yz; - roData->momcRoot.xxx = (cudatype) momcRoot.xxx; - roData->momcRoot.xyy = (cudatype) momcRoot.xyy; - roData->momcRoot.xxy = (cudatype) momcRoot.xxy; - roData->momcRoot.yyy = (cudatype) momcRoot.yyy; - roData->momcRoot.xxz = (cudatype) momcRoot.xxz; - roData->momcRoot.yyz = (cudatype) momcRoot.yyz; - roData->momcRoot.xyz = (cudatype) momcRoot.xyz; - roData->momcRoot.xxxx = (cudatype) momcRoot.xxxx; - roData->momcRoot.xyyy = (cudatype) momcRoot.xyyy; - roData->momcRoot.xxxy = (cudatype) momcRoot.xxxy; - roData->momcRoot.yyyy = (cudatype) momcRoot.yyyy; - roData->momcRoot.xxxz = (cudatype) momcRoot.xxxz; - roData->momcRoot.yyyz = (cudatype) momcRoot.yyyz; - roData->momcRoot.xxyy = (cudatype) momcRoot.xxyy; - roData->momcRoot.xxyz = (cudatype) momcRoot.xxyz; - roData->momcRoot.xyyz = (cudatype) momcRoot.xyyz; - roData->momcRoot.zz = (cudatype) momcRoot.zz; - roData->momcRoot.xzz = (cudatype) momcRoot.xzz; - roData->momcRoot.yzz = (cudatype) momcRoot.yzz; - roData->momcRoot.zzz = (cudatype) momcRoot.zzz; - roData->momcRoot.xxzz = (cudatype) momcRoot.xxzz; - roData->momcRoot.xyzz = (cudatype) momcRoot.xyzz; - roData->momcRoot.xzzz = (cudatype) momcRoot.xzzz; - roData->momcRoot.yyzz = (cudatype) momcRoot.yyzz; - roData->momcRoot.yzzz = (cudatype) momcRoot.yzzz; - roData->momcRoot.zzzz = (cudatype) momcRoot.zzzz; -#else - roData->mm.xx = (cudatype) mm.xx; - roData->mm.xy = (cudatype) mm.xy; - roData->mm.xz = (cudatype) mm.xz; - roData->mm.yy = (cudatype) mm.yy; - roData->mm.yz = (cudatype) mm.yz; - roData->mm.zz = (cudatype) mm.zz; -#endif - roData->mm.totalMass = (cudatype) mm.totalMass; - roData->mm.cmx = (cudatype) mm.cm.x; - roData->mm.cmy = (cudatype) mm.cm.y; - roData->mm.cmz = (cudatype) mm.cm.z; - roData->n = nActive; - roData->fEwCut = (cudatype) fEwCut; - roData->nReps = nReplicas ; - roData->nEwReps = (int) ceil(fEwCut); - roData->nEwhLoop = nEwhLoop; - roData->L = (cudatype) L; - roData->alpha = alpha; - roData->alpha2 = (cudatype) alpha*alpha; - roData->k1 = (cudatype) M_PI/(alpha*alpha*L*L*L); - roData->ka = (cudatype) 2.0*alpha/sqrt(M_PI); - roData->fEwCut2 = (cudatype) fEwCut*fEwCut*L*L; -/* - Break between Taylor expansion for small r and multipole expansion. - This value is for double precision. - roData->fInner2 = (cudatype) 1.2e-3*L*L; - The following is for single precision. The CUDA version currently uses - erff() and erfcf(). If these ever get changed then the following line - needs to be changed accordingly. - */ - roData->fInner2 = (cudatype) 1.1e-2*L*L; - - CkArrayIndex1D myIndex = CkArrayIndex1D(thisIndex); - cbEwaldGPU = new CkCallback(CkIndex_TreePiece::EwaldGPUComplete(), myIndex, - thisArrayID); - - int myLocalIndex; - for(myLocalIndex = 0; this != dm->registeredTreePieces[myLocalIndex].treePiece; - myLocalIndex++); - CkAssert(myLocalIndex < dm->registeredTreePieces.length()); - - EwaldHost(this->d_localParts, this->d_localVars, - h_idata, this->stream, (void *) cbEwaldGPU, myLocalIndex, largephase); -#endif -} - -/// @brief Callback for EwaldGPU. Clean up device memory + host buffer and call finishBucket -void TreePiece::EwaldGPUComplete() { - /* when not using CUDA, definition is required because - EwaldGPUComplete is an entry method - */ -#ifdef SPCUDA - int largephase = largePhase(); - EwaldHostMemoryFree(h_idata, largephase); - free(h_idata); - delete cbEwaldGPU; - - /* indicate completion of ewald */ - - for (int i=0; i>>( +#ifndef CUDA_NO_KERNELS + int numBlocks = (numParticles + (THREADS_PER_BLOCK - 1)) / THREADS_PER_BLOCK; + ZeroVars<<>>( (VariablePartData *) *d_varParts, numParticles); - cudaChk(cudaPeekAtLastError()); +#endif HAPI_TRACE_END(CUDA_XFER_LOCAL); @@ -162,7 +146,7 @@ void DataManagerTransferRemoteChunk(void *moments, size_t sMoments, cudaStream_t stream, void *callback) { -#ifdef CUDA_VERBOSE_KERNEL_ENQUEUE +#ifdef CUDA_VERBOSE_OPS printf("(%d) DM REMOTE CHUNK moments %zu partcores %zu\n", CmiMyPe(), sMoments, @@ -172,10 +156,14 @@ void DataManagerTransferRemoteChunk(void *moments, size_t sMoments, HAPI_TRACE_BEGIN(); - cudaChk(cudaMalloc(d_remoteMoments, sMoments)); - cudaChk(cudaMalloc(d_remoteParts, sRemoteParts)); + const char* funcTag = "DataManagerTransferRemoteChunk"; + cudaChk(gpuPoolMalloc(d_remoteMoments, sMoments, stream, funcTag)); + cudaChk(gpuPoolMalloc(d_remoteParts, sRemoteParts, stream, funcTag)); cudaChk(cudaMemcpyAsync(*d_remoteMoments, moments, sMoments, cudaMemcpyHostToDevice, stream)); cudaChk(cudaMemcpyAsync(*d_remoteParts, remoteParts, sRemoteParts, cudaMemcpyHostToDevice, stream)); + // Synchronize stream so that d_remoteMoments and d_remoteParts will be + // available to other streams + cudaChk(cudaStreamSynchronize(stream)); HAPI_TRACE_END(CUDA_XFER_REMOTE); @@ -184,30 +172,23 @@ void DataManagerTransferRemoteChunk(void *moments, size_t sMoments, /************** Gravity *****************/ -/// @brief Initiate a local gravity calculation on the GPU, via an interaction -/// list calculation between nodes, or do a local tree walk + +/// @brief Initiate a local gravity calculation on the GPU from the DataManager +/// The calculation is done for all local particles at once /// @param data CudaRequest object containing parameters for the calculation -void TreePieceCellListDataTransferLocal(CudaRequest *data){ +void DataManagerLocalTreeWalk(CudaRequest *data){ cudaStream_t stream = data->stream; - CudaDevPtr devPtr; - TreePieceDataTransferBasic(data, &devPtr); - -#ifdef CUDA_VERBOSE_KERNEL_ENQUEUE - printf("(%d) TRANSFER LOCAL CELL\n", CmiMyPe()); -#endif -#ifdef CUDA_NOTIFY_DATA_TRANSFER_DONE - printf("TRANSFER LOCAL CELL KERNELSELECT buffers:\nlocal_particles: (0x%x)\nlocal_particle_vars: (0x%x)\nremote_moments: (0x%x)\nil_cell: (0x%x)\n", - data->d_localParts, - data->d_localVars, - data->d_localMoments, - devPtr.d_list - ); +#ifdef CUDA_VERBOSE_OPS + printf("(%d) DM LOCAL TREE WALK numParts: %d\n", + CmiMyPe(), + data->lastParticle - data->firstParticle + ); #endif HAPI_TRACE_BEGIN(); -#ifndef CUDA_NO_KERNELS #ifdef GPU_LOCAL_TREE_WALK +#ifndef CUDA_NO_KERNELS gpuLocalTreeWalk<<<(data->lastParticle - data->firstParticle + 1) / THREADS_PER_BLOCK + 1, dim3(THREADS_PER_BLOCK), 0, stream>>> ( data->d_localMoments, @@ -223,108 +204,77 @@ void TreePieceCellListDataTransferLocal(CudaRequest *data){ data->fperiodY, data->fperiodZ ); - #else - dim3 dimensions = THREADS_PER_BLOCK; - #ifdef CUDA_2D_TB_KERNEL - dimensions = dim3(NODES_PER_BLOCK, PARTS_PER_BLOCK); - #endif - nodeGravityComputation<<numBucketsPlusOne-1), dimensions, 0, stream>>> ( - data->d_localParts, - data->d_localVars, - data->d_localMoments, - (ILCell *)devPtr.d_list, - devPtr.d_bucketMarkers, - devPtr.d_bucketStarts, - devPtr.d_bucketSizes, - data->fperiod - ); #endif #endif - TreePieceDataTransferBasicCleanup(&devPtr); - cudaChk(cudaPeekAtLastError()); - HAPI_TRACE_END(CUDA_GRAV_LOCAL); + HAPI_TRACE_END(CUDA_GRAV_TREE_LOCAL); hapiAddCallback(stream, data->cb); } -/// @brief Initiate a remote gravity calculation on the GPU between tree nodes +/// @brief Initiate a local particle-node interaction list calculation on the +/// GPU from a PE /// @param data CudaRequest object containing parameters for the calculation -void TreePieceCellListDataTransferRemote(CudaRequest *data){ +void PEListNodeListDataTransferLocal(CudaRequest *data){ cudaStream_t stream = data->stream; - CudaDevPtr devPtr; - TreePieceDataTransferBasic(data, &devPtr); -#ifdef CUDA_VERBOSE_KERNEL_ENQUEUE - printf("(%d) TRANSFER REMOTE CELL\n", CmiMyPe()); -#endif - -#ifdef CUDA_NOTIFY_DATA_TRANSFER_DONE - printf("TRANSFER REMOTE CELL KERNELSELECT buffers:\nlocal_particles: (0x%x)\nlocal_particle_vars: (0x%x)\nremote_moments: (0x%x)\nil_cell: (0x%x)\n", - data->d_localParts, - data->d_localVars, - data->d_remoteMoments, - devPtr.d_list +#ifdef CUDA_VERBOSE_OPS + printf("(%d) PE NODE LIST LOCAL numInteractions: %d\n", + CmiMyPe(), + data->numInteractions ); #endif - HAPI_TRACE_BEGIN(); + if (data->numInteractions > 0) { + HAPI_TRACE_BEGIN(); + CudaDevPtr devPtr = {nullptr, nullptr, nullptr, nullptr}; + const char* funcTag = "PEListNodeListDataTransferLocal"; + DataTransferBasic(data, &devPtr, funcTag); + + dim3 dimensions = dim3(NODES_PER_BLOCK, PARTS_PER_BLOCK); #ifndef CUDA_NO_KERNELS - dim3 dimensions = THREADS_PER_BLOCK; -#ifdef CUDA_2D_TB_KERNEL - dimensions = dim3(NODES_PER_BLOCK, PARTS_PER_BLOCK); -#endif - nodeGravityComputation<<numBucketsPlusOne-1, dimensions, 0, stream>>> ( - data->d_localParts, - data->d_localVars, - data->d_remoteMoments, - (ILCell *)devPtr.d_list, - devPtr.d_bucketMarkers, - devPtr.d_bucketStarts, - devPtr.d_bucketSizes, - data->fperiod - ); + nodeGravityComputation<<numBucketsPlusOne-1), dimensions, 0, stream>>> ( + data->d_localParts, + data->d_localVars, + data->d_localMoments, + (ILCell *)devPtr.d_list, + devPtr.d_bucketMarkers, + devPtr.d_bucketStarts, + devPtr.d_bucketSizes, + data->fperiod + ); #endif - TreePieceDataTransferBasicCleanup(&devPtr); - cudaChk(cudaPeekAtLastError()); - HAPI_TRACE_END(CUDA_GRAV_REMOTE); + + DataTransferBasicCleanup(&devPtr, stream, funcTag); + HAPI_TRACE_END(CUDA_GRAV_NODELIST_LOCAL); + } hapiAddCallback(stream, data->cb); } -/// @brief Initiate a remote resume gravity calculation on the GPU between tree nodes +/// @brief Initiate a local particle-particle interaction list calculation on the +/// GPU from a PE /// @param data CudaRequest object containing parameters for the calculation -void TreePieceCellListDataTransferRemoteResume(CudaRequest *data){ +void PEListPartListDataTransferLocal(CudaRequest *data){ cudaStream_t stream = data->stream; - CudaDevPtr devPtr; - void *d_missedNodes; - TreePieceDataTransferBasic(data, &devPtr); -#ifdef CUDA_VERBOSE_KERNEL_ENQUEUE - printf("(%d) TRANSFER REMOTE RESUME CELL\n", CmiMyPe()); +#ifdef CUDA_VERBOSE_OPS + printf("(%d) PE PART LIST LOCAL numInteractions: %d\n", + CmiMyPe(), + data->numInteractions + ); #endif - cudaChk(cudaMalloc(&d_missedNodes, data->sMissed)); - cudaChk(cudaMemcpyAsync(d_missedNodes, data->missedNodes, data->sMissed, cudaMemcpyHostToDevice, stream)); + if (data->numInteractions > 0) { + HAPI_TRACE_BEGIN(); + CudaDevPtr devPtr = {nullptr, nullptr, nullptr, nullptr}; + const char* funcTag = "PEListPartListDataTransferLocal"; + DataTransferBasic(data, &devPtr, funcTag); -#ifdef CUDA_NOTIFY_DATA_TRANSFER_DONE - printf("TRANSFER REMOTE RESUME CELL KERNELSELECT buffers:\nlocal_particles: (0x%x)\nlocal_particle_vars: (0x%x)\nmissed_moments (0x%x)\nil_cell (0x%x)\n", - data->d_localParts, - data->d_localVars, - d_missedNodes, - devPtr.d_list - ); -#endif - - HAPI_TRACE_BEGIN(); #ifndef CUDA_NO_KERNELS - dim3 dimensions = THREADS_PER_BLOCK; -#ifdef CUDA_2D_TB_KERNEL - dimensions = dim3(NODES_PER_BLOCK, PARTS_PER_BLOCK); -#endif - nodeGravityComputation<<numBucketsPlusOne-1, dimensions, 0, stream>>> ( + particleGravityComputation<<numBucketsPlusOne-1, dim3(NODES_PER_BLOCK_PART, PARTS_PER_BLOCK_PART), 0, stream>>> ( data->d_localParts, data->d_localVars, - (CudaMultipoleMoments *)d_missedNodes, + data->d_localParts, (ILCell *)devPtr.d_list, devPtr.d_bucketMarkers, devPtr.d_bucketStarts, @@ -332,241 +282,191 @@ void TreePieceCellListDataTransferRemoteResume(CudaRequest *data){ data->fperiod ); #endif - TreePieceDataTransferBasicCleanup(&devPtr); - cudaChk(cudaFree(d_missedNodes)); - cudaChk(cudaPeekAtLastError()); - HAPI_TRACE_END(CUDA_REMOTE_RESUME); + DataTransferBasicCleanup(&devPtr, stream, funcTag); + HAPI_TRACE_END(CUDA_GRAV_PARTLIST_LOCAL); + } hapiAddCallback(stream, data->cb); } -/// @brief Initiate a small phase local gravity calculation on the GPU between particles +/// @brief Initiate a remote particle-node interaction list calculation on the +/// GPU from a PE /// @param data CudaRequest object containing parameters for the calculation -void TreePiecePartListDataTransferLocalSmallPhase(CudaRequest *data, CompactPartData *particles, int len){ +void PEListNodeListDataTransferRemote(CudaRequest *data){ cudaStream_t stream = data->stream; - CudaDevPtr devPtr; - TreePieceDataTransferBasic(data, &devPtr); - size_t size = (len) * sizeof(CompactPartData); - void* bufferHostBuffer; - void* d_smallParts; +#ifdef CUDA_VERBOSE_OPS + printf("(%d) PE NODE LIST REMOTE numInteractions: %d\n", + CmiMyPe(), + data->numInteractions + ); +#endif -#ifdef CUDA_NOTIFY_DATA_TRANSFER_DONE - printf("TreePiecePartListDataTransferLocalSmallPhase KERNELSELECT buffers:\nlocal_particles: (0x%x)\nlocal_particle_vars: (0x%x)\nil_cell: (0x%x)\n", + if (data->numInteractions > 0) { + HAPI_TRACE_BEGIN(); + CudaDevPtr devPtr = {nullptr, nullptr, nullptr, nullptr}; + const char* funcTag = "PEListNodeListDataTransferRemote"; + DataTransferBasic(data, &devPtr, funcTag); + + dim3 dimensions = dim3(NODES_PER_BLOCK, PARTS_PER_BLOCK); +#ifndef CUDA_NO_KERNELS + nodeGravityComputation<<numBucketsPlusOne-1, dimensions, 0, stream>>> ( data->d_localParts, data->d_localVars, - devPtr.d_list - ); -#endif - -#ifdef CUDA_VERBOSE_KERNEL_ENQUEUE - printf("(%d) TRANSFER LOCAL SMALL PHASE %zu\n", - CmiMyPe(), - size + data->d_remoteMoments, + (ILCell *)devPtr.d_list, + devPtr.d_bucketMarkers, + devPtr.d_bucketStarts, + devPtr.d_bucketSizes, + data->fperiod ); #endif + DataTransferBasicCleanup(&devPtr, stream, funcTag); + HAPI_TRACE_END(CUDA_GRAV_NODELIST_REMOTE); + } - HAPI_TRACE_BEGIN(); - allocatePinnedHostMemory(&bufferHostBuffer, size); -#ifdef CUDA_PRINT_ERRORS - printf("TPPartSmallPhase 0: %s\n", cudaGetErrorString( cudaGetLastError() ) ); -#endif - memcpy(bufferHostBuffer, particles, size); - cudaChk(cudaMalloc(&d_smallParts, size)); - cudaChk(cudaMemcpyAsync(d_smallParts, bufferHostBuffer, size, cudaMemcpyHostToDevice, stream)); - -#ifndef CUDA_NO_KERNELS -#ifdef CUDA_2D_TB_KERNEL - particleGravityComputation<<numBucketsPlusOne-1, dim3(NODES_PER_BLOCK_PART, PARTS_PER_BLOCK_PART), 0, stream>>> ( - data->d_localParts, - data->d_localVars, - (CompactPartData *)d_smallParts, - (ILCell *)devPtr.d_list, - devPtr.d_bucketMarkers, - devPtr.d_bucketStarts, - devPtr.d_bucketSizes, - data->fperiod - ); -#else - particleGravityComputation<<numBucketsPlusOne-1, THREADS_PER_BLOCK, 0, stream>>> ( - data->d_localParts, - data->d_localVars, - (CompactPartData *)d_smallParts, - (ILCell *)devPtr.d_list, - devPtr.d_bucketMarkers, - devPtr.d_bucketStarts, - devPtr.d_bucketSizes, - data->fperiod - ); -#endif -#endif - TreePieceDataTransferBasicCleanup(&devPtr); - cudaChk(cudaPeekAtLastError()); - HAPI_TRACE_END(CUDA_PART_GRAV_LOCAL_SMALL); - cudaChk(cudaFree(d_smallParts)); hapiAddCallback(stream, data->cb); } -/// @brief Initiate a local gravity calculation on the GPU between particles +/// @brief Initiate a remote-resume particle-node interaction list calculation on +/// the GPU from a PE /// @param data CudaRequest object containing parameters for the calculation -void TreePiecePartListDataTransferLocal(CudaRequest *data){ +void PEListNodeListDataTransferRemoteResume(CudaRequest *data){ cudaStream_t stream = data->stream; - CudaDevPtr devPtr; - TreePieceDataTransferBasic(data, &devPtr); - -#ifdef CUDA_VERBOSE_KERNEL_ENQUEUE - printf("(%d) TRANSFER LOCAL LARGEPHASE PART\n", CmiMyPe()); -#endif -#ifdef CUDA_NOTIFY_DATA_TRANSFER_DONE - printf("TreePiecePartListDataTransferLocal buffers:\nlocal_particles: (0x%x)\nlocal_particle_vars: (0x%x)\nil_cell: (0x%x)\n", - data->d_localParts, - data->d_localVars, - devPtr.d_list +#ifdef CUDA_VERBOSE_OPS + printf("(%d) PE NODE LIST REMOTE RESUME numInteractions: %d\n", + CmiMyPe(), + data->numInteractions ); #endif - HAPI_TRACE_BEGIN(); + if (data->numInteractions > 0) { + HAPI_TRACE_BEGIN(); + CudaDevPtr devPtr = {nullptr, nullptr, nullptr, nullptr}; + void* d_missedNodes = nullptr; + const char* funcTag = "PEListNodeListDataTransferRemoteResume"; + DataTransferBasic(data, &devPtr, funcTag); + + cudaChk(gpuPoolMalloc(&d_missedNodes, data->sMissed, stream, funcTag)); + cudaChk(cudaMemcpyAsync(d_missedNodes, data->missedNodes, data->sMissed, cudaMemcpyHostToDevice, stream)); + + dim3 dimensions = dim3(NODES_PER_BLOCK, PARTS_PER_BLOCK); #ifndef CUDA_NO_KERNELS -#ifdef CUDA_2D_TB_KERNEL - particleGravityComputation<<numBucketsPlusOne-1, dim3(NODES_PER_BLOCK_PART, PARTS_PER_BLOCK_PART), 0, stream>>> ( - data->d_localParts, - data->d_localVars, - data->d_localParts, - (ILCell *)devPtr.d_list, - devPtr.d_bucketMarkers, - devPtr.d_bucketStarts, - devPtr.d_bucketSizes, - data->fperiod - ); -#else - particleGravityComputation<<numBucketsPlusOne-1, THREADS_PER_BLOCK, 0, stream>>> ( - data->d_localParts, - data->d_localVars, - data->d_localParts, - (ILPart *)devPtr.d_list, - devPtr.d_bucketMarkers, - devPtr.d_bucketStarts, - devPtr.d_bucketSizes, - data->fperiod - ); -#endif + nodeGravityComputation<<numBucketsPlusOne-1, dimensions, 0, stream>>> ( + data->d_localParts, + data->d_localVars, + (CudaMultipoleMoments *)d_missedNodes, + (ILCell *)devPtr.d_list, + devPtr.d_bucketMarkers, + devPtr.d_bucketStarts, + devPtr.d_bucketSizes, + data->fperiod + ); #endif - TreePieceDataTransferBasicCleanup(&devPtr); - cudaChk(cudaPeekAtLastError()); - HAPI_TRACE_END(CUDA_PART_GRAV_LOCAL); + + DataTransferBasicCleanup(&devPtr, stream, funcTag); + if (d_missedNodes != nullptr) { + cudaChk(gpuPoolFree(d_missedNodes, stream, funcTag)); + } + HAPI_TRACE_END(CUDA_GRAV_NODELIST_REMOTE_RESUME); + } hapiAddCallback(stream, data->cb); } -/// @brief Initiate a remote gravity calculation on the GPU between particles +/// @brief Initiate a remote particle-particle interaction list calculation on the +/// GPU from a PE /// @param data CudaRequest object containing parameters for the calculation -void TreePiecePartListDataTransferRemote(CudaRequest *data){ +void PEListPartListDataTransferRemote(CudaRequest *data){ cudaStream_t stream = data->stream; + +#ifdef CUDA_VERBOSE_OPS + printf("(%d) PE PART LIST REMOTE numInteractions: %d\n", + CmiMyPe(), + data->numInteractions + ); +#endif + CudaDevPtr devPtr; - TreePieceDataTransferBasic(data, &devPtr); + const char* funcTag = "PEListPartListDataTransferRemote"; + if (data->numInteractions > 0) { + DataTransferBasic(data, &devPtr, funcTag); #ifdef CUDA_VERBOSE_KERNEL_ENQUEUE - printf("(%d) TRANSFER REMOTE PART\n", CmiMyPe()); + printf("(%d) TRANSFER REMOTE PART\n", CmiMyPe()); #endif #ifdef CUDA_NOTIFY_DATA_TRANSFER_DONE - printf("TreePiecePartListDataTransferRemote KERNELSELECT buffers:\nlocal_particles: (0x%x)\nlocal_particle_vars: (0x%x)\nil_cell: (0x%x) (0x%x)\n", - data->d_localParts, - data->d_localVars, - data->list, - devPtr.d_list - ); + printf("TreePiecePartListDataTransferRemote KERNELSELECT buffers:\nlocal_particles: (0x%x)\nlocal_particle_vars: (0x%x)\nil_cell: (0x%x) (0x%x)\n", + data->d_localParts, + data->d_localVars, + data->list, + devPtr.d_list + ); #endif - HAPI_TRACE_BEGIN(); + HAPI_TRACE_BEGIN(); #ifndef CUDA_NO_KERNELS -#ifdef CUDA_2D_TB_KERNEL - particleGravityComputation<<numBucketsPlusOne-1, dim3(NODES_PER_BLOCK_PART, PARTS_PER_BLOCK_PART), 0, stream>>> ( - data->d_localParts, - data->d_localVars, - data->d_remoteParts, - (ILCell *)devPtr.d_list, - devPtr.d_bucketMarkers, - devPtr.d_bucketStarts, - devPtr.d_bucketSizes, - data->fperiod - ); -#else - particleGravityComputation<<numBucketsPlusOne-1, THREADS_PER_BLOCK, 0, stream>>> ( - data->d_localParts, - data->d_localVars, - data->d_remoteParts, - (ILPart *)devPtr.d_list, - devPtr.d_bucketMarkers, - devPtr.d_bucketStarts, - devPtr.d_bucketSizes, - data->fperiod - ); -#endif + particleGravityComputation<<numBucketsPlusOne-1, dim3(NODES_PER_BLOCK_PART, PARTS_PER_BLOCK_PART), 0, stream>>> ( + data->d_localParts, + data->d_localVars, + data->d_remoteParts, + (ILCell *)devPtr.d_list, + devPtr.d_bucketMarkers, + devPtr.d_bucketStarts, + devPtr.d_bucketSizes, + data->fperiod + ); #endif - TreePieceDataTransferBasicCleanup(&devPtr); - cudaChk(cudaPeekAtLastError()); - HAPI_TRACE_END(CUDA_PART_GRAV_REMOTE); + DataTransferBasicCleanup(&devPtr, stream, funcTag); + HAPI_TRACE_END(CUDA_GRAV_PARTLIST_REMOTE); + } hapiAddCallback(stream, data->cb); } -/// @brief Initiate a remote gravity calculation on the GPU between particles +/// @brief Initiate a remote-resume particle-particle interaction list calculation on +/// the GPU from a PE /// @param data CudaRequest object containing parameters for the calculation -void TreePiecePartListDataTransferRemoteResume(CudaRequest *data){ +void PEListPartListDataTransferRemoteResume(CudaRequest *data){ cudaStream_t stream = data->stream; - CudaDevPtr devPtr; - void* d_missedParts; - TreePieceDataTransferBasic(data, &devPtr); -#ifdef CUDA_VERBOSE_KERNEL_ENQUEUE - printf("(%d) TRANSFER REMOTE RESUME PART\n", CmiMyPe()); +#ifdef CUDA_VERBOSE_OPS + printf("(%d) PE PART LIST REMOTE RESUME numInteractions: %d\n", + CmiMyPe(), + data->numInteractions + ); #endif - cudaChk(cudaMalloc(&d_missedParts, data->sMissed)); - cudaChk(cudaMemcpyAsync(d_missedParts, data->missedParts, data->sMissed, cudaMemcpyHostToDevice, stream)); + if (data->numInteractions > 0) { + HAPI_TRACE_BEGIN(); + CudaDevPtr devPtr = {nullptr, nullptr, nullptr, nullptr}; + void* d_missedParts = nullptr; + const char* funcTag = "PEListPartListDataTransferRemoteResume"; + DataTransferBasic(data, &devPtr, funcTag); -#ifdef CUDA_NOTIFY_DATA_TRANSFER_DONE - printf("TreePiecePartListDataTransferRemoteResume KERNELSELECT buffers:\nlocal_particles: (0x%x)\nlocal_particle_vars: (0x%x)\nmissed_parts (0x%x)\nil_cell: (0x%x) (0x%x)\n", - data->d_localParts, - data->d_localVars, - (CompactPartData *)d_missedParts, - data->list, - devPtr.d_list - ); -#endif + cudaChk(gpuPoolMalloc(&d_missedParts, data->sMissed, stream, funcTag)); + cudaChk(cudaMemcpyAsync(d_missedParts, data->missedParts, data->sMissed, cudaMemcpyHostToDevice, stream)); - HAPI_TRACE_BEGIN(); #ifndef CUDA_NO_KERNELS -#ifdef CUDA_2D_TB_KERNEL - particleGravityComputation<<numBucketsPlusOne-1, dim3(NODES_PER_BLOCK_PART, PARTS_PER_BLOCK_PART), 0, stream>>> ( - data->d_localParts, - data->d_localVars, - (CompactPartData *)d_missedParts, - (ILCell *)devPtr.d_list, - devPtr.d_bucketMarkers, - devPtr.d_bucketStarts, - devPtr.d_bucketSizes, - data->fperiod - ); -#else - particleGravityComputation<<numBucketsPlusOne-1, THREADS_PER_BLOCK, 0, stream>>> ( - data->d_localParts, - data->d_localVars, - (CompactPartData *)d_missedParts, - (ILPart *)devPtr.d_list, - devPtr.d_bucketMarkers, - devPtr.d_bucketStarts, - devPtr.d_bucketSizes, - data->fperiod - ); -#endif + particleGravityComputation<<numBucketsPlusOne-1, dim3(NODES_PER_BLOCK_PART, PARTS_PER_BLOCK_PART), 0, stream>>> ( + data->d_localParts, + data->d_localVars, + (CompactPartData *)d_missedParts, + (ILCell *)devPtr.d_list, + devPtr.d_bucketMarkers, + devPtr.d_bucketStarts, + devPtr.d_bucketSizes, + data->fperiod + ); #endif - TreePieceDataTransferBasicCleanup(&devPtr); - cudaChk(cudaFree(d_missedParts)); - cudaChk(cudaPeekAtLastError()); - HAPI_TRACE_END(CUDA_PART_GRAV_REMOTE); + DataTransferBasicCleanup(&devPtr, stream, funcTag); + if (d_missedParts != nullptr) { + cudaChk(gpuPoolFree(d_missedParts, stream, funcTag)); + } + HAPI_TRACE_END(CUDA_GRAV_PARTLIST_REMOTE_RESUME); + } hapiAddCallback(stream, data->cb); } @@ -575,26 +475,35 @@ void TreePiecePartListDataTransferRemoteResume(CudaRequest *data){ /// device memory /// @param data CudaRequest object containing parameters for the calculation /// @param ptr CudaDevPtr object that stores handles to device memory -void TreePieceDataTransferBasic(CudaRequest *data, CudaDevPtr *ptr){ +/// @param functionTag String literal identifying the calling function context. +void DataTransferBasic(CudaRequest *data, CudaDevPtr *ptr, const char* functionTag){ cudaStream_t stream = data->stream; + // Initialize all pointers to NULL first to prevent freeing uninitialized pointers + if (ptr != nullptr) { + ptr->d_list = nullptr; + ptr->d_bucketMarkers = nullptr; + ptr->d_bucketStarts = nullptr; + ptr->d_bucketSizes = nullptr; + } + int numBucketsPlusOne = data->numBucketsPlusOne; int numBuckets = numBucketsPlusOne-1; size_t listSize = (data->numInteractions) * sizeof(ILCell); size_t markerSize = (numBucketsPlusOne) * sizeof(int); size_t startSize = (numBuckets) * sizeof(int); - cudaChk(cudaMalloc(&ptr->d_list, listSize)); - cudaChk(cudaMalloc(&ptr->d_bucketMarkers, markerSize)); - cudaChk(cudaMalloc(&ptr->d_bucketStarts, startSize)); - cudaChk(cudaMalloc(&ptr->d_bucketSizes, startSize)); + cudaChk(gpuPoolMalloc(&ptr->d_list, listSize, stream, functionTag)); + cudaChk(gpuPoolMalloc(&ptr->d_bucketMarkers, markerSize, stream, functionTag)); + cudaChk(gpuPoolMalloc(&ptr->d_bucketStarts, startSize, stream, functionTag)); + cudaChk(gpuPoolMalloc(&ptr->d_bucketSizes, startSize, stream, functionTag)); cudaChk(cudaMemcpyAsync(ptr->d_list, data->list, listSize, cudaMemcpyHostToDevice, stream)); cudaChk(cudaMemcpyAsync(ptr->d_bucketMarkers, data->bucketMarkers, markerSize, cudaMemcpyHostToDevice, stream)); cudaChk(cudaMemcpyAsync(ptr->d_bucketStarts, data->bucketStarts, startSize, cudaMemcpyHostToDevice, stream)); cudaChk(cudaMemcpyAsync(ptr->d_bucketSizes, data->bucketSizes, startSize, cudaMemcpyHostToDevice, stream)); -#ifdef CUDA_VERBOSE_KERNEL_ENQUEUE - printf("(%d) TRANSFER BASIC %zu bucket_markers %zu bucket_starts %zu\n", +#ifdef CUDA_VERBOSE_OBS + printf("(%d) PE TRANSFER BASIC %zu bucket_markers %zu bucket_starts %zu\n", CmiMyPe(), listSize, markerSize, @@ -605,11 +514,24 @@ void TreePieceDataTransferBasic(CudaRequest *data, CudaDevPtr *ptr){ /// @brief Free device memory used for interaction list and bucket data /// @param ptr CudaDevPtr object that stores handles to device memory -void TreePieceDataTransferBasicCleanup(CudaDevPtr *ptr){ - cudaChk(cudaFree(ptr->d_list)); - cudaChk(cudaFree(ptr->d_bucketMarkers)); - cudaChk(cudaFree(ptr->d_bucketStarts)); - cudaChk(cudaFree(ptr->d_bucketSizes)); +/// @param stream CUDA stream for async free operations +/// @param functionTag String literal identifying the calling function context. +void DataTransferBasicCleanup(CudaDevPtr *ptr, cudaStream_t stream, const char* functionTag){ + if (ptr == nullptr) return; + + // Only free non-NULL pointers to prevent illegal memory access + if (ptr->d_list != nullptr) { + cudaChk(gpuPoolFree(ptr->d_list, stream, functionTag)); + } + if (ptr->d_bucketMarkers != nullptr) { + cudaChk(gpuPoolFree(ptr->d_bucketMarkers, stream, functionTag)); + } + if (ptr->d_bucketStarts != nullptr) { + cudaChk(gpuPoolFree(ptr->d_bucketStarts, stream, functionTag)); + } + if (ptr->d_bucketSizes != nullptr) { + cudaChk(gpuPoolFree(ptr->d_bucketSizes, stream, functionTag)); + } } /** @brief Transfer forces from the GPU back to the host. Also schedules @@ -622,7 +544,6 @@ void TreePieceDataTransferBasicCleanup(CudaDevPtr *ptr){ */ void TransferParticleVarsBack(VariablePartData *hostBuffer, size_t size, void *d_varParts, cudaStream_t stream, void *cb){ - HAPI_TRACE_BEGIN(); cudaChk(cudaMemcpyAsync(hostBuffer, d_varParts, size, cudaMemcpyDeviceToHost, stream)); HAPI_TRACE_END(CUDA_XFER_BACK); @@ -718,7 +639,7 @@ __device__ __forceinline__ void stackPop(int &sp) { --sp; } -const int stackDepth = 64; +const int stackDepth = 128; //__launch_bounds__(1024,1) __global__ void gpuLocalTreeWalk( @@ -957,9 +878,7 @@ __global__ void gpuLocalTreeWalk( */ // 2d thread blocks -#ifdef CUDA_2D_TB_KERNEL #define TRANSLATE(x,y) (y*NODES_PER_BLOCK+x) -#ifndef CUDA_2D_FLAT __device__ __forceinline__ void ldg_moments(CudaMultipoleMoments &m, CudaMultipoleMoments *ptr) { m.radius = __ldg(&(ptr->radius)); @@ -1015,7 +934,6 @@ __global__ void nodeGravityComputation( int *bucketStarts, int *bucketSizes, cudatype fperiod){ - // __shared__ CudaVector3D acc[THREADS_PER_BLOCK]; // __shared__ cudatype pot[THREADS_PER_BLOCK]; // __shared__ cudatype idt2[THREADS_PER_BLOCK]; @@ -1204,337 +1122,6 @@ __global__ void nodeGravityComputation( }// end for each PARTICLE group } -#else -__global__ void nodeGravityComputation( - CompactPartData *particleCores, - VariablePartData *particleVars, - CudaMultipoleMoments *moments, - ILCell *ils, - int *ilmarks, - int *bucketStarts, - int *bucketSizes, - cudatype fperiod){ - - __shared__ cudatype accx[THREADS_PER_BLOCK]; - __shared__ cudatype accy[THREADS_PER_BLOCK]; - __shared__ cudatype accz[THREADS_PER_BLOCK]; - __shared__ cudatype pot[THREADS_PER_BLOCK]; - __shared__ cudatype idt2[THREADS_PER_BLOCK]; - //__shared__ cudatype mr[NODES_PER_BLOCK]; - __shared__ cudatype ms[NODES_PER_BLOCK]; - __shared__ cudatype mt[NODES_PER_BLOCK]; - __shared__ cudatype mcmx[NODES_PER_BLOCK]; - __shared__ cudatype mcmy[NODES_PER_BLOCK]; - __shared__ cudatype mcmz[NODES_PER_BLOCK]; - __shared__ cudatype mxx[NODES_PER_BLOCK]; - __shared__ cudatype mxy[NODES_PER_BLOCK]; - __shared__ cudatype mxz[NODES_PER_BLOCK]; - __shared__ cudatype myy[NODES_PER_BLOCK]; - __shared__ cudatype myz[NODES_PER_BLOCK]; - __shared__ cudatype mzz[NODES_PER_BLOCK]; - __shared__ int offsetID[NODES_PER_BLOCK]; - __shared__ CompactPartData shared_particle_cores[PARTS_PER_BLOCK]; - - int start = ilmarks[blockIdx.x]; - int end = ilmarks[blockIdx.x+1]; - int bucketStart = bucketStarts[blockIdx.x]; - int bucketSize = bucketSizes[blockIdx.x]; - - /* - __shared__ int start; - __shared__ int end; - __shared__ int bucketStart; - __shared__ int bucketSize; - */ - - int tx, ty; - -/* - if(threadIdx.x == 0 && threadIdx.y == 0){ - start = ilmarks[blockIdx.x]; - end = ilmarks[blockIdx.x+1]; - bucketStart = bucketStarts[blockIdx.x]; - bucketSize = bucketSizes[blockIdx.x]; - } - __syncthreads(); - */ - - int xstart; - int ystart; - tx = threadIdx.x; - ty = threadIdx.y; - - for(ystart = 0; ystart < bucketSize; ystart += PARTS_PER_BLOCK){ - - - int my_particle_idx = ystart + ty; - if(tx == 0 && my_particle_idx < bucketSize){ - shared_particle_cores[ty] = particleCores[bucketStart+my_particle_idx]; - } - - __syncthreads(); // wait for leader threads to finish using acc's, pot's of other threads - accx[TRANSLATE(tx,ty)] = 0.0; - accy[TRANSLATE(tx,ty)] = 0.0; - accz[TRANSLATE(tx,ty)] = 0.0; - pot[TRANSLATE(tx,ty)] = 0.0; - idt2[TRANSLATE(tx,ty)] = 0.0; - - - for(xstart = start; xstart < end; xstart += NODES_PER_BLOCK){ - int my_cell_idx = xstart + tx; - ILCell ilc; - - __syncthreads(); // wait for all threads to finish using - // previous iteration's nodes before reloading - - if(ty == 0 && my_cell_idx < end){ - ilc = ils[my_cell_idx]; - //mr[tx] = moments[ilc.index].radius; - ms[tx] = moments[ilc.index].soft; - mt[tx] = moments[ilc.index].totalMass; - mcmx[tx] = moments[ilc.index].cm.x; - mcmy[tx] = moments[ilc.index].cm.y; - mcmz[tx] = moments[ilc.index].cm.z; - mxx[tx] = moments[ilc.index].xx; - mxy[tx] = moments[ilc.index].xy; - mxz[tx] = moments[ilc.index].xz; - myy[tx] = moments[ilc.index].yy; - myz[tx] = moments[ilc.index].yz; - mzz[tx] = moments[ilc.index].zz; - offsetID[tx] = ilc.offsetID; - } - - __syncthreads(); // wait for nodes to be loaded before using them - - if(my_particle_idx < bucketSize && my_cell_idx < end){ // INTERACT - CudaVector3D r; - cudatype rsq; - cudatype twoh, a, b, c, d; - - r.x = shared_particle_cores[ty].position.x - - ((((offsetID[tx] >> 22) & 0x7)-3)*fperiod + mcmx[tx]); - r.y = shared_particle_cores[ty].position.y - - ((((offsetID[tx] >> 25) & 0x7)-3)*fperiod + mcmy[tx]); - r.z = shared_particle_cores[ty].position.z - - ((((offsetID[tx] >> 28) & 0x7)-3)*fperiod + mcmz[tx]); - - rsq = r.x*r.x + r.y*r.y + r.z*r.z; - twoh = ms[tx] + shared_particle_cores[ty].soft; - if(rsq != 0){ - cudatype dir = 1.0/sqrt(rsq); - // SPLINEQ(dir, rsq, twoh, a, b, c, d); - // expansion of function below: - cudatype u,dih; - if (rsq < twoh*twoh) { - dih = 2.0/twoh; - u = dih/dir; - if (u < 1.0) { - a = dih*(7.0/5.0 - 2.0/3.0*u*u + 3.0/10.0*u*u*u*u - - 1.0/10.0*u*u*u*u*u); - b = dih*dih*dih*(4.0/3.0 - 6.0/5.0*u*u + 1.0/2.0*u*u*u); - c = dih*dih*dih*dih*dih*(12.0/5.0 - 3.0/2.0*u); - d = 3.0/2.0*dih*dih*dih*dih*dih*dih*dir; - } - else { - a = -1.0/15.0*dir + dih*(8.0/5.0 - 4.0/3.0*u*u + u*u*u - - 3.0/10.0*u*u*u*u + 1.0/30.0*u*u*u*u*u); - b = -1.0/15.0*dir*dir*dir + dih*dih*dih*(8.0/3.0 - 3.0*u - + 6.0/5.0*u*u - 1.0/6.0*u*u*u); - c = -1.0/5.0*dir*dir*dir*dir*dir + 3.0*dih*dih*dih*dih*dir - + dih*dih*dih*dih*dih*(-12.0/5.0 + 1.0/2.0*u); - d = -dir*dir*dir*dir*dir*dir*dir - + 3.0*dih*dih*dih*dih*dir*dir*dir - - 1.0/2.0*dih*dih*dih*dih*dih*dih*dir; - } - } - else { - a = dir; - b = a*a*a; - c = 3.0*b*a*a; - d = 5.0*c*a*a; - } - - cudatype qirx = mxx[tx]*r.x + mxy[tx]*r.y + mxz[tx]*r.z; - cudatype qiry = mxy[tx]*r.x + myy[tx]*r.y + myz[tx]*r.z; - cudatype qirz = mxz[tx]*r.x + myz[tx]*r.y + mzz[tx]*r.z; - cudatype qir = 0.5*(qirx*r.x + qiry*r.y + qirz*r.z); - cudatype tr = 0.5*(mxx[tx] + myy[tx] + mzz[tx]); - cudatype qir3 = b*mt[tx] + d*qir - c*tr; - - pot[TRANSLATE(tx, ty)] -= mt[tx] * a + c*qir - b*tr; - - accx[TRANSLATE(tx, ty)] -= qir3*r.x - c*qirx; - accy[TRANSLATE(tx, ty)] -= qir3*r.y - c*qiry; - accz[TRANSLATE(tx, ty)] -= qir3*r.z - c*qirz; - idt2[TRANSLATE(tx, ty)] = fmax(idt2[TRANSLATE(tx, ty)], - (shared_particle_cores[ty].mass + mt[tx]) * b); - }// end if rsq != 0 - }// end INTERACT - }// end for each NODE group - - __syncthreads(); // wait for all threads to finish before results become available - - cudatype sumx, sumy, sumz, poten, idt2max; - sumx = sumy = sumz = poten = idt2max = 0.0; - // accumulate forces, potential in global memory data structure - if(tx == 0 && my_particle_idx < bucketSize){ - for(int i = 0; i < NODES_PER_BLOCK; i++){ - sumx += accx[TRANSLATE(i,ty)]; - sumy += accy[TRANSLATE(i,ty)]; - sumz += accz[TRANSLATE(i,ty)]; - poten += pot[TRANSLATE(i,ty)]; - idt2max = fmax(idt2[TRANSLATE(i,ty)], idt2max); - } - particleVars[bucketStart+my_particle_idx].a.x += sumx; - particleVars[bucketStart+my_particle_idx].a.y += sumy; - particleVars[bucketStart+my_particle_idx].a.z += sumz; - particleVars[bucketStart+my_particle_idx].potential += poten; - particleVars[bucketStart+my_particle_idx].dtGrav = fmax(idt2max, particleVars[bucketStart+my_particle_idx].dtGrav); - } - - }// end for each PARTICLE group -} -#endif -#else -__global__ void nodeGravityComputation( - CompactPartData *particleCores, - VariablePartData *particleVars, - CudaMultipoleMoments *moments, - ILCell *ils, - int *ilmarks, - int *bucketStarts, - int *bucketSizes, - cudatype fperiod){ - - // each thread has its own storage for these - __shared__ CudaVector3D acc[THREADS_PER_BLOCK]; - __shared__ cudatype pot[THREADS_PER_BLOCK]; - __shared__ cudatype idt2[THREADS_PER_BLOCK]; - __shared__ CudaMultipoleMoments m[THREADS_PER_BLOCK]; - - __shared__ CompactPartData shared_particle_core; - - - // each block is given a bucket to compute - // each thread in the block computes an interaction of a particle with a node - // threads must iterate through the interaction lists and sync. - // then, block leader (first rank in each block) reduces the forces and commits - // values to global memory. - int bucket = blockIdx.x; - int start = ilmarks[bucket]; - int end = ilmarks[bucket+1]; - int bucketSize = bucketSizes[bucket]; - int bucketStart = bucketStarts[bucket]; - int thread = threadIdx.x; - - CudaVector3D r; - cudatype rsq; - cudatype twoh, a, b, c, d; - - for(int particle = 0; particle < bucketSize; particle++){ - if(thread == 0){ - // load shared_particle_core - shared_particle_core = particleCores[bucketStart+particle]; - } - __syncthreads(); - - acc[thread].x = 0; - acc[thread].y = 0; - acc[thread].z = 0; - pot[thread] = 0; - idt2[thread] = 0; - - for(int node = start+thread; node < end; node+=THREADS_PER_BLOCK){ - ILCell ilc = ils[node]; - m[thread] = moments[ilc.index]; - int offsetID = ilc.offsetID; - - r.x = shared_particle_core.position.x - - ((((offsetID >> 22) & 0x7)-3)*fperiod + m[thread].cm.x); - r.y = shared_particle_core.position.y - - ((((offsetID >> 25) & 0x7)-3)*fperiod + m[thread].cm.y); - r.z = shared_particle_core.position.z - - ((((offsetID >> 28) & 0x7)-3)*fperiod + m[thread].cm.z); - - rsq = r.x*r.x + r.y*r.y + r.z*r.z; - twoh = m[thread].soft + shared_particle_core.soft; - if(rsq != 0){ - cudatype dir = 1.0/sqrt(rsq); - // SPLINEQ(dir, rsq, twoh, a, b, c, d); - // expansion of function below: - cudatype u,dih; - if (rsq < twoh*twoh) { - dih = 2.0/twoh; - u = dih/dir; - if (u < 1.0) { - a = dih*(7.0/5.0 - 2.0/3.0*u*u + 3.0/10.0*u*u*u*u - - 1.0/10.0*u*u*u*u*u); - b = dih*dih*dih*(4.0/3.0 - 6.0/5.0*u*u + 1.0/2.0*u*u*u); - c = dih*dih*dih*dih*dih*(12.0/5.0 - 3.0/2.0*u); - d = 3.0/2.0*dih*dih*dih*dih*dih*dih*dir; - } - else { - a = -1.0/15.0*dir + dih*(8.0/5.0 - 4.0/3.0*u*u + u*u*u - - 3.0/10.0*u*u*u*u + 1.0/30.0*u*u*u*u*u); - b = -1.0/15.0*dir*dir*dir + dih*dih*dih*(8.0/3.0 - 3.0*u - + 6.0/5.0*u*u - 1.0/6.0*u*u*u); - c = -1.0/5.0*dir*dir*dir*dir*dir + 3.0*dih*dih*dih*dih*dir - + dih*dih*dih*dih*dih*(-12.0/5.0 + 1.0/2.0*u); - d = -dir*dir*dir*dir*dir*dir*dir - + 3.0*dih*dih*dih*dih*dir*dir*dir - - 1.0/2.0*dih*dih*dih*dih*dih*dih*dir; - } - } - else { - a = dir; - b = a*a*a; - c = 3.0*b*a*a; - d = 5.0*c*a*a; - } - - cudatype qirx = m[thread].xx*r.x + m[thread].xy*r.y + m[thread].xz*r.z; - cudatype qiry = m[thread].xy*r.x + m[thread].yy*r.y + m[thread].yz*r.z; - cudatype qirz = m[thread].xz*r.x + m[thread].yz*r.y + m[thread].zz*r.z; - cudatype qir = 0.5*(qirx*r.x + qiry*r.y + qirz*r.z); - cudatype tr = 0.5*(m[thread].xx + m[thread].yy + m[thread].zz); - cudatype qir3 = b*m[thread].totalMass + d*qir - c*tr; - - pot[thread] -= m[thread].totalMass * a + c*qir - b*tr; - - acc[thread].x -= qir3*r.x - c*qirx; - acc[thread].y -= qir3*r.y - c*qiry; - acc[thread].z -= qir3*r.z - c*qirz; - idt2[thread] = fmax(idt2[thread], (shared_particle_core.mass + m[thread].totalMass)*b); - }// end if rsq != 0 - }// for each node in list - __syncthreads(); - // at this point, the total force on particle is distributed among - // all active threads; - // reduce. - // TODO: make this a parallel reduction - - cudatype sumx, sumy, sumz, poten, idt2max; - sumx = sumy = sumz = poten = idt2max = 0.0; - if(thread == 0){ - for(int i = 0; i < THREADS_PER_BLOCK; i++){ - sumx += acc[i].x; - sumy += acc[i].y; - sumz += acc[i].z; - poten += pot[i]; - idt2max = fmax(idt2[i], idt2max); - } - particleVars[bucketStart+particle].a.x += sumx; - particleVars[bucketStart+particle].a.y += sumy; - particleVars[bucketStart+particle].a.z += sumz; - particleVars[bucketStart+particle].potential += poten; - particleVars[bucketStart+particle].dtGrav = fmax(idt2max, particleVars[bucketStart+particle].dtGrav); - } - }// for each particle in bucket - - -} -#endif - /** * @brief interaction between source particles and buckets of particles. * @param particleCores Read-only properties of target particles. @@ -1556,7 +1143,6 @@ __device__ __forceinline__ void ldg_cPartData(CompactPartData &m, CompactPartDat m.position.z = __ldg(&(ptr->position.z)); } -#ifdef CUDA_2D_TB_KERNEL #define TRANSLATE_PART(x,y) (y*NODES_PER_BLOCK_PART+x) //__launch_bounds__(maxThreadsPerBlock, minBlocksPerMultiprocessor) //maxThreadsPerBlock needs to be a multiple of 128, this can be used @@ -1572,7 +1158,7 @@ __global__ void particleGravityComputation( int *bucketStarts, int *bucketSizes, cudatype fperiod){ - + //__shared__ CudaVector3D acc[THREADS_PER_BLOCK_PART]; // __shared__ cudatype pot[THREADS_PER_BLOCK_PART]; // __shared__ cudatype idt2[THREADS_PER_BLOCK_PART]; @@ -1749,228 +1335,52 @@ __global__ void particleGravityComputation( }// end for each PARTICLE group } -#else -__launch_bounds__(896, 1) -__global__ void particleGravityComputation( - CompactPartData *targetCores, - VariablePartData *targetVars, - CompactPartData *sourceCores, - ILPart *ils, - int *ilmarks, - int *bucketStarts, - int *bucketSizes, - cudatype fperiod){ - - // each thread has its own storage for these - __shared__ CudaVector3D acc[THREADS_PER_BLOCK]; - __shared__ cudatype pot[THREADS_PER_BLOCK]; - __shared__ cudatype idt2[THREADS_PER_BLOCK]; - __shared__ CompactPartData source_cores[THREADS_PER_BLOCK]; - - __shared__ CompactPartData shared_target_core; - - - // each block is given a bucket to compute - // each thread in the block computes an interaction of a particle with a node - // threads must iterate through the interaction lists and sync. - // then, block leader (first rank in each block) reduces the forces and commits - // values to global memory. - int bucket = blockIdx.x; - int start = ilmarks[bucket]; - int end = ilmarks[bucket+1]; - int bucketSize = bucketSizes[bucket]; - int bucketStart = bucketStarts[bucket]; - int thread = threadIdx.x; - - CudaVector3D r; - cudatype rsq; - cudatype twoh, a, b; - - for(int target = 0; target < bucketSize; target++){ - if(thread == 0){ - shared_target_core = targetCores[bucketStart+target]; - } - __syncthreads(); - - acc[thread].x = 0; - acc[thread].y = 0; - acc[thread].z = 0; - pot[thread] = 0; - idt2[thread] = 0; - - for(int source = start+thread; source < end; source += THREADS_PER_BLOCK){ - ILPart ilp = ils[source]; - int oid = ilp.off; - int num = ilp.num; - int ilpindex = ilp.index; - - for(int particle = 0; particle < num; particle++){ - source_cores[thread] = sourceCores[ilpindex+particle]; - - r.x = (((oid >> 22) & 0x7)-3)*fperiod + - source_cores[thread].position.x - - shared_target_core.position.x; - - r.y = (((oid >> 25) & 0x7)-3)*fperiod + - source_cores[thread].position.y - - shared_target_core.position.y; - - r.z = (((oid >> 28) & 0x7)-3)*fperiod + - source_cores[thread].position.z - - shared_target_core.position.z; - - rsq = r.x*r.x + r.y*r.y + r.z*r.z; - twoh = source_cores[thread].soft + shared_target_core.soft; - if(rsq != 0){ - cudatype r1, u,dih,dir; - r1 = sqrt(rsq); - if (r1 < (twoh)) { - dih = 2.0/(twoh); - u = r1*dih; - if (u < 1.0) { - a = dih*(7.0/5.0 - 2.0/3.0*u*u + 3.0/10.0*u*u*u*u - - 1.0/10.0*u*u*u*u*u); - b = dih*dih*dih*(4.0/3.0 - 6.0/5.0*u*u - + 1.0/2.0*u*u*u); - } - else { - dir = 1.0/r1; - a = -1.0/15.0*dir + dih*(8.0/5.0 - 4.0/3.0*u*u + - u*u*u - 3.0/10.0*u*u*u*u + 1.0/30.0*u*u*u*u*u); - b = -1.0/15.0*dir*dir*dir + - dih*dih*dih*(8.0/3.0 - 3.0*u + - 6.0/5.0*u*u - 1.0/6.0*u*u*u); - } - } - else { - a = 1.0/r1; - b = a*a*a; - } - - pot[thread] -= source_cores[thread].mass * a; - - acc[thread].x += r.x*b*source_cores[thread].mass; - acc[thread].y += r.y*b*source_cores[thread].mass; - acc[thread].z += r.z*b*source_cores[thread].mass; - idt2[thread] = fmax(idt2[thread], (shared_target_core.mass + source_cores[thread].mass) * b); - }// if rsq != 0 - }// for each particle in source bucket - }// for each source bucket - - __syncthreads(); - - cudatype sumx, sumy, sumz, poten, idt2max; - sumx = sumy = sumz = poten = idt2max = 0.0; - if(thread == 0){ - for(int i = 0; i < THREADS_PER_BLOCK; i++){ - sumx += acc[i].x; - sumy += acc[i].y; - sumz += acc[i].z; - poten += pot[i]; - idt2max = fmax(idt2[i], idt2max); - } - targetVars[bucketStart+target].a.x += sumx; - targetVars[bucketStart+target].a.y += sumy; - targetVars[bucketStart+target].a.z += sumz; - targetVars[bucketStart+target].potential += poten; - targetVars[bucketStart+target].dtGrav = fmax(idt2max, targetVars[bucketStart+target].dtGrav); - } - - }// for each target part -} -#endif - -__global__ void EwaldKernel(CompactPartData *particleCores, VariablePartData *particleVars, int *markers, int largephase, int First, int Last); extern unsigned int timerHandle; -void EwaldHostMemorySetup(EwaldData *h_idata, int nParticles, int nEwhLoop, int largephase) { - if(largephase) - allocatePinnedHostMemory((void **)&(h_idata->EwaldMarkers), nParticles*sizeof(int)); - else - h_idata->EwaldMarkers = NULL; - allocatePinnedHostMemory((void **)&(h_idata->ewt), nEwhLoop*sizeof(EwtData)); - allocatePinnedHostMemory((void **)&(h_idata->cachedData), sizeof(EwaldReadOnlyData)); -} - -void EwaldHostMemoryFree(EwaldData *h_idata, int largephase) { - if(largephase) - freePinnedHostMemory(h_idata->EwaldMarkers); - freePinnedHostMemory(h_idata->ewt); - freePinnedHostMemory(h_idata->cachedData); -} - -/** @brief Set up CUDA kernels to perform Ewald sum. - * @param d_localParts Local particle data on device - * @param d_localVars Local particle accelerations on device - * @param h_idata Host data buffers - * @param stream CUDA stream to perform GPU operations over - * @param cb Callback - * @param myIndex Chare index on this node that called this request. - * @param largephase Whether to perform large or small phase calculation - * - * The "top" and "bottom" Ewlad kernels have been combined: - * "top" for the real space loop, - * "bottom" for the k-space loop. - * - */ -void EwaldHost(CompactPartData *d_localParts, VariablePartData *d_localVars, - EwaldData *h_idata, cudaStream_t stream, void *cb, int myIndex, int largephase) -{ - int n = h_idata->cachedData->n; - int numBlocks = (int) ceilf((float)n/BLOCK_SIZE); - int nEwhLoop = h_idata->cachedData->nEwhLoop; - assert(nEwhLoop <= NEWH); - - size_t size; - if(largephase) size = n * sizeof(int); - else size = 0; +/// @brief Launch the Ewald calculation on the GPU from the DataManager +/// @param d_localParts Pointer to read-only particle data on the GPU +/// @param d_localVars Pointer to modifiable particle data on the GPU +/// @param _ewt Pointer to EwaldReadOnlyData on the host +/// @param _cachedData Pointer to cachedData on the host +/// @param numParts Number of particles to do the calculation for +/// @param stream The CUDA stream to launch the calculation on +/// @param cb Callback function after the kernel finishes +void DataManagerEwald(void *d_localParts, void *d_localVars, void *_ewt, void *_cachedData, int numParts, cudaStream_t stream, void *cb) { + // The following rounds up. + int numBlocks = (numParts + (THREADS_PER_BLOCK - 1)) / THREADS_PER_BLOCK; + +#ifdef CUDA_VERBOSE_OPS + CkPrintf("(%d) DM EWALD numParts: %d numBlocks: %d\n", + CmiMyNode(), + numParts, numBlocks + ); +#endif HAPI_TRACE_BEGIN(); - int *d_EwaldMarkers; - cudaChk(cudaMalloc(&d_EwaldMarkers, size)); - - cudaMemcpyAsync(d_EwaldMarkers, h_idata->EwaldMarkers, size, cudaMemcpyHostToDevice, stream); - cudaMemcpyToSymbolAsync(cachedData, h_idata->cachedData, sizeof(EwaldReadOnlyData), 0, cudaMemcpyHostToDevice, stream); - cudaMemcpyToSymbolAsync(ewt, h_idata->ewt, nEwhLoop * sizeof(EwtData), 0, cudaMemcpyHostToDevice, stream); - - if (largephase) - EwaldKernel<<>>(d_localParts, - d_localVars, - d_EwaldMarkers, 1, - h_idata->EwaldRange[0], h_idata->EwaldRange[1]); - else - EwaldKernel<<>>(d_localParts, - d_localVars, - NULL, 0, - h_idata->EwaldRange[0], h_idata->EwaldRange[1]); - HAPI_TRACE_END(CUDA_EWALD); -#ifdef CUDA_VERBOSE_KERNEL_ENQUEUE - printf("[%d] in EwaldHost, enqueued EwaldKernel\n", myIndex); + cudaMemcpyToSymbolAsync(ewt, _ewt, ((EwaldReadOnlyData *)_cachedData)->nEwhLoop * sizeof(EwtData), 0, cudaMemcpyHostToDevice, stream); + cudaMemcpyToSymbolAsync(cachedData, _cachedData, sizeof(EwaldReadOnlyData), 0, cudaMemcpyHostToDevice, stream); + +#ifndef CUDA_NO_KERNELS + EwaldKernel<<>>((CompactPartData *)d_localParts, + (VariablePartData *)d_localVars, + numParts); #endif + HAPI_TRACE_END(CUDA_EWALD); - cudaChk(cudaPeekAtLastError()); hapiAddCallback(stream, cb); - cudaChk(cudaFree(d_EwaldMarkers)); } __global__ void EwaldKernel(CompactPartData *particleCores, VariablePartData *particleVars, - int *markers, int largephase, - int First, int Last) { + int nParts) { ///////////////////////////////////// ////////////// Ewald TOP //////////// ///////////////////////////////////// int id; - if(largephase){ - id = blockIdx.x * BLOCK_SIZE + threadIdx.x; - if(id > Last) return; - id = markers[id]; - }else{ - id = First + blockIdx.x * BLOCK_SIZE + threadIdx.x; - if(id > Last) return; - } + id = blockIdx.x * THREADS_PER_BLOCK + threadIdx.x; + if(id >= nParts) return; CompactPartData *p; @@ -2194,7 +1604,7 @@ __global__ void EwaldKernel(CompactPartData *particleCores, // initialize accelerations and potentials to zero __global__ void ZeroVars(VariablePartData *particleVars, int nVars) { int id; - id = blockIdx.x * BLOCK_SIZE + threadIdx.x; + id = blockIdx.x * THREADS_PER_BLOCK + threadIdx.x; if(id >= nVars) return; particleVars[id].a.x = 0.0; diff --git a/HostCUDA.h b/HostCUDA.h index f5246b07..4df16a72 100644 --- a/HostCUDA.h +++ b/HostCUDA.h @@ -13,14 +13,12 @@ #define WARP_INDEX (threadIdx.x >> 5) #endif //GPU_LOCAL_TREE_WALK -#ifdef CUDA_2D_TB_KERNEL #define PARTS_PER_BLOCK 16 #define NODES_PER_BLOCK (THREADS_PER_BLOCK/PARTS_PER_BLOCK) #define THREADS_PER_BLOCK_PART 128 #define PARTS_PER_BLOCK_PART 16 #define NODES_PER_BLOCK_PART (THREADS_PER_BLOCK_PART/PARTS_PER_BLOCK_PART) -#endif // FIXME - find appropriate values #define NUM_INIT_MOMENT_INTERACTIONS_PER_BUCKET 100 @@ -29,8 +27,7 @@ /** @brief Data and parameters for requesting gravity calculations on * the GPU. */ typedef struct _CudaRequest{ - /// CUDA stream to handle memory operations and kernel launches for this request - /// Allocation and deallocation of the stream is handled by the DataManager + /// Handle to the stream where this request should be sent cudaStream_t stream; /// for accessing device memory @@ -101,7 +98,7 @@ void freePinnedHostMemory(void *); void DataManagerTransferLocalTree(void *moments, size_t sMoments, void *compactParts, size_t sCompactParts, - void *varParts, size_t sVarParts, + size_t sVarParts, void **d_localMoments, void **d_compactParts, void **d_varParts, cudaStream_t stream, int numParticles, void *callback); @@ -111,17 +108,16 @@ void DataManagerTransferRemoteChunk(void *moments, size_t sMoments, cudaStream_t stream, void *callback); -void TransferParticleVarsBack(VariablePartData *hostBuffer, size_t size, void *d_varParts, cudaStream_t stream, void *cb); - -void TreePieceCellListDataTransferLocal(CudaRequest *data); -void TreePieceCellListDataTransferRemote(CudaRequest *data); -void TreePieceCellListDataTransferRemoteResume(CudaRequest *data); +void DataManagerLocalTreeWalk(CudaRequest *data); +void PEListPartListDataTransferLocal(CudaRequest *data); +void PEListNodeListDataTransferLocal(CudaRequest *data); +void PEListNodeListDataTransferRemote(CudaRequest *data); +void PEListNodeListDataTransferRemoteResume(CudaRequest *data); +void PEListPartListDataTransferRemote(CudaRequest *data); +void PEListPartListDataTransferRemoteResume(CudaRequest *data); +void DataManagerEwald(void *d_localParts, void *d_localVars, void *_ewt, void *_cachedData, int numParts, cudaStream_t stream, void *cb); - -void TreePiecePartListDataTransferLocal(CudaRequest *data); -void TreePiecePartListDataTransferLocalSmallPhase(CudaRequest *data, CompactPartData *parts, int len); -void TreePiecePartListDataTransferRemote(CudaRequest *data); -void TreePiecePartListDataTransferRemoteResume(CudaRequest *data); +void TransferParticleVarsBack(VariablePartData *hostBuffer, size_t size, void *d_varParts, cudaStream_t stream, void *cb); void DummyKernel(void *cb); diff --git a/MEMORY_POOL_CHANGES.txt b/MEMORY_POOL_CHANGES.txt new file mode 100644 index 00000000..082bb2f6 --- /dev/null +++ b/MEMORY_POOL_CHANGES.txt @@ -0,0 +1,174 @@ +Major overhaul of GPU and host memory pool infrastructure + +This commit implements comprehensive enhancements to ChaNGa's memory pool +system, introducing smart memory management with trimming, adaptive refill, +detailed analytics, and dual allocation APIs (pool-based and direct). + +================================================================================ +GPU POOL ENHANCEMENTS +================================================================================ + +API Reorganization: +- Renamed all GPU pool functions from poolXxx to gpuPoolXxx for clarity +- Split allocation APIs into pool-based (gpuPoolMalloc/Free) and direct + (gpuMalloc/Free) variants +- Direct APIs bypass pool but maintain logging and analytics +- Updated macros: gpuMallocHelper → gpuPoolMalloc, added gpuMalloc for direct + +Configuration & Warmup: +- Changed GPU pool limit from unlimited to 5GB to prevent excessive memory + usage when multiple ranks share a GPU device +- Tuned warmup allocation counts (reduced from ultra-aggressive to moderate) +- Warmup now allocates 64KB-1GB range with balanced counts per bucket + +Implementation: +- Added logging helpers (logGpuAlloc, logGpuFree) in anonymous namespace +- Separated raw operations (gpuPoolMallocRaw/FreeRaw) from logged APIs +- Implemented type-safe template wrappers for both pool and direct allocation +- Maintained stream-aware allocation support (cudaMallocAsync/cudaFreeAsync) + +================================================================================ +HOST POOL MAJOR ENHANCEMENTS +================================================================================ + +Smart Trimming with Decay Rate: +- Implemented hostPoolTrim() with capacity-based intelligent trimming +- Trim prioritizes large buckets first while respecting usage-based protection +- Recent usage tracked per bucket (recentHits counter) +- Decay mechanism: recentHits decrements each timestep, allowing inactive + buckets to become trimmable while active buckets maintain protection +- Configurable target capacity (GB) and minimum per-bucket capacity (MB) + +Adaptive Refill Mechanism: +- Implemented hostPoolAdaptiveRefill() to maintain pool health +- Tracks misses per bucket during timestep +- Automatically refills under-provisioned buckets between timesteps +- Refills up to 3 blocks for hot buckets that had misses +- Prevents blocking on frequently-used allocation sizes + +Bucket Granularity Improvements: +- Added fine-grained tiny buckets: 256B, 512B (previously started at 1KB) +- Small allocations (≤64KB): Fine granularity (256B, 512B, 1KB, 2KB, 4KB, ...) +- Medium allocations (64KB-1MB): 64KB granularity +- Medium-large (1MB-64MB): 1MB granularity +- Large allocations (>64MB): 50MB granularity to prevent bucket proliferation +- Drastically reduces unique bucket count for large allocations + +Warmup Tuning: +- Redesigned warmup with 19 bucket sizes from 256B to 64MB +- Conservative warmup counts based on empirical profiling +- High warmup for tiny/small (1000-1500 blocks for ≤8KB) +- Moderate for medium (150-250 blocks for 16KB-128KB) +- Minimal for large (1 block for ≥1MB) +- Warmup targets used by trim logic to preserve critical buckets + +Comprehensive Analytics: +- Per-bucket tracking: total allocations, cumulative misses, recent usage +- Cumulative metrics: total allocated/freed, peak memory, current memory +- Per-timestep metrics: allocs, frees, hits, misses, blocks freed to OS +- Historical timestep snapshots for growth trend analysis +- Pool accounting validation (owned blocks = free + live) + +Detailed Reporting & Analysis: +- hostPoolReportStats(): Comprehensive pool status with per-bucket breakdown + Shows: capacity usage, hit rates, memory usage, per-bucket distribution +- hostPoolAnalyzeWarmupEffectiveness(): Compares warmup targets vs actual usage + Identifies over-provisioned and under-provisioned buckets +- hostPoolAnalyzeGrowth(): Analyzes memory growth trends over time + Detects concerning patterns (bucket proliferation, memory leaks) +- All reports include formatted output with bucket sizes (B/KB/MB) + +Direct Allocation API: +- Added hostMalloc/hostFree that bypass pool but maintain logging +- Useful for large, long-lived buffers that shouldn't consume pool capacity +- Consistent logging interface with pool-based allocations +- Updated local tree buffer allocations to use direct API + +API Standardization: +- Renamed macros: hostMallocHelper → hostPoolMalloc, hostFreeHelper → hostPoolFree +- Added hostMalloc/hostFree macros for direct allocation +- Separate implementation functions (*Impl) vs public macros +- Type-safe template wrappers for pool and direct variants +- Consistent parameter ordering across all APIs + +================================================================================ +DATAMANAGER INTEGRATION +================================================================================ + +New Configuration Parameters: +- dHostPoolTargetCapacityGB: Target pool capacity in GB (default: 5.0) +- nHostPoolMinCapacityPerBucketMB: Min capacity per bucket during trim (default: 50MB) +- bHostPoolDebug: Enable detailed pool diagnostics and reporting (default: 0) + +New Entry Methods: +- trimHostPool(): Collective trim operation across all DataManagers +- refillHostPool(): Collective adaptive refill across all DataManagers + +Local Tree Buffer Updates: +- Changed to direct allocation (hostMalloc) for large local tree buffers +- bufLocalMoments, bufLocalParts, bufLocalVars use direct API +- Reduces pool pressure from large, long-lived allocations + +Bug Fixes: +- Fixed typo: treePiecesBuff → treePiecesBufferFilled = 0 +- Added cudaDeviceSynchronize() after particle update cleanup + +================================================================================ +INTEGRATION & USAGE UPDATES +================================================================================ + +HostCUDA.cu Changes: +- Updated allocatePinnedHostMemory/freePinnedHostMemory to use pool API +- Remote chunk transfers use gpuPoolMalloc (transient allocations) +- Local tree transfers use gpuMalloc (long-lived allocations) +- Interaction list buffers use gpuPoolMalloc (per-bucket allocations) +- Missed nodes/parts use gpuPoolMalloc (resume computation buffers) + +ParallelGravity Integration: +- Added command-line parameters: hptarget, hpminbucket, hpdebug +- Integrated pool maintenance in Main::doSimulation() after each BigStep: + 1. Optional diagnostics (if bHostPoolDebug): report stats, analyze growth + 2. Trim host pool to target capacity + 3. Refill hot buckets adaptively + 4. Flush memory logs +- Warmup effectiveness analysis at step 5 (if debug enabled) +- Growth analysis every 10 steps (if debug enabled) + +Memory Log Infrastructure: +- Enhanced GPU and CPU logging with detailed event tracking +- Logs capture alloc/free operations with source location +- Support for MEMLOG_FREE_SKIP, MEMLOG_ALLOC_FAIL, MEMLOG_FREE_FAIL +- Sequential flushing across nodes to avoid file corruption + +================================================================================ +INITIALIZATION UPDATES +================================================================================ + +- Renamed poolInit → gpuPoolInit, poolInitWithWarmup → gpuPoolInit (integrated) +- Renamed hostMemoryPoolInit → hostPoolInit +- Renamed poolDestroy → gpuPoolDestroy +- Renamed hostMemoryPoolCleanup → hostPoolCleanup +- Updated ProjectionsControl to use new initialization functions +- Static initialization ensures pools initialized once per node + +================================================================================ +TESTING & VALIDATION +================================================================================ + +Pool Accounting: +- Validates total owned blocks = free pool blocks + live (checked out) blocks +- Warns on accounting mismatch to detect potential leaks +- Tracks per-bucket usage percentages + +Performance Tracking: +- Total allocation time tracking (ms average) +- Per-timestep allocation/free counts +- Hit rate monitoring (% of allocations served from pool) +- Timestep miss tracking for adaptive refill + +Growth Detection: +- Historical snapshots track: buckets, free blocks, bytes, live blocks +- Detects bucket proliferation (>10 new buckets) +- Detects memory growth (>1GB increase) +- Detects live block accumulation (potential leak indicator) + diff --git a/MemoryPool.cpp b/MemoryPool.cpp new file mode 100644 index 00000000..150400aa --- /dev/null +++ b/MemoryPool.cpp @@ -0,0 +1,1763 @@ +#ifdef CUDA + +#include +#include "MemoryPool.h" +#include "ParallelGravity.h" +#include "DataManager.h" +#include "memlog.h" +#include "ParallelGravity.decl.h" +#include "formatted_string.h" +#include +#include + + +// ============================================================================ +// GPU MEMORY POOL +// ============================================================================ + + +// ---------------------------------------------------------------------------- +// GPU Pool Initialization & Warmup +// ---------------------------------------------------------------------------- + +/** + * @brief Initialize the GPU memory pool + * + * Configures the default device memory pool with a 5GB limit to prevent + * excessive memory consumption when multiple ranks share a GPU device. + */ +void gpuPoolInit() { + const int maxCudaDevice = 16; + int nCudaDeviceCount = 0; + cudaGetDeviceCount(&nCudaDeviceCount); + CkAssert(nCudaDeviceCount <= maxCudaDevice); + + int dev = 0; + cudaGetDevice(&dev); + + // Use PE-based coordination: only the lowest-numbered PE initializes each device + // This ensures one initialization per device across all processes on the node + static int initializedDevices[maxCudaDevice] = {0}; + + if (initializedDevices[dev]) { + return; + } + + // Atomically mark this device as being initialized by this PE + initializedDevices[dev] = 1; + + // Get the CUDA memory pool for this device + cudaMemPool_t pool = nullptr; + cudaDeviceGetDefaultMemPool(&pool, dev); + + // Set 5GB limit (5 * 1024^3 bytes) to prevent excessive memory usage + // with multiple ranks per GPU device + unsigned long long threshold = 5ULL * 1024 * 1024 * 1024; + cudaError_t setResult = cudaMemPoolSetAttribute(pool, cudaMemPoolAttrReleaseThreshold, &threshold); + if (setResult != cudaSuccess) { + CkAbort("gpuPoolInit: cudaMemPoolSetAttribute failed: %s", cudaGetErrorString(setResult)); + } + + // Allocation sizes based on typical ChaNGa workload patterns + // Extended range from 64 KB to 1 GB for comprehensive coverage + size_t sizes[] = { + 64 * 1024, // 64 KB + 128 * 1024, // 128 KB + 256 * 1024, // 256 KB + 512 * 1024, // 512 KB + 1 * 1024 * 1024, // 1 MB + 2 * 1024 * 1024, // 2 MB + 4 * 1024 * 1024, // 4 MB + 8 * 1024 * 1024, // 8 MB + 16 * 1024 * 1024, // 16 MB + 32 * 1024 * 1024, // 32 MB + 64 * 1024 * 1024, // 64 MB + 128 * 1024 * 1024, // 128 MB + 256 * 1024 * 1024, // 256 MB + 512 * 1024 * 1024, // 512 MB + 1024 * 1024 * 1024 // 1 GB + }; + + // Total warmup will pre-allocate and free significant memory + // to fully populate the pool's free-lists + int counts[] = { + 80, // 64 KB + 80, // 128 KB + 60, // 256 KB + 60, // 512 KB + 60, // 1 MB + 60, // 2 MB + 40, // 4 MB + 40, // 8 MB + 30, // 16 MB + 20, // 32 MB + 10, // 64 MB + 5, // 128 MB + 5, // 256 MB + 2, // 512 MB + 1 // 1 GB + }; + + int numSizes = sizeof(sizes) / sizeof(sizes[0]); + + cudaError_t result = gpuPoolWarmup(sizes, counts, numSizes, 0); + if (result == cudaSuccess) { + // Only print from one nodegroup instance to avoid spam in production + if (CkMyNode() == 0) { + CkPrintf("GPU memory pool warmed up\n"); + } + } else { + // Print warnings from all instances so we don't miss failures + CkPrintf("WARNING: GPU memory pool warmup failed: %s\n", cudaGetErrorString(result)); + } +} + +/** + * @brief Pre-warm the GPU memory pool + * + * Allocates and immediately frees blocks to populate the pool's free-list. + * This avoids allocation overhead during the critical computation path. + * + * @param sizes Array of allocation sizes in bytes + * @param counts Array of how many allocations per size + * @param numSizes Number of different sizes + * @param stream CUDA stream to use (0 for default) + * @return cudaError_t from allocation operations, or cudaSuccess + */ +cudaError_t gpuPoolWarmup(const size_t* sizes, const int* counts, int numSizes, cudaStream_t stream) { + if (sizes == nullptr || counts == nullptr || numSizes <= 0) { + return cudaErrorInvalidValue; + } + + // Use default stream if none provided + bool createdStream = false; + if (stream == nullptr) { + cudaError_t createResult = cudaStreamCreate(&stream); + if (createResult != cudaSuccess) { + return createResult; + } + createdStream = true; + } + + // Allocate and free blocks to populate the pool + cudaError_t finalResult = cudaSuccess; + for (int i = 0; i < numSizes; i++) { + size_t size = sizes[i]; + int count = counts[i]; + + for (int j = 0; j < count; j++) { + void* ptr = nullptr; + + // Allocate block + cudaError_t result = cudaMallocAsync(&ptr, size, stream); + if (result != cudaSuccess) { + finalResult = result; + break; + } + + // Immediately free it back to the pool + cudaFreeAsync(ptr, stream); + } + + if (finalResult != cudaSuccess) { + break; + } + } + + // Synchronize to ensure all warmup operations complete + if (finalResult == cudaSuccess) { + finalResult = cudaStreamSynchronize(stream); + } + + // Clean up stream if we created it + if (createdStream) { + cudaStreamDestroy(stream); + } + + return finalResult; +} + +// ---------------------------------------------------------------------------- +// GPU Pool Raw Operations (Internal - use macros for public API) +// ---------------------------------------------------------------------------- + +/** + * @brief Allocate GPU memory from the pool (raw, internal, stream-aware) + * + * Uses cudaMallocAsync for stream-ordered allocation when a stream is provided, + * falling back to synchronous cudaMalloc when stream is NULL. + * + * @param devPtr Output pointer to allocated device memory + * @param size Size of memory to allocate in bytes + * @param stream CUDA stream for ordering (NULL for synchronous) + * @return cudaError_t result from allocation + */ +cudaError_t gpuPoolMallocRaw(void** devPtr, size_t size, cudaStream_t stream) { + if (devPtr == nullptr) { + return cudaErrorInvalidValue; + } + + // Initialize to NULL in case of failure + *devPtr = nullptr; + + if (size == 0) { + return cudaSuccess; + } + + cudaError_t result; + if (stream != nullptr) { + // Stream-ordered async allocation + result = cudaMallocAsync(devPtr, size, stream); + } else { + // Fallback to synchronous allocation + result = cudaMalloc(devPtr, size); + } + + // Ensure pointer is NULL on failure + if (result != cudaSuccess) { + *devPtr = nullptr; + } + + return result; +} + +/** + * @brief Free GPU memory back to the pool (raw, internal, stream-aware) + * + * Uses cudaFreeAsync for stream-ordered deallocation when a stream is provided, + * falling back to synchronous cudaFree when stream is NULL. + * + * @param ptr Device pointer to free + * @param stream CUDA stream for ordering (NULL for synchronous) + * @return cudaError_t result from free operation + */ +cudaError_t gpuPoolFreeRaw(void* ptr, cudaStream_t stream) { + if (ptr == nullptr) { + return cudaSuccess; + } + + if (stream != nullptr) { + // Stream-ordered async free + cudaError_t result = cudaFreeAsync(ptr, stream); + // Synchronize to ensure operations complete before memory is returned to pool. + // Without this, cudaMallocAsync can return the same memory while operations + // using it are still in flight, causing race conditions and non-determinism. + if (result == cudaSuccess) { + cudaStreamSynchronize(stream); + } + return result; + } + + // Fallback to synchronous free + return cudaFree(ptr); +} + +// ---------------------------------------------------------------------------- +// GPU Pool Logging Helpers (Private) +// ---------------------------------------------------------------------------- + +namespace { + /** + * @brief Log GPU allocation event if logging is enabled + * + * @param result Result code from allocation operation + * @param devPtr Pointer to allocated device memory (or NULL on failure) + * @param size Size of allocation in bytes + * @param timestamp Wall-clock timestamp of operation + * @param file Source file where allocation occurred + * @param line Line number where allocation occurred + * @param tag Pointer identifier tag (e.g. "devBuf") + * @param functionTag Function context tag + */ + inline void logGpuAlloc(cudaError_t result, void** devPtr, size_t size, + double timestamp, const char* file, int line, + const char* tag, const char* functionTag) { + DataManager* dm = (DataManager*)CkLocalNodeBranch(dataManagerID); + if (dm && dm->memLog && dm->bGpuMemLogger) { + MemLogOpType opType = (result == cudaSuccess) ? MEMLOG_ALLOC : MEMLOG_ALLOC_FAIL; + uintptr_t address = (result == cudaSuccess && *devPtr != NULL) ? (uintptr_t)(*devPtr) : 0; + MemLogEvent event(CkMyNode(), opType, size, address, timestamp, file, line, tag, functionTag); + + CmiLock(dm->lockMemLog); + dm->memLog->meTab.push_back(event); + CmiUnlock(dm->lockMemLog); + } + } + + /** + * @brief Log GPU free event if logging is enabled + * + * @param result Result code from free operation + * @param devPtr Device pointer being freed + * @param timestamp Wall-clock timestamp of operation + * @param file Source file where free occurred + * @param line Line number where free occurred + * @param tag Pointer identifier tag + * @param functionTag Function context tag + */ + inline void logGpuFree(cudaError_t result, void* devPtr, double timestamp, + const char* file, int line, const char* tag, + const char* functionTag) { + DataManager* dm = (DataManager*)CkLocalNodeBranch(dataManagerID); + if (dm && dm->memLog && dm->bGpuMemLogger) { + MemLogOpType opType; + if (devPtr == NULL) { + opType = MEMLOG_FREE_SKIP; + } else { + opType = (result == cudaSuccess) ? MEMLOG_FREE : MEMLOG_FREE_FAIL; + } + uintptr_t address = (uintptr_t)devPtr; + MemLogEvent event(CkMyNode(), opType, 0, address, timestamp, file, line, tag, functionTag); + + CmiLock(dm->lockMemLog); + dm->memLog->meTab.push_back(event); + CmiUnlock(dm->lockMemLog); + } + } +} + +// ---------------------------------------------------------------------------- +// GPU Pool API with Logging - Pool-based Allocation +// ---------------------------------------------------------------------------- + +/** + * @brief GPU memory allocation from pool with logging + * + * Allocates from the pool and logs the event if logging is enabled. + * Use the gpuPoolMalloc macro instead of calling this directly. + * + * @param devPtr Output pointer to allocated device memory + * @param size Size to allocate in bytes + * @param stream CUDA stream for ordering + * @param tag Pointer identifier tag + * @param functionTag Function context tag + * @param file Source file + * @param line Line number + * @return cudaError_t from allocation operation + */ +cudaError_t gpuPoolMallocImpl(void** devPtr, size_t size, cudaStream_t stream, + const char* tag, const char* functionTag, + const char* file, int line) { + cudaError_t result = gpuPoolMallocRaw(devPtr, size, stream); + logGpuAlloc(result, devPtr, size, CkWallTimer(), file, line, tag, functionTag); + return result; +} + +/** + * @brief GPU memory deallocation to pool with logging + * + * Frees to the pool and logs the event if logging is enabled. + * Use the gpuPoolFree macro instead of calling this directly. + * + * @param devPtr Device pointer to free + * @param stream CUDA stream for ordering + * @param tag Pointer identifier tag + * @param functionTag Function context tag + * @param file Source file + * @param line Line number + * @return cudaError_t from free operation + */ +cudaError_t gpuPoolFreeImpl(void* devPtr, cudaStream_t stream, + const char* tag, const char* functionTag, + const char* file, int line) { + if (devPtr == NULL) { + logGpuFree(cudaSuccess, devPtr, CkWallTimer(), file, line, tag, functionTag); + return cudaSuccess; + } + + cudaError_t result = gpuPoolFreeRaw(devPtr, stream); + logGpuFree(result, devPtr, CkWallTimer(), file, line, tag, functionTag); + return result; +} + +// ---------------------------------------------------------------------------- +// GPU API with Logging - Direct Allocation (Bypass Pool) +// ---------------------------------------------------------------------------- + +/** + * @brief Direct GPU memory allocation with logging (bypasses pool) + * + * Allocates directly using cudaMalloc and logs the event if logging is enabled. + * Use the gpuMalloc macro instead of calling this directly. + * + * @param devPtr Output pointer to allocated device memory + * @param size Size to allocate in bytes + * @param stream CUDA stream (currently unused, reserved for future use) + * @param tag Pointer identifier tag + * @param functionTag Function context tag + * @param file Source file + * @param line Line number + * @return cudaError_t from allocation operation + */ +cudaError_t gpuMallocImpl(void** devPtr, size_t size, cudaStream_t stream, + const char* tag, const char* functionTag, + const char* file, int line) { + if (devPtr == nullptr) { + return cudaErrorInvalidValue; + } + + if (size == 0) { + *devPtr = nullptr; + return cudaSuccess; + } + + // Use cudaMalloc directly (ignores stream for now) + cudaError_t result = cudaMalloc(devPtr, size); + logGpuAlloc(result, devPtr, size, CkWallTimer(), file, line, tag, functionTag); + return result; +} + +/** + * @brief Direct GPU memory deallocation with logging (bypasses pool) + * + * Frees directly using cudaFree and logs the event if logging is enabled. + * Use the gpuFree macro instead of calling this directly. + * + * @param devPtr Device pointer to free + * @param stream CUDA stream (currently unused, reserved for future use) + * @param tag Pointer identifier tag + * @param functionTag Function context tag + * @param file Source file + * @param line Line number + * @return cudaError_t from free operation + */ +cudaError_t gpuFreeImpl(void* devPtr, cudaStream_t stream, + const char* tag, const char* functionTag, + const char* file, int line) { + if (devPtr == NULL) { + logGpuFree(cudaSuccess, devPtr, CkWallTimer(), file, line, tag, functionTag); + return cudaSuccess; + } + + // Direct free bypassing pool + cudaError_t result = cudaFree(devPtr); + logGpuFree(result, devPtr, CkWallTimer(), file, line, tag, functionTag); + return result; +} + +// ============================================================================ +// Host Memory Pool (Pinned Memory) +// ============================================================================ + +// ---------------------------------------------------------------------------- +// Host Pool Internals (Private State and Helpers) +// ---------------------------------------------------------------------------- + +namespace { + // --- Host Pool Global State --- + + static CmiNodeLock g_hostPoolLock = nullptr; + + // Initialization lock - used to protect lazy initialization of g_hostPoolLock + // Function-local static ensures thread-safe initialization (C++11+) + static CmiNodeLock& getInitLock() { + static CmiNodeLock initLock = CmiCreateLock(); + return initLock; + } + static std::unordered_map> g_hostFreePool; + static std::unordered_map g_hostPointerSizes; + + // Adaptive refill tracking + static std::unordered_map g_bucketMissCount; + static int g_timestepCount = 0; + + // Per-bucket usage tracking for capacity-based trimming + struct BucketUsageInfo { + int lastUsedTimestep = 0; + int recentHits = 0; // Hits in last 10 timesteps + }; + static std::unordered_map g_bucketUsage; + + // --- Host Pool Analytics --- + + struct PoolAnalytics { + // Cumulative data (never reset) + size_t totalAllocated = 0; // Total memory allocated from OS + size_t totalFreed = 0; // Total memory freed back to OS + size_t peakMemory = 0; // Peak memory usage + size_t currentMemory = 0; // Current memory usage + int totalAllocations = 0; // Total allocation calls + int totalFrees = 0; // Total free calls + int totalHits = 0; // Pool hits + int totalMisses = 0; // Pool misses + int totalOsAllocations = 0; // Blocks allocated from OS (misses + refill) + int totalFreedToOS = 0; // Blocks actually freed to OS (not pool) + double totalAllocTime = 0.0; // Total time spent in allocations + double totalFreeTime = 0.0; // Total time spent in frees + + // Per-timestep data (reset each timestep) + int timestepAllocations = 0; + int timestepFrees = 0; + int timestepHits = 0; + int timestepMisses = 0; + int timestepFreedToOS = 0; + + // Per-bucket cumulative tracking + std::unordered_map bucketTotalAllocations; + std::unordered_map bucketTotalMisses; // Cumulative misses per bucket + std::unordered_map bucketTotalOsAllocations; // Cumulative OS allocations per bucket + }; + static PoolAnalytics g_poolAnalytics; + + // Historical tracking for memory leak detection + struct TimestepSnapshot { + int timestep; + size_t totalBuckets; + size_t totalFreeBlocks; + size_t totalFreeBytes; + size_t liveBlocks; + size_t totalAllocated; + size_t totalFreed; + double timestamp; + }; + static std::vector g_timestepHistory; + + // --- Host Pool Helper Functions --- + + /** + * @brief Round up size to appropriate bucket using tiered granularity + * + * Tiered approach to minimize bucket proliferation: + * - Tiny (≤64KB): Fine-grained buckets (256B, 512B, 1KB, 2KB, 4KB, ..., 64KB) + * - Small (64KB-1MB): 64KB granularity + * - Medium (1MB-64MB): 1MB granularity + * - Large (>64MB): 50MB granularity + * + * This drastically reduces unique bucket count for large allocations + * while maintaining efficiency for small allocations. + * + * @param requested Requested allocation size in bytes + * @return Bucket size that will accommodate the requested size + */ + static inline size_t hostBucketSize(size_t requested) { + if (requested == 0) return 0; + + // Tiny allocations: use 256B and 512B buckets + if (requested <= 256) return 256; + if (requested <= 512) return 512; + if (requested <= 1024) return 1024; + if (requested <= 2048) return 2048; + if (requested <= 4 * 1024) return 4 * 1024; + if (requested <= 8 * 1024) return 8 * 1024; + if (requested <= 16 * 1024) return 16 * 1024; + if (requested <= 32 * 1024) return 32 * 1024; + if (requested <= 64 * 1024) return 64 * 1024; + + // Small allocations (64KB - 1MB): 64KB granularity + if (requested <= 1024 * 1024) { + const size_t step = 64 * 1024; + return ((requested + step - 1) / step) * step; + } + + // Medium allocations (1MB - 64MB): 1MB granularity + if (requested <= 64 * 1024 * 1024) { + const size_t step = 1024 * 1024; + return ((requested + step - 1) / step) * step; + } + + // Large allocations (>64MB): 50MB granularity + // This prevents bucket proliferation for very large buffers + const size_t step = 50 * 1024 * 1024; + return ((requested + step - 1) / step) * step; + } + + // --- Host Pool Logging Helpers --- + + /** + * @brief Log host allocation event if logging is enabled + * + * @param result Result code from allocation operation + * @param ptr Pointer to allocated memory (or NULL on failure) + * @param size Size of allocation in bytes + * @param timestamp Wall-clock timestamp of operation + * @param file Source file where allocation occurred + * @param line Line number where allocation occurred + * @param tag Pointer identifier tag + * @param functionTag Function context tag + */ + inline void logHostAlloc(cudaError_t result, void** ptr, size_t size, + double timestamp, const char* file, int line, + const char* tag, const char* functionTag) { + DataManager* dm = (DataManager*)CkLocalNodeBranch(dataManagerID); + if (dm && dm->cpuMemLog && dm->bCpuMemLogger) { + MemLogOpType opType = (result == cudaSuccess) ? MEMLOG_ALLOC : MEMLOG_ALLOC_FAIL; + uintptr_t address = (result == cudaSuccess && *ptr != NULL) ? (uintptr_t)(*ptr) : 0; + MemLogEvent event(CkMyNode(), opType, size, address, timestamp, file, line, tag, functionTag); + + CmiLock(dm->lockCpuMemLog); + dm->cpuMemLog->meTab.push_back(event); + CmiUnlock(dm->lockCpuMemLog); + } + } + + /** + * @brief Log host free event if logging is enabled + * + * @param result Result code from free operation + * @param ptr Pointer being freed + * @param timestamp Wall-clock timestamp of operation + * @param file Source file where free occurred + * @param line Line number where free occurred + * @param tag Pointer identifier tag + * @param functionTag Function context tag + */ + inline void logHostFree(cudaError_t result, void* ptr, double timestamp, + const char* file, int line, const char* tag, + const char* functionTag) { + DataManager* dm = (DataManager*)CkLocalNodeBranch(dataManagerID); + if (dm && dm->cpuMemLog && dm->bCpuMemLogger) { + MemLogOpType opType; + if (ptr == NULL) { + opType = MEMLOG_FREE_SKIP; + } else { + opType = (result == cudaSuccess) ? MEMLOG_FREE : MEMLOG_FREE_FAIL; + } + uintptr_t address = (uintptr_t)ptr; + MemLogEvent event(CkMyNode(), opType, 0, address, timestamp, file, line, tag, functionTag); + + CmiLock(dm->lockCpuMemLog); + dm->cpuMemLog->meTab.push_back(event); + CmiUnlock(dm->lockCpuMemLog); + } + } +} + +// ---------------------------------------------------------------------------- +// Host Pool Initialization & Warmup +// ---------------------------------------------------------------------------- + +/** + * @brief Initialize the host memory pool + * + * Creates the pool lock and pre-warms the pool with common allocation sizes + * ranging from 256B to 64MB across 19 bucket sizes. Total warmup allocates + * and immediately frees blocks to populate the pool's free-lists. + */ +// Shared warmup configuration - single source of truth +// This configuration is used by: +// 1. hostPoolInit() - for actual warmup allocation +// 2. hostPoolAnalyzeWarmupEffectiveness() - for analysis targets +// 3. hostPoolAdaptiveRefill() - for refill target buckets +// Modify only this array to change warmup behavior everywhere +namespace { + struct WarmupConfig { + size_t size; + int count; + }; + + // Warmup configuration - single place to modify bucket sizes and counts + const WarmupConfig g_warmupConfig[] = { + {256, 1000}, // 256 B + {512, 1500}, // 512 B + {1024, 1000}, // 1 KB + {2048, 1000}, // 2 KB + {4096, 1500}, // 4 KB + {8192, 1500}, // 8 KB + {16384, 500}, // 16 KB + {32768, 500}, // 32 KB + {64 * 1024, 500}, // 64 KB + {128 * 1024, 800}, // 128 KB + {256 * 1024, 300}, // 256 KB + {512 * 1024, 300}, // 512 KB + {1 * 1024 * 1024, 1}, // 1 MB + {2 * 1024 * 1024, 1}, // 2 MB + {4 * 1024 * 1024, 1}, // 4 MB + {8 * 1024 * 1024, 1}, // 8 MB + {16 * 1024 * 1024, 5}, // 16 MB + {32 * 1024 * 1024, 1}, // 32 MB + {64 * 1024 * 1024, 5} // 64 MB + }; + + const int g_numWarmupSizes = sizeof(g_warmupConfig) / sizeof(g_warmupConfig[0]); +} + +void hostPoolInit() { + // Initialize lock before concurrent access to avoid race conditions + // in lazy initialization + if (g_hostPoolLock == nullptr) { + g_hostPoolLock = CmiCreateLock(); + } + + size_t totalWarmupBytes = 0; + + for (int i = 0; i < g_numWarmupSizes; i++) { + for (int j = 0; j < g_warmupConfig[i].count; j++) { + void* ptr = nullptr; + cudaError_t res = hostPoolMallocRaw(&ptr, g_warmupConfig[i].size); + if (res != cudaSuccess) { + CkPrintf("WARNING: Host pool warmup failed for size %zu: %s\n", + g_warmupConfig[i].size, cudaGetErrorString(res)); + return; + } + hostPoolFreeRaw(ptr); + totalWarmupBytes += g_warmupConfig[i].size; + } + } + + // Only print from one nodegroup instance to avoid spam in production + if (CkMyNode() == 0) { + CkPrintf("Host memory pool warmed up\n"); + } +} + +// ---------------------------------------------------------------------------- +// Host Pool Raw Operations (Internal - use macros for public API) +// ---------------------------------------------------------------------------- + +/** + * @brief Allocate host memory from pool (raw, internal) + * + * Attempts to reuse a block from the pool. If no block is available, + * allocates a new pinned memory block from the OS at the bucket size. + * + * @param ptr Output pointer to allocated memory + * @param size Size to allocate in bytes + * @return cudaError_t from allocation operation + */ +cudaError_t hostPoolMallocRaw(void** ptr, size_t size) { + if (ptr == nullptr) return cudaErrorInvalidValue; + if (size == 0) { *ptr = nullptr; return cudaSuccess; } + + // Lazy initialization of lock (fallback if hostPoolInit() not called) + // Use double-checked locking pattern with proper synchronization + if (g_hostPoolLock == nullptr) { + CmiNodeLock& initLock = getInitLock(); + CmiLock(initLock); + + // Double-check after acquiring lock + if (g_hostPoolLock == nullptr) { + g_hostPoolLock = CmiCreateLock(); + } + CmiUnlock(initLock); + } + + size_t bucket = hostBucketSize(size); + + // Try to reuse from pool + CmiLock(g_hostPoolLock); + auto it = g_hostFreePool.find(bucket); + if (it != g_hostFreePool.end() && !it->second.empty()) { + void* p = it->second.back(); + it->second.pop_back(); + g_hostPointerSizes[p] = bucket; + + // Update analytics + g_poolAnalytics.totalHits++; + g_poolAnalytics.timestepHits++; + g_poolAnalytics.totalAllocations++; + g_poolAnalytics.timestepAllocations++; + g_poolAnalytics.bucketTotalAllocations[bucket]++; + g_poolAnalytics.currentMemory += bucket; + if (g_poolAnalytics.currentMemory > g_poolAnalytics.peakMemory) { + g_poolAnalytics.peakMemory = g_poolAnalytics.currentMemory; + } + + // Update usage tracking for capacity-based trimming + g_bucketUsage[bucket].lastUsedTimestep = g_timestepCount; + g_bucketUsage[bucket].recentHits++; + + CmiUnlock(g_hostPoolLock); + *ptr = p; + return cudaSuccess; + } + CmiUnlock(g_hostPoolLock); + + // Allocate new pinned block at bucket size (do this unlocked - it's slow) + double startTime = CkWallTimer(); + void* newPtr = nullptr; + cudaError_t res = cudaMallocHost(&newPtr, bucket); + double allocTime = CkWallTimer() - startTime; + + if (res != cudaSuccess) return res; + + // Re-acquire lock for all map and analytics updates + // This ensures thread-safety while keeping the slow cudaMallocHost unlocked + CmiLock(g_hostPoolLock); + + // Track miss for adaptive refill + g_bucketMissCount[bucket]++; + + // Update analytics + g_poolAnalytics.totalMisses++; + g_poolAnalytics.timestepMisses++; + g_poolAnalytics.totalAllocations++; + g_poolAnalytics.timestepAllocations++; + g_poolAnalytics.bucketTotalAllocations[bucket]++; + g_poolAnalytics.bucketTotalMisses[bucket]++; // Track cumulative misses per bucket + g_poolAnalytics.totalOsAllocations++; + g_poolAnalytics.bucketTotalOsAllocations[bucket]++; + g_poolAnalytics.totalAllocTime += allocTime; + g_poolAnalytics.totalAllocated += bucket; + g_poolAnalytics.currentMemory += bucket; + if (g_poolAnalytics.currentMemory > g_poolAnalytics.peakMemory) { + g_poolAnalytics.peakMemory = g_poolAnalytics.currentMemory; + } + + // Register pointer in pool tracking + g_hostPointerSizes[newPtr] = bucket; + + CmiUnlock(g_hostPoolLock); + *ptr = newPtr; + return cudaSuccess; +} + +/** + * @brief Free host memory back to pool (raw, internal) + * + * Returns the memory block to the pool for reuse. If the pointer + * is unknown or the pool is not initialized, falls back to cudaFreeHost. + * + * @param ptr Pointer to free + * @return cudaError_t from operation + */ +cudaError_t hostPoolFreeRaw(void* ptr) { + if (ptr == nullptr) return cudaSuccess; + + if (g_hostPoolLock == nullptr) { + // Pool not initialized; free directly + return cudaFreeHost(ptr); + } + + CmiLock(g_hostPoolLock); + auto it = g_hostPointerSizes.find(ptr); + if (it == g_hostPointerSizes.end()) { + CmiUnlock(g_hostPoolLock); + // Unknown pointer - free directly to be safe + return cudaFreeHost(ptr); + } + + size_t bucket = it->second; + g_hostPointerSizes.erase(it); + g_hostFreePool[bucket].push_back(ptr); + + // Update analytics + g_poolAnalytics.totalFrees++; + g_poolAnalytics.timestepFrees++; + g_poolAnalytics.currentMemory -= bucket; + + CmiUnlock(g_hostPoolLock); + return cudaSuccess; +} + +// ---------------------------------------------------------------------------- +// Host Pool API with Logging - Pool-based Allocation +// ---------------------------------------------------------------------------- + +/** + * @brief Host memory allocation from pool with logging + * + * Allocates from the host pool and logs the event if logging is enabled. + * Use the hostPoolMalloc macro instead of calling this directly. + * + * @param ptr Output pointer to allocated memory + * @param size Size to allocate in bytes + * @param tag Pointer identifier tag + * @param functionTag Function context tag + * @param file Source file + * @param line Line number + * @return cudaError_t from allocation operation + */ +cudaError_t hostPoolMallocImpl(void** ptr, size_t size, + const char* tag, const char* functionTag, + const char* file, int line) { + cudaError_t result = hostPoolMallocRaw(ptr, size); + logHostAlloc(result, ptr, size, CkWallTimer(), file, line, tag, functionTag); + return result; +} + +/** + * @brief Host memory deallocation to pool with logging + * + * Frees to the pool and logs the event if logging is enabled. + * Use the hostPoolFree macro instead of calling this directly. + * + * @param ptr Pointer to free + * @param tag Pointer identifier tag + * @param functionTag Function context tag + * @param file Source file + * @param line Line number + * @return cudaError_t from free operation + */ +cudaError_t hostPoolFreeImpl(void* ptr, + const char* tag, const char* functionTag, + const char* file, int line) { + if (ptr == NULL) { + logHostFree(cudaSuccess, ptr, CkWallTimer(), file, line, tag, functionTag); + return cudaSuccess; + } + + cudaError_t result = hostPoolFreeRaw(ptr); + logHostFree(result, ptr, CkWallTimer(), file, line, tag, functionTag); + return result; +} + +// ---------------------------------------------------------------------------- +// Host API with Logging - Direct Allocation (Bypass Pool) +// ---------------------------------------------------------------------------- + +/** + * @brief Direct host memory allocation with logging (bypasses pool) + * + * Allocates directly using cudaMallocHost and logs the event if logging is enabled. + * Use the hostMalloc macro instead of calling this directly. + * + * @param ptr Output pointer to allocated memory + * @param size Size to allocate in bytes + * @param tag Pointer identifier tag + * @param functionTag Function context tag + * @param file Source file + * @param line Line number + * @return cudaError_t from allocation operation + */ +cudaError_t hostMallocImpl(void** ptr, size_t size, + const char* tag, const char* functionTag, + const char* file, int line) { + if (ptr == nullptr) { + return cudaErrorInvalidValue; + } + + if (size == 0) { + *ptr = nullptr; + return cudaSuccess; + } + + // Use cudaMallocHost directly + cudaError_t result = cudaMallocHost(ptr, size); + logHostAlloc(result, ptr, size, CkWallTimer(), file, line, tag, functionTag); + return result; +} + +/** + * @brief Direct host memory deallocation with logging (bypasses pool) + * + * Frees directly using cudaFreeHost and logs the event if logging is enabled. + * Use the hostFree macro instead of calling this directly. + * + * @param ptr Pointer to free + * @param tag Pointer identifier tag + * @param functionTag Function context tag + * @param file Source file + * @param line Line number + * @return cudaError_t from free operation + */ +cudaError_t hostFreeImpl(void* ptr, + const char* tag, const char* functionTag, + const char* file, int line) { + if (ptr == NULL) { + logHostFree(cudaSuccess, ptr, CkWallTimer(), file, line, tag, functionTag); + return cudaSuccess; + } + + // Direct free bypassing pool + cudaError_t result = cudaFreeHost(ptr); + logHostFree(result, ptr, CkWallTimer(), file, line, tag, functionTag); + return result; +} + +// ---------------------------------------------------------------------------- +// Host Pool Maintenance & Trimming +// ---------------------------------------------------------------------------- + +/** + * @brief Trim host pool to target capacity + * + * Frees excess blocks from the pool to the OS, prioritizing large buckets + * and respecting usage-based protection for recently-used buckets. + * + * @param targetCapacityGB Target maximum pool capacity in GB + * @param minCapacityPerBucketMB Minimum capacity to preserve per active bucket in MB + * @return Total bytes freed to the OS + */ +size_t hostPoolTrim(double targetCapacityGB, size_t minCapacityPerBucketMB) { + if (g_hostPoolLock == nullptr) { + return 0; + } + + size_t targetCapacity = (size_t)(targetCapacityGB * 1024.0 * 1024.0 * 1024.0); + size_t minCapacityPerBucket = minCapacityPerBucketMB * 1024 * 1024; + + CmiLock(g_hostPoolLock); + + // Calculate current free pool capacity + size_t currentCapacity = 0; + for (const auto &kv : g_hostFreePool) { + currentCapacity += kv.first * kv.second.size(); + } + + // If under target, nothing to trim + if (currentCapacity <= targetCapacity) { + CmiUnlock(g_hostPoolLock); + return 0; + } + + // Build list of (bucket, trimmableBlocks, bucketSize) sorted by bucketSize descending + struct TrimCandidate { + size_t bucket; + int trimmableBlocks; + size_t bytesPerBlock; + }; + std::vector candidates; + + for (auto &kv : g_hostFreePool) { + size_t bucket = kv.first; + int currentBlocks = kv.second.size(); + + // Calculate protection for this bucket + int protectedBlocks = 0; + auto usageIt = g_bucketUsage.find(bucket); + int recentHits = 0; + if (usageIt != g_bucketUsage.end() && usageIt->second.recentHits > 0) { + recentHits = usageIt->second.recentHits; + // Protect based on recent usage, but limited by capacity + int capacityLimit = (int)(minCapacityPerBucket / bucket); + protectedBlocks = std::min(usageIt->second.recentHits, capacityLimit); + } + + int trimmable = std::max(0, currentBlocks - protectedBlocks); + + if (trimmable > 0) { + candidates.push_back({bucket, trimmable, bucket}); + } + } + + // Sort by bucket size descending (trim largest first) + std::sort(candidates.begin(), candidates.end(), + [](const TrimCandidate& a, const TrimCandidate& b) { + return a.bytesPerBlock > b.bytesPerBlock; + }); + + // Greedily trim blocks from largest buckets until target reached + // Collect pointers to free while holding lock, then free them unlocked + struct BlockToFree { + void* ptr; + size_t bucket; + }; + std::vector blocksToFree; + size_t totalFreed = 0; + + for (auto &candidate : candidates) { + if (currentCapacity <= targetCapacity) { + break; + } + + auto &freeBlocks = g_hostFreePool[candidate.bucket]; + int blocksToFreeCount = std::min(candidate.trimmableBlocks, + (int)((currentCapacity - targetCapacity) / candidate.bucket) + 1); + + for (int i = 0; i < blocksToFreeCount && !freeBlocks.empty(); i++) { + void* ptr = freeBlocks.back(); + freeBlocks.pop_back(); + g_hostPointerSizes.erase(ptr); + blocksToFree.push_back({ptr, candidate.bucket}); + + totalFreed += candidate.bucket; + currentCapacity -= candidate.bucket; + } + } + + // Release lock before slow cudaFreeHost operations + CmiUnlock(g_hostPoolLock); + + // Free all blocks unlocked (slow operation) + for (const auto& block : blocksToFree) { + cudaFreeHost(block.ptr); + } + + // Re-acquire lock for analytics updates + CmiLock(g_hostPoolLock); + int blocksTrimmed = blocksToFree.size(); + for (const auto& block : blocksToFree) { + g_poolAnalytics.totalFreedToOS++; + g_poolAnalytics.timestepFreedToOS++; + g_poolAnalytics.totalFreed += block.bucket; + } + CmiUnlock(g_hostPoolLock); + + return totalFreed; +} + +/** + * @brief Adaptive refill of hot buckets between timesteps + * + * Refills buckets that had misses during the previous timestep to maintain + * non-blocking behavior. Only refills if misses occurred and capacity is low. + */ +void hostPoolAdaptiveRefill() { + if (g_hostPoolLock == nullptr) return; + + CmiLock(g_hostPoolLock); + + // Use shared warmup configuration for refill targets + // Collect refill requirements while holding lock + struct RefillRequest { + size_t bucket; + int count; + }; + std::vector refillRequests; + + for (int i = 0; i < g_numWarmupSizes; i++) { + size_t bucket = g_warmupConfig[i].size; + auto missIt = g_bucketMissCount.find(bucket); + + // Only refill if we had misses and current free count is low + if (missIt != g_bucketMissCount.end() && missIt->second > 0) { + auto freeIt = g_hostFreePool.find(bucket); + int currentFree = (freeIt != g_hostFreePool.end()) ? freeIt->second.size() : 0; + + // Refill up to 3 blocks if we had misses and have < 2 free + int refillCount = std::min(3 - currentFree, missIt->second); + if (refillCount > 0) { + refillRequests.push_back({bucket, refillCount}); + } + } + } + + // Reset miss counts for next timestep + g_bucketMissCount.clear(); + g_timestepCount++; + + // Decay recent hits to allow trimming of less-used buckets + // Active buckets maintain protection through ongoing hits + // Inactive buckets gradually lose protection and become trimmable + for (auto &kv : g_bucketUsage) { + if (kv.second.recentHits > 0) { + kv.second.recentHits--; + } + } + + CmiUnlock(g_hostPoolLock); + + // Allocate blocks unlocked (slow operation) + struct AllocatedBlock { + void* ptr; + size_t bucket; + }; + std::vector allocatedBlocks; + int totalRefilled = 0; + size_t totalRefillBytes = 0; + + for (const auto& req : refillRequests) { + for (int j = 0; j < req.count; j++) { + void* ptr = nullptr; + cudaError_t res = cudaMallocHost(&ptr, req.bucket); + if (res == cudaSuccess) { + allocatedBlocks.push_back({ptr, req.bucket}); + totalRefilled++; + totalRefillBytes += req.bucket; + } else { + break; // Stop if allocation fails + } + } + } + + // Re-acquire lock to add allocated blocks to pool + if (!allocatedBlocks.empty()) { + CmiLock(g_hostPoolLock); + for (const auto& block : allocatedBlocks) { + g_hostFreePool[block.bucket].push_back(block.ptr); + // Note: Free pool blocks are not tracked in g_hostPointerSizes + // Only live (allocated) blocks are tracked there + g_poolAnalytics.totalOsAllocations++; + g_poolAnalytics.bucketTotalOsAllocations[block.bucket]++; + g_poolAnalytics.totalAllocated += block.bucket; + } + CmiUnlock(g_hostPoolLock); + } + + if (totalRefilled > 0) { + DataManager* dm = (DataManager*)CkLocalNodeBranch(dataManagerID); + if (dm && dm->bHostPoolDebug) { + CkPrintf("[HostPool] Refilled %d blocks (%.1f MB) for hot buckets\n", + totalRefilled, totalRefillBytes / (1024.0 * 1024.0)); + } + } + + // Reset per-timestep analytics after reporting/refill completes + CmiLock(g_hostPoolLock); + g_poolAnalytics.timestepAllocations = 0; + g_poolAnalytics.timestepFrees = 0; + g_poolAnalytics.timestepHits = 0; + g_poolAnalytics.timestepMisses = 0; + g_poolAnalytics.timestepFreedToOS = 0; + CmiUnlock(g_hostPoolLock); +} + +// ---------------------------------------------------------------------------- +// Host Pool Analytics & Reporting +// ---------------------------------------------------------------------------- + +/** + * @brief Report host pool statistics with enhanced analytics + * + * Prints current pool state including bucket count, free blocks, total memory, + * hit rates, per-timestep statistics, and per-bucket distribution. + * Also performs pool accounting validation and stores timestep snapshot. + * + * @param prefix String prefix for output lines (for formatting) + * @param targetCapacityGB Target pool capacity in GB (for trim reporting) + */ +void hostPoolReportStats(const char* prefix, double targetCapacityGB) { + if (g_hostPoolLock == nullptr) { + CkPrintf("%s Host pool not initialized\n", prefix); + return; + } + + CmiLock(g_hostPoolLock); + + size_t totalBlocks = 0; + size_t totalBytes = 0; + size_t numBuckets = g_hostFreePool.size(); + + for (const auto &kv : g_hostFreePool) { + size_t bucketSize = kv.first; + size_t numBlocks = kv.second.size(); + totalBlocks += numBlocks; + totalBytes += bucketSize * numBlocks; + } + + size_t liveBlocks = g_hostPointerSizes.size(); + + // Copy all analytics values while holding lock to avoid race conditions + int totalAllocations = g_poolAnalytics.totalAllocations; + int totalHits = g_poolAnalytics.totalHits; + int totalMisses = g_poolAnalytics.totalMisses; + int timestepAllocations = g_poolAnalytics.timestepAllocations; + int timestepFrees = g_poolAnalytics.timestepFrees; + int timestepHits = g_poolAnalytics.timestepHits; + int timestepMisses = g_poolAnalytics.timestepMisses; + int timestepFreedToOS = g_poolAnalytics.timestepFreedToOS; + int totalOsAllocations = g_poolAnalytics.totalOsAllocations; + size_t peakMemory = g_poolAnalytics.peakMemory; + size_t currentMemory = g_poolAnalytics.currentMemory; + double totalAllocTime = g_poolAnalytics.totalAllocTime; + size_t totalAllocated = g_poolAnalytics.totalAllocated; + size_t totalFreed = g_poolAnalytics.totalFreed; + size_t totalFreedToOS = g_poolAnalytics.totalFreedToOS; + int timestepCount = g_timestepCount; + + CmiUnlock(g_hostPoolLock); + + // Calculate capacity usage + double totalBytesGB = totalBytes / (1024.0 * 1024.0 * 1024.0); + double capacityUsagePct = (targetCapacityGB > 0) ? (totalBytesGB / targetCapacityGB * 100.0) : 0.0; + const char* capacityStatus = (totalBytesGB > targetCapacityGB) ? "OVER" : "OK"; + + // Pool summary header + CkPrintf("\n"); + CkPrintf("%s ═══════════════════════════════════════════════════════════════\n", prefix); + CkPrintf("%s HOST POOL STATUS\n", prefix); + CkPrintf("%s ═══════════════════════════════════════════════════════════════\n", prefix); + + // Capacity metrics + CkPrintf("%s Free Pool Capacity: %.2f GB / %.2f GB target (%.1f%% %s)\n", + prefix, totalBytesGB, targetCapacityGB, capacityUsagePct, capacityStatus); + CkPrintf("%s Pool Structure: %zu buckets, %zu free blocks, %zu live blocks\n", + prefix, numBuckets, totalBlocks, liveBlocks); + + // Hit rate and efficiency + double hitRate = (totalAllocations > 0) ? + (100.0 * totalHits) / totalAllocations : 0.0; + + CkPrintf("%s Hit Rate: %.1f%% (%d hits, %d misses)\n", + prefix, hitRate, totalHits, totalMisses); + + // Memory usage + CkPrintf("%s Memory Usage: Peak %.2f GB, Current %.2f GB\n", + prefix, + peakMemory / (1024.0 * 1024.0 * 1024.0), + currentMemory / (1024.0 * 1024.0 * 1024.0)); + + // This timestep activity + CkPrintf("%s This Step: %d allocs, %d frees, %d hits, %d misses, %d freed to OS\n", + prefix, timestepAllocations, timestepFrees, + timestepHits, timestepMisses, timestepFreedToOS); + + // Cumulative stats + CkPrintf("%s Cumulative: %d allocs, %.3f ms avg time, %.2f GB from OS\n", + prefix, totalAllocations, + (totalAllocTime * 1000.0) / std::max(1, totalAllocations), + totalAllocated / (1024.0 * 1024.0 * 1024.0)); + + // Pool accounting check + // Total blocks owned by pool = (allocated from OS) - (freed to OS) + // These blocks are either: (1) free in pool, or (2) currently checked out (live) + size_t totalOwnedBlocks = totalOsAllocations - totalFreedToOS; + size_t accountedBlocks = totalBlocks + liveBlocks; + + if (totalOwnedBlocks != accountedBlocks) { + size_t diff = (totalOwnedBlocks > accountedBlocks) ? + (totalOwnedBlocks - accountedBlocks) : (accountedBlocks - totalOwnedBlocks); + CkPrintf("%s WARNING: Pool accounting mismatch! Owned: %zu, Accounted: %zu (Free: %zu + Live: %zu), Diff: %zu\n", + prefix, totalOwnedBlocks, accountedBlocks, totalBlocks, liveBlocks, diff); + } + + // Per-bucket diagnostics (sorted by bucket size) + CkPrintf("%s Per-Bucket Distribution:\n", prefix); + + CmiLock(g_hostPoolLock); + + // Collect and sort buckets by size + std::vector bucketSizes; + for (const auto &kv : g_hostFreePool) { + bucketSizes.push_back(kv.first); + } + + // Track live blocks per bucket for accurate usage + std::unordered_map liveCounts; + for (const auto &kv : g_hostPointerSizes) { + liveCounts[kv.second]++; + if (g_hostFreePool.find(kv.second) == g_hostFreePool.end()) { + bucketSizes.push_back(kv.second); + } + } + std::sort(bucketSizes.begin(), bucketSizes.end()); + bucketSizes.erase(std::unique(bucketSizes.begin(), bucketSizes.end()), bucketSizes.end()); + + // Print header + CkPrintf("%s %10s | %6s | %6s | %6s | %10s | %8s\n", prefix, + "Bucket", "Free", "Live", "Misses", "TotalAllocs", "Usage"); + CkPrintf("%s %s\n", prefix, "----------------------------------------------------------------"); + + for (size_t bucket : bucketSizes) { + int freeCount = g_hostFreePool[bucket].size(); + int totalMisses = g_poolAnalytics.bucketTotalMisses.count(bucket) ? + g_poolAnalytics.bucketTotalMisses[bucket] : 0; + int totalAllocs = g_poolAnalytics.bucketTotalAllocations.count(bucket) ? + g_poolAnalytics.bucketTotalAllocations[bucket] : 0; + int liveCount = liveCounts.count(bucket) ? liveCounts[bucket] : 0; + + // Format bucket size nicely + auto bucketStr = (bucket < 1024) ? make_formatted_string("%zu B", bucket) + : (bucket < 1024 * 1024) ? make_formatted_string("%zu KB", bucket / 1024) + : make_formatted_string("%zu MB", bucket / (1024 * 1024)); + + // Calculate usage percentage based on current pool state. + // This avoids negative values when free blocks exceed historical misses + // (e.g., due to adaptive refills). + int totalCurrent = freeCount + liveCount; + double usagePct = (totalCurrent > 0) ? (100.0 * liveCount / totalCurrent) : 0.0; + + CkPrintf("%s %10s | %6d | %6d | %6d | %10d | %7.1f%%\n", prefix, + bucketStr.c_str(), freeCount, liveCount, totalMisses, totalAllocs, usagePct); + } + + CkPrintf("%s %s\n", prefix, "----------------------------------------------------------------"); + CkPrintf("%s ═══════════════════════════════════════════════════════════════\n", prefix); + CkPrintf("\n"); + + CmiUnlock(g_hostPoolLock); // Release lock after per-bucket diagnostics + + // Store timestep snapshot for historical tracking (protected by lock) + TimestepSnapshot snapshot; + snapshot.timestep = timestepCount; + snapshot.totalBuckets = numBuckets; + snapshot.totalFreeBlocks = totalBlocks; + snapshot.totalFreeBytes = totalBytes; + snapshot.liveBlocks = liveBlocks; + snapshot.totalAllocated = totalAllocated; + snapshot.totalFreed = totalFreed; + snapshot.timestamp = CkWallTimer(); + + // Re-acquire lock to update history + CmiLock(g_hostPoolLock); + g_timestepHistory.push_back(snapshot); + + // Keep only last 100 timesteps to prevent memory growth + if (g_timestepHistory.size() > 100) { + g_timestepHistory.erase(g_timestepHistory.begin()); + } + CmiUnlock(g_hostPoolLock); +} + +/** + * @brief Analyze warmup effectiveness by comparing targets vs. actual usage + * + * Identifies over-provisioned and under-provisioned buckets based on + * warmup targets vs. actual free counts and misses. Provides recommendations + * for adjusting warmup parameters. + */ +void hostPoolAnalyzeWarmupEffectiveness() { + if (g_hostPoolLock == nullptr) { + CkPrintf("[PoolAnalysis] Host pool not initialized\n"); + return; + } + + // Build warmup targets from shared configuration - always in sync + std::unordered_map warmupTargets; + for (int i = 0; i < g_numWarmupSizes; i++) { + warmupTargets[g_warmupConfig[i].size] = g_warmupConfig[i].count; + } + + CkPrintf("\n=== WARMUP EFFECTIVENESS ANALYSIS ===\n"); + + CmiLock(g_hostPoolLock); + + // Collect buckets sorted by size + std::vector bucketSizes; + for (const auto &kv : g_hostFreePool) { + bucketSizes.push_back(kv.first); + } + std::sort(bucketSizes.begin(), bucketSizes.end()); + + CkPrintf("Bucket Analysis (Target vs. Actual):\n"); + CkPrintf(" %10s | %8s | %6s | %6s | %10s | Status\n", + "Bucket", "Target", "Free", "Misses", "TotalAllocs", ""); + CkPrintf(" --------------------------------------------------------------------\n"); + + int overProvisionedCount = 0; + int underProvisionedCount = 0; + int wellProvisionedCount = 0; + + for (size_t bucket : bucketSizes) { + int target = warmupTargets.count(bucket) ? warmupTargets[bucket] : 0; + int freeCount = g_hostFreePool[bucket].size(); + int totalMisses = g_poolAnalytics.bucketTotalMisses.count(bucket) ? + g_poolAnalytics.bucketTotalMisses[bucket] : 0; + int totalAllocs = g_poolAnalytics.bucketTotalAllocations.count(bucket) ? + g_poolAnalytics.bucketTotalAllocations[bucket] : 0; + + // Format bucket size + char bucketStr[16]; + if (bucket < 1024) { + snprintf(bucketStr, sizeof(bucketStr), "%zu B", bucket); + } else if (bucket < 1024 * 1024) { + snprintf(bucketStr, sizeof(bucketStr), "%zu KB", bucket / 1024); + } else { + snprintf(bucketStr, sizeof(bucketStr), "%zu MB", bucket / (1024 * 1024)); + } + + // Determine status + const char* status = ""; + if (target == 0) { + // Bucket not in warmup - likely dynamic allocation + status = "DYNAMIC"; + } else if (freeCount > target * 1.5) { + // Over-provisioned: way more free than target + status = "OVER"; + overProvisionedCount++; + } else if (totalMisses > target * 0.1) { + // Under-provisioned: significant misses after warmup + status = "UNDER"; + underProvisionedCount++; + } else { + // Well-provisioned: close to target, few misses + status = "OK"; + wellProvisionedCount++; + } + + CkPrintf(" %10s | %8d | %6d | %6d | %10d | %s\n", + bucketStr, target, freeCount, totalMisses, totalAllocs, status); + } + + CmiUnlock(g_hostPoolLock); + + // Summary + CkPrintf("\nSummary:\n"); + CkPrintf(" Well-provisioned buckets: %d\n", wellProvisionedCount); + CkPrintf(" Over-provisioned buckets: %d (reduce warmup counts)\n", overProvisionedCount); + CkPrintf(" Under-provisioned buckets: %d (increase warmup counts)\n", underProvisionedCount); + + // Recommendations + if (overProvisionedCount > 0) { + CkPrintf("\nRecommendation: Reduce warmup counts for OVER buckets to save initial memory.\n"); + } + if (underProvisionedCount > 0) { + CkPrintf("\nRecommendation: Increase warmup counts for UNDER buckets to reduce misses.\n"); + } + if (wellProvisionedCount == bucketSizes.size()) { + CkPrintf("\nExcellent! All buckets are well-provisioned.\n"); + } + + CkPrintf("=== END WARMUP ANALYSIS ===\n\n"); +} + +/** + * @brief Analyze memory growth trends over time + * + * Reports memory growth patterns and detects potential leaks by analyzing + * historical timestep snapshots. Identifies concerning trends in bucket + * proliferation, memory growth, and live block accumulation. + */ +void hostPoolAnalyzeGrowth() { + if (g_hostPoolLock == nullptr) { + CkPrintf("[PoolAnalysis] Host pool not initialized\n"); + return; + } + + // Protect all access to g_timestepHistory with lock + CmiLock(g_hostPoolLock); + + if (g_timestepHistory.size() < 2) { + CmiUnlock(g_hostPoolLock); + CkPrintf("[PoolAnalysis] Not enough data for growth analysis (need at least 2 timesteps)\n"); + return; + } + + // Copy data we need while holding lock + size_t historySize = g_timestepHistory.size(); + TimestepSnapshot first = g_timestepHistory.front(); + TimestepSnapshot last = g_timestepHistory.back(); + + // Copy recent snapshots if available + std::vector recentSnapshots; + if (historySize >= 5) { + for (int i = historySize - 5; i < (int)historySize; i++) { + recentSnapshots.push_back(g_timestepHistory[i]); + } + } + + CmiUnlock(g_hostPoolLock); + + // Now do calculations and printing unlocked (using copied data) + CkPrintf("\n=== POOL GROWTH ANALYSIS ===\n"); + + // Calculate growth rates + double timeSpan = last.timestamp - first.timestamp; + size_t bucketGrowth = last.totalBuckets - first.totalBuckets; + size_t memoryGrowth = last.totalFreeBytes - first.totalFreeBytes; + size_t liveGrowth = last.liveBlocks - first.liveBlocks; + + CkPrintf("Time span: %.1f seconds (%d timesteps)\n", timeSpan, last.timestep - first.timestep + 1); + CkPrintf("Bucket growth: %zu buckets (%.1f per timestep)\n", + bucketGrowth, (double)bucketGrowth / (last.timestep - first.timestep + 1)); + CkPrintf("Memory growth: %.2f GB (%.3f GB per timestep)\n", + memoryGrowth / (1024.0 * 1024.0 * 1024.0), + memoryGrowth / (1024.0 * 1024.0 * 1024.0) / (last.timestep - first.timestep + 1)); + CkPrintf("Live block growth: %zu blocks (%.1f per timestep)\n", + liveGrowth, (double)liveGrowth / (last.timestep - first.timestep + 1)); + + // Detect concerning patterns + if (bucketGrowth > 10) { + CkPrintf("WARNING: High bucket proliferation detected! Consider adjusting bucket granularity.\n"); + } + + if (memoryGrowth > 1024 * 1024 * 1024) { // > 1GB + CkPrintf("WARNING: Significant memory growth detected! Check for memory leaks.\n"); + } + + if (liveGrowth > 0) { + CkPrintf("WARNING: Live block count increasing - potential memory leak!\n"); + } + + // Show recent trend (last 5 timesteps) + if (!recentSnapshots.empty()) { + CkPrintf("\nRecent trend (last 5 timesteps):\n"); + for (const auto& snap : recentSnapshots) { + CkPrintf(" Step %d: %zu buckets, %.2f GB free, %zu live\n", + snap.timestep, snap.totalBuckets, + snap.totalFreeBytes / (1024.0 * 1024.0 * 1024.0), snap.liveBlocks); + } + } + + CkPrintf("=== END ANALYSIS ===\n\n"); +} + +// ============================================================================ +// MEMORY LOGGING INFRASTRUCTURE +// ============================================================================ + +// ---------------------------------------------------------------------------- +// DataManager Log Initialization +// ---------------------------------------------------------------------------- + +/** + * @brief Initialize GPU memory log on DataManager + * + * Called collectively across all DataManagers to set the log filename + * and enable/disable GPU memory logging. + * + * @param _fileName Output filename for GPU memory log + * @param bGpuMemLoggerFlag Flag to enable (1) or disable (0) GPU logging + * @param cb Callback to contribute when initialization completes + */ +void DataManager::initMemLog(std::string _fileName, int bGpuMemLoggerFlag, const CkCallback &cb) { + CmiLock(lockMemLog); + if (memLog != nullptr) { + memLog->fileName = _fileName; + } else { + CkPrintf("WARNING PE %d: memLog is NULL in initMemLog! Cannot set filename.\n", CkMyPe()); + } + bGpuMemLogger = bGpuMemLoggerFlag; + CmiUnlock(lockMemLog); + contribute(cb); +} + +/** + * @brief Initialize CPU memory log on DataManager + * + * Called collectively across all DataManagers to set the log filename + * and enable/disable CPU memory logging. + * + * @param _fileName Output filename for CPU memory log + * @param bCpuMemLoggerFlag Flag to enable (1) or disable (0) CPU logging + * @param cb Callback to contribute when initialization completes + */ +void DataManager::initCpuMemLog(std::string _fileName, int bCpuMemLoggerFlag, const CkCallback &cb) { + CmiLock(lockCpuMemLog); + if (cpuMemLog != nullptr) { + cpuMemLog->fileName = _fileName; + } else { + CkPrintf("WARNING PE %d: cpuMemLog is NULL in initCpuMemLog! Cannot set filename.\n", CkMyPe()); + } + bCpuMemLogger = bCpuMemLoggerFlag; + CmiUnlock(lockCpuMemLog); + contribute(cb); +} + +// ---------------------------------------------------------------------------- +// Main Chare Log Initialization +// ---------------------------------------------------------------------------- + +/** + * @brief Initialize GPU memory log file on PE 0 + * + * Called from Main chare to create the GPU memory log file and write header. + * Sends configuration to all DataManagers before creating the file. + */ +void Main::initMemLog() { + std::string memLogFile = "gpumemlog.out"; + + // Send filename and flag to all DataManagers and wait + dMProxy.initMemLog(memLogFile, param.bGpuMemLogger, CkCallbackResumeThread()); + + // PE 0 creates the file and writes header + if (CkMyPe() == 0) { + FILE* fpLog = CmiFopen(memLogFile.c_str(), "w"); + fprintf(fpLog, "# ChaNGa Memory Log v1.1\n"); + fprintf(fpLog, "# NodeID OpType Size Address Timestamp File:Line PointerID FunctionTag\n"); + int close_err = CmiFclose(fpLog); + if (close_err != 0) { + CkPrintf("WARNING: PE 0 failed to close gpumemlog file: %s (Error %d)\n", memLogFile.c_str(), close_err); + } + } +} + +/** + * @brief Initialize CPU memory log file on PE 0 + * + * Called from Main chare to create the CPU memory log file and write header. + * Sends configuration to all DataManagers before creating the file. + */ +void Main::initCpuMemLog() { + std::string cpuMemLogFile = "cpumemlog.out"; + + // Send filename and flag to all DataManagers and wait + dMProxy.initCpuMemLog(cpuMemLogFile, param.bCpuMemLogger, CkCallbackResumeThread()); + + // PE 0 creates the file and writes header + if (CkMyPe() == 0) { + FILE* fpLog = CmiFopen(cpuMemLogFile.c_str(), "w"); + fprintf(fpLog, "# ChaNGa CPU Memory Log v1.0\n"); + fprintf(fpLog, "# NodeID OpType Size Address Timestamp File:Line PointerID FunctionTag\n"); + int close_err = CmiFclose(fpLog); + if (close_err != 0) { + CkPrintf("WARNING: PE 0 failed to close CPU memlog file: %s (Error %d)\n", cpuMemLogFile.c_str(), close_err); + } + } +} + +// ---------------------------------------------------------------------------- +// DataManager Log Flushing +// ---------------------------------------------------------------------------- + +/** + * @brief Flush GPU memlog sequentially across nodes + * + * Sequentially flushes GPU memory logs from each node to avoid file corruption. + * Each node triggers the next node in sequence before invoking the callback. + * + * @param cb Callback to invoke after all nodes have flushed + */ +void DataManager::flushMemLog(const CkCallback& cb) { + if (memLog) { + memLog->flush(); + } else { + CkPrintf("WARNING Node %d: memLog is NULL in flushMemLog! Skipping flush.\n", thisIndex); + } + + // Sequential flushing to avoid file corruption + if(thisIndex != CkNumNodes()-1) { + thisProxy[thisIndex + 1].flushMemLog(cb); + } else { + cb.send(); + } +} + +/** + * @brief Flush CPU memlog sequentially across nodes + * + * Sequentially flushes CPU memory logs from each node to avoid file corruption. + * Each node triggers the next node in sequence before invoking the callback. + * + * @param cb Callback to invoke after all nodes have flushed + */ +void DataManager::flushCpuMemLog(const CkCallback& cb) { + if (cpuMemLog) { + cpuMemLog->flush(); + } else { + CkPrintf("WARNING Node %d: cpuMemLog is NULL in flushCpuMemLog! Skipping flush.\n", thisIndex); + } + + // Sequential flushing to avoid file corruption + if(thisIndex != CkNumNodes()-1) { + thisProxy[thisIndex + 1].flushCpuMemLog(cb); + } else { + cb.send(); + } +} + +// ---------------------------------------------------------------------------- +// MemLog Implementation +// ---------------------------------------------------------------------------- + +/** + * @brief Write buffered memory log events to file + * + * Flushes all buffered memory log events to the log file in append mode. + * Converts operation types to human-readable strings and formats each event. + * Clears the event buffer after writing. + */ +void MemLog::flush() { + if (meTab.empty()) { + return; + } + + FILE* outfile = CmiFopen(fileName.c_str(), "a"); + if (outfile == NULL) { + CkPrintf("WARNING: Could not open memlog file '%s' for appending.\n", fileName.c_str()); + return; + } + + // Write all buffered events + for (const auto& event : meTab) { + const char* opTypeStr; + switch (event.opType) { + case MEMLOG_ALLOC: opTypeStr = "ALLOC"; break; + case MEMLOG_FREE: opTypeStr = "FREE "; break; + case MEMLOG_ALLOC_FAIL: opTypeStr = "ALLOC_F"; break; + case MEMLOG_FREE_FAIL: opTypeStr = "FREE_F "; break; + case MEMLOG_FREE_SKIP: opTypeStr = "FREE_S "; break; + default: opTypeStr = "UNKNOWN"; break; + } + + fprintf(outfile, "%d %s %zu %p %.6f %s %s %s\n", + event.nodeId, + opTypeStr, + event.size, + (void*)event.address, + event.timestamp, + event.location.c_str(), + event.pointerId.c_str(), + event.functionTag.c_str()); + } + + int result = CmiFclose(outfile); + if (result != 0) { + CkPrintf("WARNING: Failed to close memlog file '%s' properly (Error %d).\n", fileName.c_str(), result); + } + + meTab.clear(); +} + +#endif // CUDA + diff --git a/MemoryPool.h b/MemoryPool.h new file mode 100644 index 00000000..4062dbe9 --- /dev/null +++ b/MemoryPool.h @@ -0,0 +1,374 @@ +#ifndef MEMORY_POOL_H +#define MEMORY_POOL_H + +#ifdef CUDA +#include + +// ============================================================================ +// GPU MEMORY POOL +// ============================================================================ + +// ---------------------------------------------------------------------------- +// GPU Pool Initialization & Warmup +// ---------------------------------------------------------------------------- + +/** + * @brief Initialize the GPU memory pool + * + * Configures the default CUDA memory pool with a 5GB retention limit + * and pre-warms the pool with common allocation sizes. + */ +void gpuPoolInit(); + +/** + * @brief Pre-warm the GPU memory pool with common allocation sizes + * + * Allocates and immediately frees blocks to populate the pool's free-list. + * + * @param sizes Array of allocation sizes in bytes + * @param counts Array of how many allocations per size + * @param numSizes Number of different sizes to pre-allocate + * @param stream CUDA stream to use for allocations (use 0 for default) + * @return cudaError_t from allocation operations + */ +cudaError_t gpuPoolWarmup(const size_t* sizes, const int* counts, int numSizes, cudaStream_t stream); + +// ---------------------------------------------------------------------------- +// GPU Pool Raw Operations (Internal - use macros for public API) +// ---------------------------------------------------------------------------- + +/** + * @brief Allocate GPU memory from the pool (raw, internal, stream-aware) + * + * @param devPtr Output pointer to allocated device memory + * @param size Size of memory to allocate in bytes + * @param stream CUDA stream for ordering (NULL for synchronous) + * @return cudaError_t result from allocation + */ +cudaError_t gpuPoolMallocRaw(void** devPtr, size_t size, cudaStream_t stream); + +/** + * @brief Free GPU memory back to the pool (raw, internal, stream-aware) + * + * @param ptr Device pointer to free + * @param stream CUDA stream for ordering (NULL for synchronous) + * @return cudaError_t result from free operation + */ +cudaError_t gpuPoolFreeRaw(void* ptr, cudaStream_t stream); + +// ---------------------------------------------------------------------------- +// GPU Pool API with Logging - Pool-based Allocation +// ---------------------------------------------------------------------------- + +/** + * @brief GPU memory allocation from pool with logging (use gpuPoolMalloc macro) + * + * @param devPtr Output pointer to allocated device memory + * @param size Size to allocate in bytes + * @param stream CUDA stream for ordering + * @param tag Pointer identifier tag + * @param functionTag Function context tag + * @param file Source file + * @param line Line number + * @return cudaError_t from allocation operation + */ +cudaError_t gpuPoolMallocImpl(void** devPtr, size_t size, cudaStream_t stream, + const char* tag, const char* functionTag, + const char* file, int line); + +/** + * @brief GPU memory deallocation to pool with logging (use gpuPoolFree macro) + * + * @param devPtr Device pointer to free + * @param stream CUDA stream for ordering + * @param tag Pointer identifier tag + * @param functionTag Function context tag + * @param file Source file + * @param line Line number + * @return cudaError_t from free operation + */ +cudaError_t gpuPoolFreeImpl(void* devPtr, cudaStream_t stream, + const char* tag, const char* functionTag, + const char* file, int line); + +// ---------------------------------------------------------------------------- +// GPU API with Logging - Direct Allocation (Bypass Pool) +// ---------------------------------------------------------------------------- + +/** + * @brief Direct GPU memory allocation with logging (use gpuMalloc macro) + * + * Bypasses the pool and allocates directly using cudaMalloc. + * + * @param devPtr Output pointer to allocated device memory + * @param size Size to allocate in bytes + * @param stream CUDA stream (currently unused, reserved for future use) + * @param tag Pointer identifier tag + * @param functionTag Function context tag + * @param file Source file + * @param line Line number + * @return cudaError_t from allocation operation + */ +cudaError_t gpuMallocImpl(void** devPtr, size_t size, cudaStream_t stream, + const char* tag, const char* functionTag, + const char* file, int line); + +/** + * @brief Direct GPU memory deallocation with logging (use gpuFree macro) + * + * Bypasses the pool and frees directly using cudaFree. + * + * @param devPtr Device pointer to free + * @param stream CUDA stream (currently unused, reserved for future use) + * @param tag Pointer identifier tag + * @param functionTag Function context tag + * @param file Source file + * @param line Line number + * @return cudaError_t from free operation + */ +cudaError_t gpuFreeImpl(void* devPtr, cudaStream_t stream, + const char* tag, const char* functionTag, + const char* file, int line); + +// ---------------------------------------------------------------------------- +// GPU Pool Type-Safe Wrappers (Templates) +// ---------------------------------------------------------------------------- + +/** + * @brief Type-safe wrapper for GPU pool allocation + */ +template +inline cudaError_t gpuPoolMallocTyped(T** ptr, size_t size, cudaStream_t stream, + const char* pointerIdTag, const char* functionTag, + const char* file, int line) { + return gpuPoolMallocImpl(reinterpret_cast(ptr), size, stream, + pointerIdTag, functionTag, file, line); +} + +/** + * @brief Type-safe wrapper for direct GPU allocation (bypass pool) + */ +template +inline cudaError_t gpuMallocTyped(T** ptr, size_t size, cudaStream_t stream, + const char* pointerIdTag, const char* functionTag, + const char* file, int line) { + return gpuMallocImpl(reinterpret_cast(ptr), size, stream, + pointerIdTag, functionTag, file, line); +} + +// ---------------------------------------------------------------------------- +// GPU Pool Public API Macros +// ---------------------------------------------------------------------------- + +// Pool-based allocation +#define gpuPoolMalloc(ptr, size, stream, funcTag) \ + gpuPoolMallocTyped(ptr, size, stream, #ptr, funcTag, __FILE__, __LINE__) + +#define gpuPoolFree(ptr, stream, funcTag) \ + gpuPoolFreeImpl(ptr, stream, #ptr, funcTag, __FILE__, __LINE__) + +// Direct allocation (bypasses pool) +#define gpuMalloc(ptr, size, stream, funcTag) \ + gpuMallocTyped(ptr, size, stream, #ptr, funcTag, __FILE__, __LINE__) + +#define gpuFree(ptr, stream, funcTag) \ + gpuFreeImpl(ptr, stream, #ptr, funcTag, __FILE__, __LINE__) + +// ============================================================================ +// Host Memory Pool (Pinned Memory) +// ============================================================================ + +// ---------------------------------------------------------------------------- +// Host Pool Initialization & Warmup +// ---------------------------------------------------------------------------- + +/** + * @brief Initialize the host memory pool + * + * Creates the pool lock and pre-warms the pool with common allocation sizes + * ranging from 256B to 64MB across 19 bucket sizes. + */ +void hostPoolInit(); + +// ---------------------------------------------------------------------------- +// Host Pool Raw Operations (Internal - use macros for public API) +// ---------------------------------------------------------------------------- + +/** + * @brief Allocate host memory from pool (raw, internal) + * + * @param ptr Output pointer to allocated memory + * @param size Size to allocate in bytes + * @return cudaError_t from allocation operation + */ +cudaError_t hostPoolMallocRaw(void** ptr, size_t size); + +/** + * @brief Free host memory back to pool (raw, internal) + * + * @param ptr Pointer to free + * @return cudaError_t from operation + */ +cudaError_t hostPoolFreeRaw(void* ptr); + +// ---------------------------------------------------------------------------- +// Host Pool Maintenance & Trimming +// ---------------------------------------------------------------------------- + +/** + * @brief Trim host pool to target capacity + * + * Frees excess blocks from the pool to the OS, prioritizing large buckets. + * + * @param targetCapacityGB Target maximum pool capacity in GB + * @param minCapacityPerBucketMB Minimum capacity to preserve per active bucket in MB + * @return Total bytes freed to the OS + */ +size_t hostPoolTrim(double targetCapacityGB, size_t minCapacityPerBucketMB); + +/** + * @brief Adaptive refill of hot buckets between timesteps + * + * Refills buckets that had misses during the previous timestep. + */ +void hostPoolAdaptiveRefill(); + +// ---------------------------------------------------------------------------- +// Host Pool Analytics & Reporting +// ---------------------------------------------------------------------------- + +/** + * @brief Report host pool statistics with enhanced analytics + * + * @param prefix String prefix for output lines (for formatting) + * @param targetCapacityGB Target pool capacity in GB (for trim reporting) + */ +void hostPoolReportStats(const char* prefix, double targetCapacityGB); + +/** + * @brief Analyze warmup effectiveness by comparing targets vs. actual usage + */ +void hostPoolAnalyzeWarmupEffectiveness(); + +/** + * @brief Analyze memory growth trends over time + */ +void hostPoolAnalyzeGrowth(); + +// ---------------------------------------------------------------------------- +// Host Pool API with Logging - Pool-based Allocation +// ---------------------------------------------------------------------------- + +/** + * @brief Host memory allocation from pool with logging (use hostPoolMalloc macro) + * + * @param ptr Output pointer to allocated memory + * @param size Size to allocate in bytes + * @param tag Pointer identifier tag + * @param functionTag Function context tag + * @param file Source file + * @param line Line number + * @return cudaError_t from allocation operation + */ +cudaError_t hostPoolMallocImpl(void** ptr, size_t size, + const char* tag, const char* functionTag, + const char* file, int line); + +/** + * @brief Host memory deallocation to pool with logging (use hostPoolFree macro) + * + * @param ptr Pointer to free + * @param tag Pointer identifier tag + * @param functionTag Function context tag + * @param file Source file + * @param line Line number + * @return cudaError_t from free operation + */ +cudaError_t hostPoolFreeImpl(void* ptr, + const char* tag, const char* functionTag, + const char* file, int line); + +// ---------------------------------------------------------------------------- +// Host API with Logging - Direct Allocation (Bypass Pool) +// ---------------------------------------------------------------------------- + +/** + * @brief Direct host memory allocation with logging (use hostMalloc macro) + * + * Bypasses the pool and allocates directly using cudaMallocHost. + * + * @param ptr Output pointer to allocated memory + * @param size Size to allocate in bytes + * @param tag Pointer identifier tag + * @param functionTag Function context tag + * @param file Source file + * @param line Line number + * @return cudaError_t from allocation operation + */ +cudaError_t hostMallocImpl(void** ptr, size_t size, + const char* tag, const char* functionTag, + const char* file, int line); + +/** + * @brief Direct host memory deallocation with logging (use hostFree macro) + * + * Bypasses the pool and frees directly using cudaFreeHost. + * + * @param ptr Pointer to free + * @param tag Pointer identifier tag + * @param functionTag Function context tag + * @param file Source file + * @param line Line number + * @return cudaError_t from free operation + */ +cudaError_t hostFreeImpl(void* ptr, + const char* tag, const char* functionTag, + const char* file, int line); + +// ---------------------------------------------------------------------------- +// Host Pool Type-Safe Wrappers (Templates) +// ---------------------------------------------------------------------------- + +/** + * @brief Type-safe wrapper for host pool allocation + */ +template +inline cudaError_t hostPoolMallocTyped(T** ptr, size_t size, + const char* pointerIdTag, const char* functionTag, + const char* file, int line) { + return hostPoolMallocImpl(reinterpret_cast(ptr), size, + pointerIdTag, functionTag, file, line); +} + +/** + * @brief Type-safe wrapper for direct host allocation (bypass pool) + */ +template +inline cudaError_t hostMallocTyped(T** ptr, size_t size, + const char* pointerIdTag, const char* functionTag, + const char* file, int line) { + return hostMallocImpl(reinterpret_cast(ptr), size, + pointerIdTag, functionTag, file, line); +} + +// ---------------------------------------------------------------------------- +// Host Pool Public API Macros +// ---------------------------------------------------------------------------- + +// Pool-based allocation +#define hostPoolMalloc(ptr, size, funcTag) \ + hostPoolMallocTyped(ptr, size, #ptr, funcTag, __FILE__, __LINE__) + +#define hostPoolFree(ptr, funcTag) \ + hostPoolFreeImpl(ptr, #ptr, funcTag, __FILE__, __LINE__) + +// Direct allocation (bypasses pool) +#define hostMalloc(ptr, size, funcTag) \ + hostMallocTyped(ptr, size, #ptr, funcTag, __FILE__, __LINE__) + +#define hostFree(ptr, funcTag) \ + hostFreeImpl(ptr, #ptr, funcTag, __FILE__, __LINE__) + +#endif // CUDA + +#endif // MEMORY_POOL_H diff --git a/PEList.cpp b/PEList.cpp new file mode 100644 index 00000000..974a65dc --- /dev/null +++ b/PEList.cpp @@ -0,0 +1,234 @@ +#ifdef CUDA + +#include "PEList.h" +#include "DataManager.h" +#include "HostCUDA.h" + +/// @brief Each TreePiece on a given PE checks in as its tree walk completes +/// Once all TreePieces are done, launch a gravity kernel on the GPU +/// @param treePiece A reference to the TreePiece that checked in +void PEList::finishWalk() { + nCompleted++; + + // On first call, find the total number of active pieces on this PE. + // The charm++ location manager gives us this count in cTreePieces + if(nCompleted == 1) { + CkLocMgr *locMgr = treeProxy.ckLocMgr(); + locMgr->iterate(cTreePieces); + } + + // check if we have everyone + if(nCompleted < cTreePieces.count) + return; + + CkAssert(nCompleted == cTreePieces.count); // a treepiece + // checked in more + // than once? + + // bucketMarkers[i+1] is needed to determine # of IL entries per bucket + if(finalBucketMarker != -1) + bucketMarkers.push_back(finalBucketMarker); + + finishCb = new CkCallback(CkIndex_PEList::finishWalkCb(), CkMyPe(), thisProxy); + + // If the DataManager device pointer is NULL, the GPU data transfer is + // still in progress and we need to delay the kernel launch. + // launchKernel() always uses d_localMoments/d_localParts; remote also needs remote data. + DataManager *dm = dMProxy.ckLocalBranch(); + bool dataReady = dm->bLocalDataTransferred.load() && + (!bRemote || dm->bRemoteDataTransferred.load()); + if (!dataReady) + bKernelDelayed = 1; + else + launchKernel(); +} + +/// @brief Called from DataManager when remote or local transfer finishes. +/// Launch our kernel only if it was delayed AND the data we need is now ready. +/// launchKernel() always uses d_localMoments/d_localParts; remote also needs remote data. +void PEList::tryLaunchDelayedKernel() { + if (!bKernelDelayed) + return; + DataManager *dm = dMProxy.ckLocalBranch(); + bool dataReady = dm->bLocalDataTransferred.load() && + (!bRemote || dm->bRemoteDataTransferred.load()); + if (!dataReady) + return; + bKernelDelayed = 0; + launchKernel(); +} + +void PEList::finishWalkCb() { + dMProxy.ckLocalBranch()->transferParticleVarsBack(false); + reset(); +} + +/// @brief Launch the corresponding CUDA kernel, depending what type of request this was +void PEList::launchKernel() { + request = new CudaRequest; + + // Ensure required data is present on the GPU + CkAssert(dMProxy.ckLocalBranch()->d_localMoments != nullptr); + CkAssert(dMProxy.ckLocalBranch()->d_localParts != nullptr); + CkAssert(dMProxy.ckLocalBranch()->d_localVars != nullptr); + // The following checks can fail if remote prefetch does not bring in any data. + // if (bRemote) { + // CkAssert(dMProxy.ckLocalBranch()->d_remoteParts != nullptr); + // CkAssert(dMProxy.ckLocalBranch()->d_remoteMoments != nullptr); + // } + + request->d_localMoments = dMProxy.ckLocalBranch()->d_localMoments; + request->d_localParts = dMProxy.ckLocalBranch()->d_localParts; + request->d_localVars = dMProxy.ckLocalBranch()->d_localVars; + request->d_remoteParts = dMProxy.ckLocalBranch()->d_remoteParts; + request->d_remoteMoments = dMProxy.ckLocalBranch()->d_remoteMoments; + request->stream = stream; + + request->numBucketsPlusOne = bucketSizes.size()+1; + + request->node = bNode; + request->remote = bRemote; + + request->fperiod = fperiod; + + const char* funcTag = "PEList::finish"; + if (iList.size() > 0) { + hostMalloc((void**)&iListHost, iList.size()*sizeof(ILCell), funcTag); + memcpy(iListHost, iList.data(), iList.size()*sizeof(ILCell)); + hostMalloc((void**)&bucketMarkersHost, bucketMarkers.size()*sizeof(int), funcTag); + memcpy(bucketMarkersHost, bucketMarkers.data(), bucketMarkers.size()*sizeof(int)); + hostMalloc((void**)&bucketStartsHost, bucketStarts.size()*sizeof(int), funcTag); + memcpy(bucketStartsHost, bucketStarts.data(), bucketStarts.size()*sizeof(int)); + hostMalloc((void**)&bucketSizesHost, bucketSizes.size()*sizeof(int), funcTag); + memcpy(bucketSizesHost, bucketSizes.data(), bucketSizes.size()*sizeof(int)); + } + + if (missedParts.size() > 0) { + hostMalloc((void**)&missedPartsHost, missedParts.size()*sizeof(CompactPartData), funcTag); + memcpy(missedPartsHost, missedParts.data(), missedParts.size()*sizeof(CompactPartData)); + } + if (missedNodes.size() > 0) { + hostMalloc((void**)&missedNodesHost, missedNodes.size()*sizeof(CudaMultipoleMoments), funcTag); + memcpy(missedNodesHost, missedNodes.data(), missedNodes.size()*sizeof(CudaMultipoleMoments)); + } + + request->list = iListHost; + request->missedParts = missedPartsHost; + request->missedNodes = missedNodesHost; + request->sMissed = bNode ? missedNodes.size()*sizeof(CudaMultipoleMoments) : missedParts.size()*sizeof(CompactPartData); + request->bucketMarkers = bucketMarkersHost; + request->bucketStarts = bucketStartsHost; + request->bucketSizes = bucketSizesHost; + request->numInteractions = iList.size(); + + request->cb = finishCb; + + void (*transferFunc)(CudaRequest*); + if (bNode) { + transferFunc = bRemote ? PEListNodeListDataTransferRemote : PEListNodeListDataTransferLocal; + if (bResume) { + transferFunc = PEListNodeListDataTransferRemoteResume; + } + } else { + transferFunc = bRemote ? PEListPartListDataTransferRemote : PEListPartListDataTransferLocal; + if (bResume) { + transferFunc = PEListPartListDataTransferRemoteResume; + } + } + + transferFunc(request); +} + +/// @brief Collect the interaction list results each a Compute operation completes +/// @param treePiece The TreePiece that sent the operation +/// @data A CudaRequest object containing the interaction list data +void PEList::sendList(TreePiece *treePiece, CudaRequest* data) { + int numBucketsPlusOne = data->numBucketsPlusOne; + int numBuckets = numBucketsPlusOne-1; + + // bucketMarkers need an offset because we are concatenating the interaction lists + for (int i = 0; i < numBuckets; i++) { + bucketMarkers.push_back(data->bucketMarkers[i] + iList.size()); + } + finalBucketMarker = data->bucketMarkers[numBuckets] + iList.size(); + + for (int i = 0; i < numBuckets; i++) { + bucketStarts.push_back(data->bucketStarts[i]); + bucketSizes.push_back(data->bucketSizes[i]); + } + + // If we have missed parts/nodes, the indices in the interaction list + // need to be shifted because the remote data is being concatenated + if (data->missedParts) { + // Note that many TreePieces will have the same missed particles + // We are copying a lot of duplicate data to the GPU here + int numMissedParts = data->sMissed/sizeof(CompactPartData); + int missedOffset = missedParts.size(); + for (int i = 0; i < data->numInteractions; i++) { + ((ILCell *)data->list)[i].index += missedOffset; + iList.push_back(((ILCell *)data->list)[i]); + } + for (int i = 0; i < numMissedParts; i++) { + missedParts.push_back(((CompactPartData *)data->missedParts)[i]); + } + } else if (data->missedNodes) { + int numMissedNodes = data->sMissed/sizeof(CudaMultipoleMoments); + int missedOffset = missedNodes.size(); + for (int i = 0; i < data->numInteractions; i++) { + ((ILCell *)data->list)[i].index += missedOffset; + iList.push_back(((ILCell *)data->list)[i]); + } + for (int i = 0; i < numMissedNodes; i++) { + missedNodes.push_back(((CudaMultipoleMoments *)data->missedNodes)[i]); + } + } else { + for (int i = 0; i < data->numInteractions; i++) { + iList.push_back(((ILCell *)data->list)[i]); + } + } + + fperiod = data->fperiod; + + // Call finishBucket for all buckets involved in this interaction + treePiece->cudaFinishAffectedBuckets(data->affectedBuckets, numBuckets, bRemote); + + // deallocate the memory used by the incoming cudaRequest + free(data->list); + free(data->bucketMarkers); + free(data->bucketStarts); + free(data->bucketSizes); + delete[] data->affectedBuckets; +} + +/// @brief Re-initalize data arrays and clean up callback objects at the end of the step +void PEList::reset() { + const char* funcTag = "PEList::reset"; + if (iList.size() > 0) { + hostFree(iListHost, funcTag); + hostFree(bucketMarkersHost, funcTag); + hostFree(bucketStartsHost, funcTag); + hostFree(bucketSizesHost, funcTag); + } + if (missedParts.size() > 0) { + hostFree(missedPartsHost, funcTag); + } + if (missedNodes.size() > 0) { + hostFree(missedNodesHost, funcTag); + } + iList.clear(); + missedParts.clear(); + missedNodes.clear(); + bucketMarkers.clear(); + bucketStarts.clear(); + bucketSizes.clear(); + + cTreePieces.reset(); + nCompleted = 0; + bKernelDelayed = 0; + finalBucketMarker = -1; + delete finishCb; + delete request; + request = nullptr; +} + +#endif diff --git a/PEList.h b/PEList.h new file mode 100644 index 00000000..76e1093c --- /dev/null +++ b/PEList.h @@ -0,0 +1,96 @@ +#ifndef PE_LIST_H +#define PE_LIST_H + +#ifdef CUDA +#include "ParallelGravity.h" + +/// @brief PE-level group to coordinate gravity interaction list calculations +/// on the GPU +/// +/// There is one instance of this class for each type of walk (particle/node, +/// local, remote, remote-resume). Walks are done at the TreePiece level and +/// the data from each walk is consolidated via the sendList function. Once +/// all of the walks have finished, the finishWalk routine moves the +/// interaction lists to the GPU and starts the computation. +/// +/// When GPU_LOCAL_TREE_WALK is enabled, only the remote and remote-resume +/// walks are handled here. +class PEList : public CBase_PEList +{ + /// non-empty TreePieces that have completed their walk on this PE + int nCompleted; + /// Count of TreePieces with particles on this PE + NonEmptyTreePieceCounter cTreePieces; + + vector iList; + + vector bucketMarkers; + int finalBucketMarker; + vector bucketStarts; + vector bucketSizes; + + vector missedParts; + vector missedNodes; + + ILCell* iListHost; + CompactPartData* missedPartsHost; + CudaMultipoleMoments* missedNodesHost; + int *bucketMarkersHost; + int *bucketStartsHost; + int *bucketSizesHost; + + CudaRequest *request; + + /// Type of request + int bNode; + int bRemote; + int bResume; + + /// Flags whether the GPU kernel launch was delayed due to the data transfer + int bKernelDelayed; + + cudaStream_t stream; + CkCallback *finishCb; + + cudatype fperiod; + + public: + PEList(int _bNode, int _bRemote, int _bResume) { + bNode = _bNode; + bRemote = _bRemote; + bResume = _bResume; + + nCompleted = 0; + bKernelDelayed = 0; + + finalBucketMarker = -1; + request = nullptr; + cudaStreamCreate(&stream); + } + PEList(CkMigrateMessage *m) : CBase_PEList(m) { + nCompleted = 0; + bKernelDelayed = 0; + + finalBucketMarker = -1; + request = nullptr; + cudaStreamCreate(&stream); + } + ~PEList() { cudaStreamDestroy(stream); } + void pup(PUP::er &p) { + CBase_PEList::pup(p); + + p|bNode; + p|bRemote; + p|bResume; + } + + void finishWalk(); + void launchKernel(); + void finishWalkCb(); + void tryLaunchDelayedKernel(); + void sendList(TreePiece *treePiece, CudaRequest *data); + void reset(); +}; + +#endif +#endif diff --git a/ParallelGravity.ci b/ParallelGravity.ci index c3870824..8635783a 100644 --- a/ParallelGravity.ci +++ b/ParallelGravity.ci @@ -29,6 +29,15 @@ mainmodule ParallelGravity { readonly CProxy_DumpFrameData dfDataProxy; readonly CProxy_IntraNodeLBManager nodeLBMgrProxy; +#ifdef CUDA + readonly CProxy_PEList peNodeLocalListProxy; + readonly CProxy_PEList peNodeRemoteListProxy; + readonly CProxy_PEList peNodeRemoteResumeListProxy; + readonly CProxy_PEList pePartLocalListProxy; + readonly CProxy_PEList pePartRemoteListProxy; + readonly CProxy_PEList pePartRemoteResumeListProxy; +#endif + readonly bool _cache; readonly int _nocache; readonly int _cacheLineDepth; @@ -56,7 +65,6 @@ mainmodule ParallelGravity { readonly int localPartsPerReq; readonly int remotePartsPerReq; readonly int remoteResumePartsPerReq; - readonly double largePhaseThreshold; readonly int START_REG; readonly int START_IB; @@ -91,9 +99,12 @@ mainmodule ParallelGravity { message dummyMsg; message ComputeChunkMsg; - message EwaldMsg{ - bool fromInit; +#ifdef CUDA + message fillGPUMsg{ + int partIndex; + int nParts; }; +#endif message ORBSplittersMsg{ double pos[]; @@ -261,8 +272,19 @@ mainmodule ParallelGravity { entry void combineLocalTrees(CkReductionMsg *m); #ifdef CUDA entry void startLocalWalk(); +#ifdef GPU_LOCAL_TREE_WALK + entry void finishLocalWalk(); +#endif + entry void startEwaldGPU(); + entry void finishEwaldGPU(); entry void resumeRemoteChunk(); - entry void createStreams(int _numStreams, const CkCallback& cb); + entry void createStream(const CkCallback& cb); + entry void initMemLog(std::string _fileName, int bGpuMemLoggerFlag, const CkCallback &cb); + entry void initCpuMemLog(std::string _fileName, int bCpuMemLoggerFlag, const CkCallback &cb); + entry void flushMemLog(const CkCallback& cb); + entry void flushCpuMemLog(const CkCallback& cb); + entry void trimHostPool(double targetCapacityGB, const CkCallback& cb); + entry void refillHostPool(const CkCallback& cb); #endif entry void initCooling(double dGmPerCcUnit, double dComovingGmPerCcUnit, double dErgPerGmUnit, double dSecUnit, double dKpcUnit, @@ -295,12 +317,9 @@ mainmodule ParallelGravity { entry void setPeriodic(int nReplicas, Vector3D fPeriod, int bEwald, double fEwCut, double fEwhCut, int bPeriod, int bComove, double dRhoFac); - entry [notrace] void EwaldInit(); entry [notrace] void initCoolingData(const CkCallback& cb); entry [notrace] void initLWData(const CkCallback& cb); - entry void calculateEwald(EwaldMsg *m); - entry [notrace] void EwaldGPUComplete(); - entry void EwaldGPU(); + entry void calculateEwald(dummyMsg *msg); entry void velScale(double dScale, const CkCallback& cb); //entry void markPresence(const CkCallback& cb); entry void loadNChilada(const std::string& filename, const double dTuFac, @@ -485,16 +504,10 @@ mainmodule ParallelGravity { OrientedBox& box, OrientedBox& boxBall, unsigned int iParticleTypes, int64_t nSPH); -#ifdef CUDA - entry void commenceCalculateGravityLocal(intptr_t d_localMoments, - intptr_t d_localParts, - intptr_t d_localVars, - intptr_t streams, int numStreams, - size_t sMoments, size_t sCompactParts, size_t sVarParts); -#else + entry void commenceCalculateGravityLocal(); -#endif entry void calculateGravityRemote(ComputeChunkMsg *msg); + entry void cudaFinishAllBuckets(int fromEwald); entry [notrace] void calculateReSmoothLocal(); entry [notrace] void calculateMarkSmoothLocal(); @@ -523,14 +536,11 @@ mainmodule ParallelGravity { // jetley #ifdef CUDA - entry void continueStartRemoteChunk(int chunk, intptr_t d_remoteMoments, intptr_t d_remoteParts); - entry void fillGPUBuffer(intptr_t bufLocalParts, - intptr_t bufLocalMoments, - intptr_t pLocalMoments, int partIndex, int nParts, intptr_t node); + entry void fillGPUBuffer(fillGPUMsg* msg); entry void updateParticles(intptr_t data, int partIndex); -#else - entry void continueStartRemoteChunk(int chunk); + entry void flushInteractionsToGpu(intptr_t ptr_lc, intptr_t ptr_ds); #endif + entry void continueStartRemoteChunk(int chunk); entry void continueWrapUp(); // jetley - balancing before initial forces @@ -618,6 +628,14 @@ mainmodule ParallelGravity { entry [local] void finishedTPWork(); }; +#ifdef CUDA + group [migratable] PEList { + entry PEList(int _bNode, int _bRemote, int _bResume); + entry void finishWalkCb(); + entry void tryLaunchDelayedKernel(); + }; +#endif + initproc void registerStatistics(); diff --git a/ParallelGravity.cpp b/ParallelGravity.cpp index 420c791e..4c62c13f 100644 --- a/ParallelGravity.cpp +++ b/ParallelGravity.cpp @@ -50,8 +50,10 @@ #endif #ifdef CUDA +#include "PEList.h" // for default per-list parameters #include "cuda_typedef.h" +#include "MemoryPool.h" #endif extern char *optarg; @@ -89,7 +91,26 @@ CProxy_DumpFrameData dfDataProxy; /// @brief Proxy for the PETreeMerger group. CProxy_PETreeMerger peTreeMergerProxy; - +#ifdef CUDA +/// @brief Proxies for the PE groups that ship interactions to the GPU +CProxy_PEList peNodeLocalListProxy; +CProxy_PEList peNodeRemoteListProxy; +CProxy_PEList peNodeRemoteResumeListProxy; +CProxy_PEList pePartLocalListProxy; +CProxy_PEList pePartRemoteListProxy; +CProxy_PEList pePartRemoteResumeListProxy; + +/// @brief Array for quickly iterating over interaction list proxies +CProxy_PEList* PEListProxies[] = { + &peNodeLocalListProxy, + &peNodeRemoteListProxy, + &peNodeRemoteResumeListProxy, + &pePartLocalListProxy, + &pePartRemoteListProxy, + &pePartRemoteResumeListProxy +}; +const int numPEListProxies = sizeof(PEListProxies) / sizeof(PEListProxies[0]); +#endif /// @brief Use the cache (always on) bool _cache; @@ -115,10 +136,6 @@ GenericTrees useTree; CProxy_TreePiece streamingProxy; /// @brief Number of pieces into which to divide the tree. unsigned int numTreePieces; -#ifdef CUDA - /// @brief Number of CUDA streams to use - unsigned int numStreams; -#endif /// @brief Number of particles per TreePiece. Used to determine the /// number of TreePieces. unsigned int particlesPerChare; @@ -141,10 +158,6 @@ int localPartsPerReq; int remotePartsPerReq; int remoteResumePartsPerReq; -// multi-stepping particle transfer strategy -// switch threshold -double largePhaseThreshold; - cosmoType theta; ///< BH-like opening criterion cosmoType thetaMono; ///< Criterion of excepting monopole /// only cells. @@ -255,12 +268,13 @@ Main::Main(CkArgMsg* m) { traceRegisterUserEvent("Xfer Local", CUDA_XFER_LOCAL); traceRegisterUserEvent("Xfer Remote", CUDA_XFER_REMOTE); - traceRegisterUserEvent("Grav Local", CUDA_GRAV_LOCAL); - traceRegisterUserEvent("Grav Remote", CUDA_GRAV_REMOTE); - traceRegisterUserEvent("Remote Resume", CUDA_REMOTE_RESUME); - traceRegisterUserEvent("Part Gravity Local", CUDA_PART_GRAV_LOCAL); - traceRegisterUserEvent("Part Gravity Local Small", CUDA_PART_GRAV_LOCAL_SMALL); - traceRegisterUserEvent("Part Gravity Remote", CUDA_PART_GRAV_REMOTE); + traceRegisterUserEvent("Grav Tree Local", CUDA_GRAV_TREE_LOCAL); + traceRegisterUserEvent("Grav Nodelist Local", CUDA_GRAV_NODELIST_LOCAL); + traceRegisterUserEvent("Grav Partlist Local", CUDA_GRAV_PARTLIST_LOCAL); + traceRegisterUserEvent("Grav Nodelist Remote", CUDA_GRAV_NODELIST_REMOTE); + traceRegisterUserEvent("Grav Partlist Remote", CUDA_GRAV_PARTLIST_REMOTE); + traceRegisterUserEvent("Grav Nodelist Remote Resume", CUDA_GRAV_NODELIST_REMOTE_RESUME); + traceRegisterUserEvent("Grav Partlist Remote Resume", CUDA_GRAV_PARTLIST_REMOTE_RESUME); traceRegisterUserEvent("Xfer Back", CUDA_XFER_BACK); traceRegisterUserEvent("Ewald", CUDA_EWALD); #endif @@ -772,12 +786,31 @@ Main::Main(CkArgMsg* m) { prmAddParam(prm, "nTreePieces", paramInt, &numTreePieces, sizeof(int),"p", "Number of TreePieces (default: 8*procs)"); #ifdef CUDA - numStreams = 100; - prmAddParam(prm, "nStreams", paramInt, &numStreams, - sizeof(int),"str", "Number of CUDA streams (default: 100)"); param.nGpuMinParts = 1000; prmAddParam(prm, "nGpuMinParts", paramInt, ¶m.nGpuMinParts, sizeof(int),"gpup", "Min particles on rung to trigger GPU (default: 1000)"); + param.bGpuMemLogger = 0; + prmAddParam(prm, "bGpuMemLogger", paramBool, ¶m.bGpuMemLogger, + sizeof(int), "gpumemlog", + "Enable GPU memory pool logger"); + param.bCpuMemLogger = 0; + prmAddParam(prm, "bCpuMemLogger", paramBool, ¶m.bCpuMemLogger, + sizeof(int), "cpumemlog", + "Enable CPU memory pool logger"); + param.dHostPoolTargetCapacityGB = 5.0; + prmAddParam(prm, "dHostPoolTargetCapacityGB", paramDouble, + ¶m.dHostPoolTargetCapacityGB, + sizeof(double), "hptarget", + "Target capacity for host pool (per DM) in GB (default: 5.0)"); + param.nHostPoolMinCapacityPerBucketMB = 50; + prmAddParam(prm, "nHostPoolMinCapacityPerBucketMB", paramInt, + ¶m.nHostPoolMinCapacityPerBucketMB, + sizeof(int), "hpminbucket", + "Min capacity per host pool bucket size during trim in MB (default: 50)"); + param.bHostPoolDebug = 0; + prmAddParam(prm, "bHostPoolDebug", paramBool, ¶m.bHostPoolDebug, + sizeof(int), "hpdebug", + "Enable host pool debug output and statistics"); #endif particlesPerChare = 0; prmAddParam(prm, "nPartPerChare", paramInt, &particlesPerChare, @@ -910,10 +943,6 @@ Main::Main(CkArgMsg* m) { prmAddParam(prm, "remoteResumePartsPerReq", paramDouble, &remoteResumePartsPerReqDouble, sizeof(double),"remoteresumeparts", "Num. remote resume particle interactions allowed per CUDA request (in millions)"); - largePhaseThreshold = TP_LARGE_PHASE_THRESHOLD_DEFAULT; -// prmAddParam(prm, "largePhaseThreshold", paramDouble, &largePhaseThreshold, -// sizeof(double),"largephasethresh", "Ratio of active to total particles at which all particles (not just active ones) are sent to gpu in the target buffer (No source particles are sent.)"); - #endif int processSimfile = 1; @@ -1118,9 +1147,6 @@ Main::Main(CkArgMsg* m) { ckout << "remoteParts: " << remotePartsPerReqDouble << endl; ckout << "INFO: "; ckout << "remoteResumeParts: " << remoteResumePartsPerReqDouble << endl; - - ckout << "INFO: largePhaseThreshold: " << largePhaseThreshold << endl; - #endif if(prmSpecified(prm, "bGeometric")) { @@ -1392,6 +1418,15 @@ Main::Main(CkArgMsg* m) { peTreeMergerProxy = CProxy_PETreeMerger::ckNew(); dfDataProxy = CProxy_DumpFrameData::ckNew(); + +#ifdef CUDA + peNodeLocalListProxy = CProxy_PEList::ckNew(1, 0, 0); + peNodeRemoteListProxy = CProxy_PEList::ckNew(1, 1, 0); + peNodeRemoteResumeListProxy = CProxy_PEList::ckNew(1, 1, 1); + pePartLocalListProxy = CProxy_PEList::ckNew(0, 0, 0); + pePartRemoteListProxy = CProxy_PEList::ckNew(0, 1, 0); + pePartRemoteResumeListProxy = CProxy_PEList::ckNew(0, 1, 1); +#endif // create CacheManagers // Gravity particles @@ -1992,6 +2027,7 @@ void Main::buildTree(int iPhase) #else treeProxy.buildTree(bucketSize, CkCallbackResumeThread()); #endif + double tTB = CkWallTimer()-startTime; timings[iPhase].tTBuild += tTB; CkPrintf("took %g seconds.\n", tTB); @@ -2473,8 +2509,9 @@ void Main::setupICs() { double startTime; #ifdef CUDA - dMProxy.createStreams(numStreams, CkCallbackResumeThread()); + dMProxy.createStream(CkCallbackResumeThread()); #endif + treeProxy.setPeriodic(param.nReplicas, param.vPeriod, param.bEwald, param.dEwCut, param.dEwhCut, param.bPeriodic, param.csm->bComove, @@ -2593,6 +2630,13 @@ void Main::setupICs() { treeProxy.initRand(param.iRandomSeed, CkCallbackResumeThread()); } +#ifdef CUDA + if(param.bGpuMemLogger) + initMemLog(); // Initialize GPU memory logging + if(param.bCpuMemLogger) + initCpuMemLog(); // Initialize CPU memory logging +#endif + if(param.bStarForm) { initStarLog(); if(param.feedback->sn.bUseStoch) @@ -2961,6 +3005,13 @@ Main::restart(CkCheckpointStatusMsg *msg) } dMProxy.resetReadOnly(param, CkCallbackResumeThread()); +#ifdef CUDA + dMProxy.createStream(CkCallbackResumeThread()); + if(param.bGpuMemLogger) + initMemLog(); // Initialize GPU memory logging + if(param.bCpuMemLogger) + initCpuMemLog(); // Initialize CPU memory logging +#endif if (bUseCkLoopPar) { CkPrintf("Using CkLoop %d\n", param.bUseCkLoopPar); } else { @@ -3144,6 +3195,37 @@ Main::doSimulation() << endl; writeTimings(iStep); +#ifdef CUDA + + // Host pool diagnostics and maintenance + if (param.bHostPoolDebug) { + auto prefix = make_formatted_string("[Step %d]", iStep); + hostPoolReportStats(prefix.c_str(), param.dHostPoolTargetCapacityGB); + + // Periodic growth analysis every 10 steps + if (iStep % 10 == 0) { + hostPoolAnalyzeGrowth(); + } + + // Early warmup effectiveness check + if (iStep == 5) { + hostPoolAnalyzeWarmupEffectiveness(); + } + } + + // Trim host pool to reclaim memory (always runs) + dMProxy.trimHostPool(param.dHostPoolTargetCapacityGB, CkCallbackResumeThread()); + + // Refill host pool hot buckets to prepare for next timestep + dMProxy.refillHostPool(CkCallbackResumeThread()); + + // Flush memory logs + if(param.bGpuMemLogger) + dMProxy[0].flushMemLog(CkCallbackResumeThread()); + if(param.bCpuMemLogger) + dMProxy[0].flushCpuMemLog(CkCallbackResumeThread()); +#endif + if(iStep%param.iOrbitOutInterval == 0) { outputBlackHoles(dTime); } diff --git a/ParallelGravity.h b/ParallelGravity.h index 4657d419..3d8b3779 100644 --- a/ParallelGravity.h +++ b/ParallelGravity.h @@ -169,6 +169,19 @@ extern int nIOProcessor; extern CProxy_DumpFrameData dfDataProxy; extern CProxy_PETreeMerger peTreeMergerProxy; + +#ifdef CUDA +extern CProxy_PEList peNodeLocalListProxy; +extern CProxy_PEList peNodeRemoteListProxy; +extern CProxy_PEList peNodeRemoteResumeListProxy; +extern CProxy_PEList pePartLocalListProxy; +extern CProxy_PEList pePartRemoteListProxy; +extern CProxy_PEList pePartRemoteResumeListProxy; + +extern CProxy_PEList* PEListProxies[]; +extern const int numPEListProxies; +#endif + extern CProxy_CkCacheManager cacheGravPart; extern CProxy_CkCacheManager cacheSmoothPart; extern CProxy_NodeCache cacheNode; @@ -202,8 +215,6 @@ extern int localPartsPerReq; extern int remotePartsPerReq; extern int remoteResumePartsPerReq; -extern double largePhaseThreshold; - extern cosmoType theta; extern cosmoType thetaMono; @@ -332,12 +343,13 @@ struct BucketMsg : public CkMcastBaseMsg, public CMessage_BucketMsg { }; #endif -/// Associated with calls to calculateEwald -/// Indicates whether the function was called by initEwald -struct EwaldMsg: public CMessage_EwaldMsg { - bool fromInit; +#ifdef CUDA +struct fillGPUMsg: public CMessage_fillGPUMsg { + int partIndex; + int nParts; }; - +#endif + /// Class to count added and deleted particles class CountSetPart { @@ -619,6 +631,8 @@ class Main : public CBase_Main { void initLWData(); void initStarLog(); void initHMStarLog(); + void initMemLog(); + void initCpuMemLog(); int ReadASCII(char *extension, int nDataPerLine, double *dDataOut); void restartGas(); void doSph(int activeRung, int bNeedDensity = 1); @@ -645,14 +659,6 @@ class Main : public CBase_Main { void restartNSIDM(); }; -/* IBM brain damage */ -#undef hz -/// @brief Coefficients for the Fourier space part of the Ewald sum. -typedef struct ewaldTable { - double hx,hy,hz; - double hCfac,hSfac; -} EWT; - // jetley class MissRecord; class State; @@ -905,6 +911,11 @@ class TreePiece : public CBase_TreePiece { /// Return the pointer to the particles on this TreePiece. GravityParticle *getParticles(){return myParticles;} + /// Access method for number of active particles on this TreePiece + int getNumActiveParticles(){ + return myNumActiveParticles; + } + #ifdef CUDA // this variable holds the number of buckets active at @@ -917,23 +928,13 @@ class TreePiece : public CBase_TreePiece { // in the list of interations to the sent to the gpu, we flush // the list int numActiveBuckets; - int myNumActiveParticles; - // First and Last indices of GPU particle - int FirstGPUParticleIndex; - int LastGPUParticleIndex; - int NumberOfGPUParticles; - BucketActiveInfo *bucketActiveInfo; - - // For accessing GPU memory - CudaMultipoleMoments *d_localMoments; - CudaMultipoleMoments *d_remoteMoments; - CompactPartData *d_localParts; - CompactPartData *d_remoteParts; - VariablePartData *d_localVars; - size_t sMoments; - size_t sCompactParts; - size_t sVarParts; - cudaStream_t stream; + /// Specifies that the Treepiece has filled the buffer with + /// particle data to be sent to the GPU. This implies that + /// this TreePiece's buckets now know where their particles + /// are in the GPU particle array. + int bGPUBufferFilled; + + int getParentPE() { return CkMyPe(); } int getNumBuckets(){ return numBuckets; @@ -944,80 +945,21 @@ class TreePiece : public CBase_TreePiece { // depending on fraction of active particles to their // total count. int getDMNumParticles(){ - if(largePhase()){ - return myNumParticles; - } - else{ - return myNumActiveParticles; - } - } - - int getNumActiveParticles(){ - return myNumActiveParticles; - } - - void calculateNumActiveParticles(){ - myNumActiveParticles = 0; - for(int i = 1; i <= myNumParticles; i++){ - if(myParticles[i].rung >= activeRung){ - myNumActiveParticles++; - } - } - } - - bool largePhase(){ - return (1.0*myNumActiveParticles/myNumParticles) >= largePhaseThreshold; + return myNumParticles; } void getDMParticles(CompactPartData *fillArray, int &fillIndex){ - NumberOfGPUParticles = 0; - FirstGPUParticleIndex = fillIndex;//This is for the GPU Ewald - if(largePhase()){ - for(int b = 0; b < numBuckets; b++){ - GenericTreeNode *bucket = bucketList[b]; - int buckstart = bucket->firstParticle; - int buckend = bucket->lastParticle; - GravityParticle *buckparts = bucket->particlePointer; - bucket->bucketArrayIndex = fillIndex; - for(int i = buckstart; i <= buckend; i++){ - fillArray[fillIndex] = buckparts[i-buckstart]; - fillIndex++; - } - } - } - else{ - for(int b = 0; b < numBuckets; b++){ - GenericTreeNode *bucket = bucketList[b]; - if(bucket->rungs < activeRung){ - continue; - } - BucketActiveInfo *binfo = &(bucketActiveInfo[b]); - - int buckstart = bucket->firstParticle; - int buckend = bucket->lastParticle; - GravityParticle *buckparts = bucket->particlePointer; - - binfo->start = fillIndex; - for(int i = buckstart; i <= buckend; i++){ - if(buckparts[i-buckstart].rung >= activeRung){ - fillArray[fillIndex] = buckparts[i-buckstart]; - fillIndex++; - } - } - binfo->size = fillIndex-binfo->start; - } - } - //This is for the GPU Ewald - if(FirstGPUParticleIndex == fillIndex){ - //This means no particle is on GPU - FirstGPUParticleIndex = -1; - LastGPUParticleIndex = -1; - NumberOfGPUParticles = 0; - } - else{ - LastGPUParticleIndex = fillIndex - 1; - NumberOfGPUParticles = LastGPUParticleIndex - FirstGPUParticleIndex + 1; - } + for(int b = 0; b < numBuckets; b++){ + GenericTreeNode *bucket = bucketList[b]; + int buckstart = bucket->firstParticle; + int buckend = bucket->lastParticle; + GravityParticle *buckparts = bucket->particlePointer; + bucket->bucketArrayIndex = fillIndex; + for(int i = buckstart; i <= buckend; i++){ + fillArray[fillIndex] = buckparts[i-buckstart]; + fillIndex++; + } + } } bool isActive(int partNum){ @@ -1039,14 +981,11 @@ class TreePiece : public CBase_TreePiece { #endif #ifdef CUDA - void continueStartRemoteChunk(int chunk, intptr_t d_remoteMoments, intptr_t d_remoteParts); - void fillGPUBuffer(intptr_t bufLocalParts, - intptr_t bufLocalMoments, - intptr_t pLocalMoments, int partIndex, int nParts, intptr_t node); + void fillGPUBuffer(fillGPUMsg *msg); void updateParticles(intptr_t data, int partIndex); -#else - void continueStartRemoteChunk(int chunk); + void flushInteractionsToGpu(intptr_t ptr_lc, intptr_t ptr_ds); #endif + void continueStartRemoteChunk(int chunk); void continueWrapUp(); #if INTERLIST_VER > 0 @@ -1104,7 +1043,9 @@ class TreePiece : public CBase_TreePiece { CkCallback after_dd_callback; /// Total number of particles contained in this chare unsigned int myNumParticles; - unsigned int numActiveParticles; + /// Number of particles on this TreePiece needing force calculations + int myNumActiveParticles; + /// Array with the particles in this chare GravityParticle* myParticles; int nbor_msgs_count_; @@ -1228,18 +1169,14 @@ class TreePiece : public CBase_TreePiece { int bComove; /// Background density of the Universe double dRhoFac; +public: Vector3D fPeriod; int nReplicas; int bEwald; /* Perform Ewald */ double fEwCut; double dEwhCut; - EWT *ewt; - int nMaxEwhLoop; - int nEwhLoop; -#ifdef HEXADECAPOLE - MOMC momcRoot; /* complete moments of root */ -#endif +private: int bGasCooling; #ifndef COOLING_NONE clDerivsData *CoolData; @@ -1356,13 +1293,6 @@ class TreePiece : public CBase_TreePiece { double totalTime; public: -#ifdef SPCUDA - EwaldData *h_idata; - CkCallback *cbEwaldGPU; -#endif - void EwaldGPU(); - void EwaldGPUComplete(); - #if COSMO_DEBUG > 1 || defined CHANGA_REFACTOR_WALKCHECK || defined CHANGA_REFACTOR_WALKCHECK_INTERLIST ///This function checks the correctness of the treewalk void checkWalkCorrectness(); @@ -1416,8 +1346,7 @@ class TreePiece : public CBase_TreePiece { /// Initialize all the buckets for the tree walk - /// @TODO: Eliminate this redundant copy! - void initBuckets(); + int initBuckets(); template void initBucketsSmooth(Tsmooth tSmooth); void smoothNextBucket(); @@ -1432,6 +1361,8 @@ class TreePiece : public CBase_TreePiece { * to trigger nextBucket() which will loop over all the buckets. */ void doAllBuckets(); + void cudaFinishAllBuckets(int fromEwald); + void cudaFinishAffectedBuckets(int *affectedBuckets, int numBuckets, int bRemote); void reconstructNodeLookup(GenericTreeNode *node); //void rebuildSFCTree(GenericTreeNode *node,GenericTreeNode *parent,int *); @@ -1490,8 +1421,6 @@ class TreePiece : public CBase_TreePiece { // temporarely set to -1, it will updated after the tree is built numChunks=-1; prefetchRoots = NULL; - ewt = NULL; - nMaxEwhLoop = 100; incomingParticlesMsg.clear(); incomingParticlesArrived = 0; @@ -1501,6 +1430,7 @@ class TreePiece : public CBase_TreePiece { mySPHParticles = NULL; myStarParticles = NULL; myNumParticles = myNumSPH = myNumStar = 0; + myNumActiveParticles = 0; nStore = nStoreSPH = nStoreStar = 0; bBucketsInited = false; myTreeParticles = -1; @@ -1530,8 +1460,6 @@ class TreePiece : public CBase_TreePiece { nPartCacheEntries = 0; completedActiveWalks = 0; prefetchRoots = NULL; - //remaining Chunk = NULL; - ewt = NULL; root = NULL; pTreeNodes = NULL; @@ -1579,7 +1507,6 @@ class TreePiece : public CBase_TreePiece { delete[] nodeInterRemote; delete[] particleInterRemote; delete[] bucketReqs; - delete[] ewt; deleteTree(); @@ -1597,10 +1524,8 @@ class TreePiece : public CBase_TreePiece { double fEwCut, double fEwhCut, int bPeriod, int bComove, double dRhoFac); void BucketEwald(GenericTreeNode *req, int nReps,double fEwCut); - void EwaldInit(); - void ewaldCPU(EwaldMsg *msg); - void calculateEwald(EwaldMsg *m); - void calculateEwaldUsingCkLoop(int yield_num); + void calculateEwald(dummyMsg *msg); + void calculateEwaldUsingCkLoop(int yield_num); void callBucketEwald(int id); void doParallelNextBucketWork(int id, LoopParData* lpdata); void initCoolingData(const CkCallback& cb); @@ -1853,15 +1778,8 @@ class TreePiece : public CBase_TreePiece { void calculateGravityLocal(); /// Do some minor preparation for the local walk then /// calculateGravityLocal(). -#ifdef CUDA - void commenceCalculateGravityLocal(intptr_t d_localMoments, - intptr_t d_localParts, - intptr_t d_localVars, - intptr_t streams, int numStreams, - size_t sMoments, size_t sCompactParts, size_t sVarParts); -#else + void commenceCalculateGravityLocal(); -#endif /// Entry point for the remote computation: for each bucket compute the /// force that its particles see due to the other particles NOT hosted diff --git a/State.h b/State.h index 261a6cee..78a592dd 100644 --- a/State.h +++ b/State.h @@ -89,23 +89,7 @@ class GenericList{ CkAssert(bucketStart >= 0); } - void getActiveBucketParameters(TreePiece *tp, - int bucket, - int &bucketStart, int &bucketSize){ - //std::map&lpref){ - // bucket is listed in this offload - GenericTreeNode *bucketNode = tp->bucketList[bucket]; - BucketActiveInfo *binfo = &(tp->bucketActiveInfo[bucket]); - - //bucketSize = bucketNode->lastParticle - bucketNode->firstParticle + 1; - //bucketStart = bucketNode->bucketArrayIndex; - bucketSize = tp->bucketActiveInfo[bucket].size; - bucketStart = tp->bucketActiveInfo[bucket].start; - CkAssert(bucketStart >= 0); - } - void push_back(int b, T &ilc, DoubleWalkState *state, TreePiece *tp); - }; @@ -168,6 +152,9 @@ class DoubleWalkState : public State { std::unordered_map nodeMap; std::unordered_map partMap; + // TODO do these need to be shut off (no restriction)? + // The PELists still need to collect the entire interaction list in memory + // before sending it to the GPU, so this isn't saving us from running OOM bool nodeOffloadReady(){ return nodeLists.totalNumInteractions >= nodeThreshold; } @@ -191,25 +178,6 @@ class DoubleWalkState : public State { partMap.reserve(100); #endif } - -#ifdef HAPI_INSTRUMENT_WRS - void nodeListConstructionTimeStart(){ - nodeListTime = CmiWallTimer(); - } - - double nodeListConstructionTimeStop(){ - return CmiWallTimer()-nodeListTime; - } - - void partListConstructionTimeStart(){ - partListTime = CmiWallTimer(); - } - - double partListConstructionTimeStop(){ - return CmiWallTimer()-partListTime; - } - -#endif }; #endif // INTERLIST_VER diff --git a/TreePiece.cpp b/TreePiece.cpp index 6d53b33d..19d14e03 100644 --- a/TreePiece.cpp +++ b/TreePiece.cpp @@ -44,6 +44,10 @@ #include "ckmulticast.h" #endif +#ifdef CUDA +#include "PEList.h" +#endif + using namespace std; using namespace SFC; using namespace TreeStuff; @@ -99,9 +103,6 @@ void TreePiece::setPeriodic(int nRepsPar, // Number of replicas in bPeriodic = bPeriodPar; bComove = bComovePar; dRhoFac = dRhoFacPar; - if(ewt == NULL) { - ewt = new EWT[nMaxEwhLoop]; - } } // Scale velocities (needed to convert to canonical momenta for @@ -1959,7 +1960,7 @@ void TreePiece::countActive(int activeRung, const CkCallback& cb) { } } } - numActiveParticles = nActive[0]; + myNumActiveParticles = nActive[0]; contribute(2*sizeof(int64_t), nActive, CkReduction::sum_long, cb); } @@ -3665,9 +3666,11 @@ int reEncodeOffset(int reqID, int offsetID) /** * Initialize all particles for gravity force calculation. * This includes zeroing out the acceleration and potential. + * @returns Number of active particles on this TreePiece. */ -void TreePiece::initBuckets() { +int TreePiece::initBuckets() { int ewaldCondition = (bEwald ? 0 : 1); + int nActive = 0; for (unsigned int j=0; jfirstParticle; i <= node->lastParticle; ++i) { if (myParticles[i].rung >= activeRung) { + nActive++; myParticles[i].treeAcceleration = 0; myParticles[i].potential = 0; myParticles[i].dtGrav = 0; @@ -3707,6 +3711,7 @@ void TreePiece::initBuckets() { #if COSMO_DEBUG > 1 || defined CHANGA_REFACTOR_WALKCHECK || defined CHANGA_REFACTOR_WALKCHECK_INTERLIST bucketcheckList.resize(numBuckets); #endif + return nActive; } void TreePiece::startNextBucket() { @@ -3809,7 +3814,9 @@ void TreePiece::finishBucket(int iBucket) { CkAssert(remaining >= 0); #ifdef COSMO_PRINT - CkPrintf("[%d] Is finished %d? finished=%d, %d still missing!\n",thisIndex,iBucket,req->finished, remaining); + CkPrintf("[%d %d] Is finished %d? finished=%d, %d still missing; remote %d local %d!\n", + CkMyPe(), thisIndex,iBucket,req->finished, remaining, + sRemoteGravityState->counterArrays[0][iBucket], sLocalGravityState->counterArrays[0][iBucket]); #endif // XXX finished means Ewald is done. @@ -3817,10 +3824,18 @@ void TreePiece::finishBucket(int iBucket) { sLocalGravityState->myNumParticlesPending -= 1; #ifdef COSMO_PRINT_BK - CkPrintf("[%d] Finished bucket %d, %d particles remaining\n",thisIndex,iBucket, sLocalGravityState->myNumParticlesPending); + CkPrintf("[%d %d] Finished bucket %d, %d buckets remaining\n",CkMyPe(), thisIndex,iBucket, sLocalGravityState->myNumParticlesPending); #endif if(sLocalGravityState->myNumParticlesPending == 0) { +#ifdef CUDA + if (!bUseCpu) { + for (int i = 0; i < numPEListProxies; i++) { + PEListProxies[i]->ckLocalBranch()->finishWalk(); + } + } +#endif + if(verbosity>1){ #if COSMO_STATS > 0 CkPrintf("[%d] TreePiece %d finished with bucket %d , openCriterions:%lld\n",CkMyPe(),thisIndex,iBucket,numOpenCriterionCalls); @@ -3829,70 +3844,37 @@ void TreePiece::finishBucket(int iBucket) { #endif } - if (bUseCpu) - { - continueWrapUp(); - } -#if defined CUDA - else { - // in cuda version, must wait till particle accels. - // are copied back from gpu; can markwalkdone only - // after this, otherwise there is a race condition - // between the start of the next iteration, wherein - // the treepiece registers itself afresh with the - // data manager. if during this time the particles - // haven't been copied (and updateParticles hasn't - // been called), the registeredTreePieces list will - // not be reset, so that the data manager gets - // confused. - dm->transferParticleVarsBack(); - bUseCpu = 1; - } -#endif + if (bUseCpu) continueWrapUp(); } } } #ifdef CUDA /// @brief Fill GPU buffer with particle data -/// @param bufLocalParts GPU buffer for particles -/// @param bufLocalMoments GPU buffer for Moments -/// @param pLocalMoments pointer to vector of Moments to be copied into the GPU -/// buffer -/// @param partIndex index into bufLocalParts at which to copy this -/// TreePieces particles -/// @param nParts total number of particles to be transfered to the -/// GPU (pass through) -/// @param node Root node of tree walk (pass through) -void TreePiece::fillGPUBuffer(intptr_t bufLocalParts, - intptr_t bufLocalMoments, - intptr_t pLocalMoments, int partIndex, int nParts, intptr_t node) +/// @param msg A struct containing info on where to write in the shared buffer +void TreePiece::fillGPUBuffer(fillGPUMsg *msg) { - CompactPartData *aLocalParts = (CompactPartData *)bufLocalParts; - CudaMultipoleMoments *aLocalMoments = (CudaMultipoleMoments *)pLocalMoments; + int partIndex = msg->partIndex; + int nParts = msg->nParts; + + CompactPartData *aLocalParts = dm->getBufLocalParts(); + CudaMultipoleMoments *aLocalMoments = dm->getLocalMoments(); getDMParticles(aLocalParts, partIndex); #ifdef GPU_LOCAL_TREE_WALK // set the bucketStart and bucketSize for each bucket Node - if (largePhase()) { - for (int j = 0; j < numBuckets; ++j) { - GenericTreeNode *bucketNode = bucketList[j]; - int id = bucketNode->nodeArrayIndex; - aLocalMoments[id].bucketStart = bucketNode->bucketArrayIndex; - aLocalMoments[id].bucketSize = bucketNode->lastParticle - - bucketNode->firstParticle + 1; - } - } else { - for (int j = 0; j < numBuckets; ++j) { - GenericTreeNode *bucketNode = bucketList[j]; - int id = bucketNode->nodeArrayIndex; - aLocalMoments[id].bucketStart = bucketActiveInfo[id].start; - aLocalMoments[id].bucketSize = bucketActiveInfo[id].size; - } + for (int j = 0; j < numBuckets; ++j) { + GenericTreeNode *bucketNode = bucketList[j]; + int id = bucketNode->nodeArrayIndex; + CkAssert(id >= 0); + aLocalMoments[id].bucketStart = bucketNode->bucketArrayIndex; + aLocalMoments[id].bucketSize = bucketNode->lastParticle + - bucketNode->firstParticle + 1; } // tell each particle which node it belongs to for (int j = 0; j < numBuckets; ++j) { GenericTreeNode *bucketNode = bucketList[j]; int id = bucketNode->nodeArrayIndex; + CkAssert(id >= 0); int start = aLocalMoments[id].bucketStart; int end = start + aLocalMoments[id].bucketSize; for (int k = start; k < end; k ++) { @@ -3900,31 +3882,35 @@ void TreePiece::fillGPUBuffer(intptr_t bufLocalParts, } } #endif - dm->transferLocalToGPU(nParts, (GenericTreeNode *)node); + delete msg; + bGPUBufferFilled = 1; + dm->transferLocalToGPU(nParts); } /// @brief update particle accelerations with GPU results void TreePiece::updateParticles(intptr_t data, int partIndex) { VariablePartData *deviceParticles = ((UpdateParticlesStruct *)data)->buf; - for(int j = 1; j <= myNumParticles; j++){ - if(isActive(j)){ -#ifndef CUDA_NO_ACC_UPDATES - myParticles[j].treeAcceleration.x += deviceParticles[partIndex].a.x; - myParticles[j].treeAcceleration.y += deviceParticles[partIndex].a.y; - myParticles[j].treeAcceleration.z += deviceParticles[partIndex].a.z; - myParticles[j].potential += deviceParticles[partIndex].potential; - myParticles[j].dtGrav = fmax(myParticles[j].dtGrav, - deviceParticles[partIndex].dtGrav); -#endif - if(!largePhase()) partIndex++; - } - if(largePhase()) partIndex++; + if(getNumActiveParticles() > 0) { + for(int j = 1; j <= myNumParticles; j++){ + if(isActive(j)){ + #ifndef CUDA_NO_ACC_UPDATES + myParticles[j].treeAcceleration.x += deviceParticles[partIndex].a.x; + myParticles[j].treeAcceleration.y += deviceParticles[partIndex].a.y; + myParticles[j].treeAcceleration.z += deviceParticles[partIndex].a.z; + myParticles[j].potential += deviceParticles[partIndex].potential; + myParticles[j].dtGrav = fmax(myParticles[j].dtGrav, + deviceParticles[partIndex].dtGrav); + #endif + } + partIndex++; } - - dm->updateParticlesFreeMemory((UpdateParticlesStruct *)data); - continueWrapUp(); } + dm->updateParticlesFreeMemory((UpdateParticlesStruct *)data); + if(getNumActiveParticles() > 0) + continueWrapUp(); // Only appropriate if a walk was + // actually started +} #endif void TreePiece::continueWrapUp(){ @@ -3966,28 +3952,39 @@ void TreePiece::doAllBuckets(){ thisProxy[thisIndex].nextBucket(msg); } -#ifdef GPU_LOCAL_TREE_WALK - else { - ListCompute *listcompute = (ListCompute *) sGravity; - DoubleWalkState *state = (DoubleWalkState *)sLocalGravityState; +} - listcompute->sendLocalTreeWalkTriggerToGpu(state, this, activeRung, 0, numBuckets); - // - // Set up the book keeping flags - bool useckloop = false; - for (int i = 0; i < numBuckets; i ++) { - sLocalGravityState->currentBucket = i; - GenericTreeNode *target = bucketList[i]; - if(target->rungs >= activeRung){ - doBookKeepingForTargetActive(i, i+1, -1, !useckloop, sLocalGravityState); - } else { - i += doBookKeepingForTargetInactive(-1, !useckloop, sLocalGravityState) - 1; - } - } - listcompute->resetCudaNodeState(state); - listcompute->resetCudaPartState(state); +/// @brief Call finishBucket for all buckets on this node +/// Used by local tree walk and Ewald GPU operations +/// @param fromEwald Flags whether this function was called after an Ewald calculation +void TreePiece::cudaFinishAllBuckets(int fromEwald){ + ListCompute *listcompute = (ListCompute *) sGravity; + DoubleWalkState *state = (DoubleWalkState *)sLocalGravityState; + + for (int i = 0; i < numBuckets; ++i) { + if (fromEwald) bucketReqs[i].finished = 1; + else state->counterArrays[0][i]--; + finishBucket(i); + } +} + +/// @brief Call finishBucket for a subset of buckets on this node +/// Used by PEList after an interaction list has been sent to the GPU +/// @param affectedBuckets Indices of buckets involved in the interaction +/// @param numBuckets Size of the affectedBuckets array +/// @param bRemote Flag to indicate whether to update remote gravity bookkeeping +void TreePiece::cudaFinishAffectedBuckets(int *affectedBuckets, int numBuckets, int bRemote) { + DoubleWalkState *state; + if (bRemote) state = (DoubleWalkState *)sRemoteGravityState; + else state = (DoubleWalkState *)sLocalGravityState; + + int bucket; + + for (int i = 0; i < numBuckets; ++i) { + bucket = affectedBuckets[i]; + state->counterArrays[0][bucket]--; + finishBucket(bucket); } -#endif } void TreePiece::nextBucket(dummyMsg *msg){ @@ -4149,14 +4146,8 @@ void TreePiece::nextBucket(dummyMsg *msg){ // these are local or remote interactions, which it learns from // sGravity. sGravity->init(NULL, activeRung, sLocal); - if(ds->nodeLists.totalNumInteractions > 0){ - lc->sendNodeInteractionsToGpu(ds, this); - lc->resetCudaNodeState(ds); - } - if(ds->particleLists.totalNumInteractions > 0){ - lc->sendPartInteractionsToGpu(ds, this); - lc->resetCudaPartState(ds); - } + thisProxy[thisIndex].flushInteractionsToGpu((intptr_t) lc, + (intptr_t) ds); } } #endif @@ -4196,54 +4187,40 @@ void TreePiece::calculateGravityLocal() { doAllBuckets(); } -void TreePiece::ewaldCPU(EwaldMsg *msg) { - bool useckloop = false; - int yield_num = _yieldPeriod; - - if (bUseCkLoopPar && otherIdlePesAvail()) { - useckloop = true; - // This value was chosen to be32*Nodesize so that we have enough buckets for - // all the PEs in the node and also giving some extra for load balance. - yield_num = 3 * CkMyNodeSize(); - calculateEwaldUsingCkLoop(yield_num); - } - - unsigned int i=0; - while (i < yield_num && ewaldCurrentBucket < numBuckets) { - if (!useckloop) { - BucketEwald(bucketList[ewaldCurrentBucket], nReplicas, fEwCut); +/// @brief Start the ewald calculation on this TreePiece +void TreePiece::calculateEwald(dummyMsg *msg) { + if (bUseCpu) + { + bool useckloop = false; + int yield_num = _yieldPeriod; + + if (bUseCkLoopPar && otherIdlePesAvail()) { + useckloop = true; + // This value was chosen to be32*Nodesize so that we have enough buckets for + // all the PEs in the node and also giving some extra for load balance. + yield_num = 3 * CkMyNodeSize(); + calculateEwaldUsingCkLoop(yield_num); } - bucketReqs[ewaldCurrentBucket].finished = 1; - finishBucket(ewaldCurrentBucket); + unsigned int i=0; + while (i < yield_num && ewaldCurrentBucket < numBuckets) { + if (!useckloop) { + BucketEwald(bucketList[ewaldCurrentBucket], nReplicas, fEwCut); + } - ewaldCurrentBucket++; - i++; - } + bucketReqs[ewaldCurrentBucket].finished = 1; + finishBucket(ewaldCurrentBucket); - if (ewaldCurrentBucketfromInit){ - thisProxy[thisIndex].EwaldGPU(); + if (ewaldCurrentBucketnodeLists.totalNumInteractions > 0){ - lc->sendNodeInteractionsToGpu(ds, this); - lc->resetCudaNodeState(ds); - } - if(ds->particleLists.totalNumInteractions > 0){ - lc->sendPartInteractionsToGpu(ds, this); - lc->resetCudaPartState(ds); - } + thisProxy[thisIndex].flushInteractionsToGpu((intptr_t) lc, + (intptr_t) ds); } } @@ -4595,18 +4566,16 @@ void TreePiece::calculateGravityRemote(ComputeChunkMsg *msg) { if(lc && ds){ - if(ds->nodeLists.totalNumInteractions > 0){ - lc->sendNodeInteractionsToGpu(ds, this); - lc->resetCudaNodeState(ds); - } - if(ds->particleLists.totalNumInteractions > 0){ - lc->sendPartInteractionsToGpu(ds, this); - lc->resetCudaPartState(ds); - } + thisProxy[thisIndex].flushInteractionsToGpu((intptr_t) lc, + (intptr_t) ds); } } #endif + if (particleInterRemote == NULL) { + CkPrintf("ERROR [%d] TP %d particleInterRemote NULL at chunk %d (calculateGravityRemote)\n", CkMyPe(), thisIndex, msg->chunkNum); + CkAbort("particleInterRemote NULL"); + } cacheGravPart[CkMyPe()].finishedChunk(msg->chunkNum, particleInterRemote[msg->chunkNum]); #ifdef CHECK_WALK_COMPLETIONS CkPrintf("[%d] finishedChunk TreePiece::calculateGravityRemote\n", thisIndex); @@ -4636,7 +4605,7 @@ int TreePiece::doBookKeepingForTargetActive(int curbucket, int end, // initialize in startiteration: +1 // for request sent to gpu: +1 // initial walk completed: -1 (this is where we are in the code currently) - // for each request completed: -1 (this happens in cudaCallback) + // for each request completed: -1 (this happens in cudaFinishAllBuckets) gravityState->counterArrays[0][j]--; #if COSMO_PRINT_BK > 1 CkPrintf("[%d] bucket %d numAddReq: %d,%d\n", thisIndex, j, sRemoteGravityState->counterArrays[0][j], sLocalGravityState->counterArrays[0][j]); @@ -5045,6 +5014,9 @@ void TreePiece::startGravity(int am, // the active mask for multistepping theta = myTheta; thetaMono = theta*theta*theta*theta; +#ifdef CUDA + bGPUBufferFilled = 0; +#endif starttime = CmiWallTimer(); int oldNumChunks = numChunks; @@ -5055,6 +5027,25 @@ void TreePiece::startGravity(int am, // the active mask for multistepping if (numChunks == 0 && myNumParticles == 0) numChunks = 1; int dummy; + // Allocate particleInterRemote/nodeInterRemote before PEList::finishWalk() + // is invoked (e.g. in the myNumParticles==0 path via finishedTPWork). + // They must be valid if any callback arrives. PUP sets them to NULL on + // unpack; allocation must happen before any cache/Compute path can touch them. + if (oldNumChunks != numChunks) { + delete[] nodeInterRemote; + delete[] particleInterRemote; + nodeInterRemote = new u_int64_t[numChunks]; + particleInterRemote = new u_int64_t[numChunks]; + } + if (nodeInterRemote == NULL) + nodeInterRemote = new u_int64_t[numChunks]; + if (particleInterRemote == NULL) + particleInterRemote = new u_int64_t[numChunks]; + for (int i = 0; i < numChunks; ++i) { + nodeInterRemote[i] = 0; + particleInterRemote[i] = 0; + } + cacheNode.ckLocalBranch()->cacheSync(numChunks, idxMax, localIndex); cacheGravPart.ckLocalBranch()->cacheSync(numChunks, idxMax, dummy); @@ -5069,31 +5060,16 @@ void TreePiece::startGravity(int am, // the active mask for multistepping CkCallback cbf = CkCallback(CkIndex_TreePiece::finishWalk(), pieces); gravityProxy[thisIndex].ckLocal()->contribute(cbf); bBucketsInited = true; + myNumActiveParticles = 0; return; } - - // allocate and zero out statistics counters - if (oldNumChunks != numChunks ) { - delete[] nodeInterRemote; - delete[] particleInterRemote; - nodeInterRemote = new u_int64_t[numChunks]; - particleInterRemote = new u_int64_t[numChunks]; - } - if(nodeInterRemote == NULL) - nodeInterRemote = new u_int64_t[numChunks]; - if(particleInterRemote == NULL) - particleInterRemote = new u_int64_t[numChunks]; #if COSMO_STATS > 0 nodesOpenedLocal = 0; nodesOpenedRemote = 0; numOpenCriterionCalls=0; #endif nodeInterLocal = 0; - for (int i=0; i1) @@ -5112,8 +5088,35 @@ void TreePiece::startGravity(int am, // the active mask for multistepping traceUserBracketEvent(START_REG, starttime, CmiWallTimer()); starttime = CmiWallTimer(); - initBuckets(); + myNumActiveParticles = initBuckets(); traceUserBracketEvent(START_IB, starttime, CmiWallTimer()); + if(getNumActiveParticles() == 0) { + // No active particles on this TreePiece; Short-circuit the tree + // walk. This is similar to the code above for no particles on + // the TreePiece. + for (int i=0; i< numChunks; ++i) { + cacheGravPart.ckLocalBranch()->finishedChunk(i, 0); + } + nodeLBMgrProxy.ckLocalBranch()->finishedTPWork(); +#ifdef CUDA + if (!bUseCpu) { + // This trees particles and nodes still need to be loaded + // onto the GPU + dm->serializeLocalTree(); + // Let DM know we have "completed" the prefetch as well + dm->donePrefetch(0); // 0 is the chunk number; assuming + // only one chunk + // Every PE must call TransferParticleVarsBack, even if no particle data + for (int i = 0; i < numPEListProxies; i++) { + PEListProxies[i]->ckLocalBranch()->finishWalk(); + } + } +#endif + CkCallback cbf = CkCallback(CkIndex_TreePiece::finishWalk(), pieces); + gravityProxy[thisIndex].ckLocal()->contribute(cbf); + bBucketsInited = true; + return; + } starttime = CmiWallTimer(); prefetchReq.reset(); @@ -5193,7 +5196,6 @@ void TreePiece::startGravity(int am, // the active mask for multistepping #ifdef CUDA numActiveBuckets = 0; - calculateNumActiveParticles(); if (!bUseCpu) { @@ -5220,26 +5222,7 @@ void TreePiece::startGravity(int am, // the active mask for multistepping // ditto lstate->nodes = NULL; // allocate space for local particles - // if this is a small phase; we do not transfer - // all particles owned by this processor in small - // phases and so refer to particles inside an - // auxiliary array shipped with computation requests. - // this auxiliary array is 'particles', below: - if(largePhase()){ - lstate->particles = NULL; - } - else{ - // allocate an amount of space that - // depends on the rung - lstate->particles = new CkVec(AVG_SOURCE_PARTICLES_PER_ACTIVE*myNumActiveParticles); - lstate->particles->length() = 0; - // need to allocate memory for data structure that stores bucket - // active info (where in the gpu's target particle memory this - // bucket starts, and its size; strictly speaking, we don't need - // the size attribute.) - // XXX - no need to allocate/delete every iteration - bucketActiveInfo = new BucketActiveInfo[numBuckets]; - } + lstate->particles = NULL; } { DoubleWalkState *state = (DoubleWalkState *)sInterListStateRemoteResume; @@ -5312,7 +5295,14 @@ void TreePiece::startGravity(int am, // the active mask for multistepping #endif traceUserBracketEvent(START_PW, starttime, CmiWallTimer()); - if (bEwald) thisProxy[thisIndex].EwaldInit(); + if (bEwald) { + dummyMsg *msg = new (8*sizeof(int)) dummyMsg; + // Make priority lower than gravity or smooth. + *((int *)CkPriorityPtr(msg)) = 3*numTreePieces * numChunks + thisIndex + 1; + CkSetQueueing(msg,CK_QUEUEING_IFIFO); + thisProxy[thisIndex].calculateEwald(msg); + } + #if defined CUDA // ask datamanager to serialize local trees // prefetch can occur concurrently with this, @@ -5320,13 +5310,11 @@ void TreePiece::startGravity(int am, // the active mask for multistepping // afterwards. if (!bUseCpu) { dm->serializeLocalTree(); - } else { - thisProxy[thisIndex].commenceCalculateGravityLocal(0, 0, 0, 0, 0, 0, 0, 0); - } -#else - thisProxy[thisIndex].commenceCalculateGravityLocal(); + } else #endif - + { + thisProxy[thisIndex].commenceCalculateGravityLocal(); + } #ifdef CHANGA_PRINT_MEMUSAGE int piecesPerPe = numTreePieces/CmiNumPes(); @@ -5375,27 +5363,7 @@ void TreePiece::initiatePrefetch(int chunk){ } /// @brief Entry method wrapper for calculateGravityLocal -/// If using the GPU, this TreePiece is assigned a cudaStream and given -/// handles to device memory -#ifdef CUDA -void TreePiece::commenceCalculateGravityLocal(intptr_t d_localMoments, - intptr_t d_localParts, - intptr_t d_localVars, - intptr_t streams, int numStreams, - size_t sMoments, size_t sCompactParts, size_t sVarParts) { - if (!bUseCpu) { - this->d_localMoments = (CudaMultipoleMoments *)d_localMoments; - this->d_localParts = (CompactPartData *)d_localParts; - this->d_localVars = (VariablePartData *)d_localVars; - this->stream = ((cudaStream_t *)streams)[thisIndex % numStreams]; - this->sMoments = sMoments; - this->sCompactParts = sCompactParts; - this->sVarParts = sVarParts; - } - -#else void TreePiece::commenceCalculateGravityLocal(){ -#endif #if INTERLIST_VER > 0 // must set placedRoots to false before starting local comp. DoubleWalkState *lstate = (DoubleWalkState *)sLocalGravityState; @@ -5416,9 +5384,11 @@ void TreePiece::startRemoteChunk() { // dm counts until all treepieces have acknowledged prefetch completion // it then flattens the tree on the processor, sends it to the device // and sends messages to each of the registered treepieces to continueStartRemoteChunk() + // note that continueStartRemoteChunk needs both the local and remote data + // serialization to complete first dm->donePrefetch(sPrefetchState->currentBucket); } else { - continueStartRemoteChunk(sPrefetchState->currentBucket, 0, 0); + continueStartRemoteChunk(sPrefetchState->currentBucket); } #else continueStartRemoteChunk(sPrefetchState->currentBucket); @@ -5428,15 +5398,7 @@ void TreePiece::startRemoteChunk() { /// @brief Main work of StartRemoteChunk() /// Schedule a TreePiece::calculateGravityRemote() then start /// prefetching for the next chunk. -#ifdef CUDA -void TreePiece::continueStartRemoteChunk(int chunk, intptr_t d_remoteMoments, intptr_t d_remoteParts){ - if (!bUseCpu) { - this->d_remoteMoments = (CudaMultipoleMoments *)d_remoteMoments; - this->d_remoteParts = (CompactPartData *)d_remoteParts; - } -#else void TreePiece::continueStartRemoteChunk(int chunk){ -#endif // FIXME - can value of chunk be different from current Prefetch? ComputeChunkMsg *msg = new (8*sizeof(int)) ComputeChunkMsg(sPrefetchState->currentBucket); *(int*)CkPriorityPtr(msg) = numTreePieces * numChunks + thisIndex + 1; @@ -5475,13 +5437,13 @@ void TreePiece::continueStartRemoteChunk(int chunk){ /// @param activeRung Rung to use. void TreePiece::setTreePieceLoad(int activeRung) { double dLoadExp; - nPrevActiveParts = numActiveParticles; + nPrevActiveParts = myNumActiveParticles; if (havePhaseData(activeRung)) { dLoadExp = savedPhaseLoad[activeRung]; } else if (havePhaseData(0)) { float ratio = 1.0; if(myNumParticles != 0){ - ratio = numActiveParticles/(float)myNumParticles; + ratio = myNumActiveParticles/(float)myNumParticles; } dLoadExp = ratio * savedPhaseLoad[0]; @@ -5519,23 +5481,23 @@ void TreePiece::startlb(const CkCallback &cb, int activeRung, bool bDoLB){ // We need to recount the number of active particles since DD has // moved particles around if(activeRung == 0){ // Everybody is active; no need to count - numActiveParticles = myNumParticles; + myNumActiveParticles = myNumParticles; } else if(activeRung == PHASE_FEEDBACK) { // Also no need to recount - numActiveParticles = myNumSPH + myNumStar; + myNumActiveParticles = myNumSPH + myNumStar; } else{ - numActiveParticles = 0; + myNumActiveParticles = 0; for(unsigned int i = 1; i <= myNumParticles; ++i) { if(myParticles[i].rung >= activeRung) { - numActiveParticles++; + myNumActiveParticles++; } } } LDObjHandle myHandle = myRec->getLdHandle(); - TaggedVector3D tv(savedCentroid, myHandle, numActiveParticles, myNumParticles, + TaggedVector3D tv(savedCentroid, myHandle, getNumActiveParticles(), myNumParticles, iActiveRungLB, iPrevRungLB); tv.tp = thisIndex; tv.tag = thisIndex; @@ -5815,6 +5777,7 @@ void TreePiece::pup(PUP::er& p) { p | myNumParticles; p | nTotalSPH; p | myNumSPH; + p | myNumActiveParticles; p | nTotalDark; p | nTotalStar; p | myNumStar; @@ -5879,10 +5842,6 @@ void TreePiece::pup(PUP::er& p) { p | bPeriodic; p | bComove; p | dRhoFac; - p | nMaxEwhLoop; - if (p.isUnpacking() && bEwald) { - ewt = new EWT[nMaxEwhLoop]; - } p | numBuckets; #if INTERLIST_VER > 0 @@ -6373,11 +6332,6 @@ void TreePiece::freeWalkObjects(){ delete state->particles; delete rstate->nodes; delete rstate->particles; - if(!largePhase()){ - DoubleWalkState *lstate = (DoubleWalkState *) sLocalGravityState; - delete lstate->particles; - delete [] bucketActiveInfo; - } } #endif // remote-no-resume state @@ -6434,9 +6388,10 @@ void TreePiece::markWalkDone() { // there may be outstanding requests by other pieces. We need to // wait for all walks to complete before freeing data structures. #ifdef CHECK_WALK_COMPLETIONS - CkPrintf("[%d] inside markWalkDone, completedActiveWalks: %d, activeWalks: %d, contrib finishWalk\n", thisIndex, completedActiveWalks, activeWalks.size()); + CkPrintf("[%d] inside markWalkDone, completedActiveWalks: %d, activeWalks: %d, contrib finishWalk\n", thisIndex, completedActiveWalks, (int) activeWalks.size()); #endif nodeLBMgrProxy.ckLocalBranch()->finishedTPWork(); + CkAssert(myNumParticles > 0 && getNumActiveParticles() > 0); CkCallback cb = CkCallback(CkIndex_TreePiece::finishWalk(), pieces); gravityProxy[thisIndex].ckLocal()->contribute(cb); } @@ -6457,13 +6412,14 @@ void TreePiece::finishWalk() #endif #ifdef HAPI_TRACE - CkPrintf("[%d] (%d) HAPI_TRACE localnode: %lld\n", thisIndex, activeRung, localNodeInteractions); - CkPrintf("[%d] (%d) HAPI_TRACE remotenode: %lld\n", thisIndex, activeRung, remoteNodeInteractions); - CkPrintf("[%d] (%d) HAPI_TRACE remoteresumenode: %lld\n", thisIndex, activeRung, remoteResumeNodeInteractions); - CkPrintf("[%d] (%d) HAPI_TRACE localpart: %lld\n", thisIndex, activeRung, localPartInteractions); - CkPrintf("[%d] (%d) HAPI_TRACE remotepart: %lld\n", thisIndex, activeRung, remotePartInteractions); - CkPrintf("[%d] (%d) HAPI_TRACE remoteresumepart: %lld\n", thisIndex, activeRung, remoteResumePartInteractions); - + if(verbosity > 1) { + CkPrintf("[%d] (%d) HAPI_TRACE localnode: %lld\n", thisIndex, activeRung, localNodeInteractions); + CkPrintf("[%d] (%d) HAPI_TRACE remotenode: %lld\n", thisIndex, activeRung, remoteNodeInteractions); + CkPrintf("[%d] (%d) HAPI_TRACE remoteresumenode: %lld\n", thisIndex, activeRung, remoteResumeNodeInteractions); + CkPrintf("[%d] (%d) HAPI_TRACE localpart: %lld\n", thisIndex, activeRung, localPartInteractions); + CkPrintf("[%d] (%d) HAPI_TRACE remotepart: %lld\n", thisIndex, activeRung, remotePartInteractions); + CkPrintf("[%d] (%d) HAPI_TRACE remoteresumepart: %lld\n", thisIndex, activeRung, remoteResumePartInteractions); + } #endif gravityProxy[thisIndex].ckLocal()->contribute(cbGravity); diff --git a/cuda.mk.in b/cuda.mk.in index b420b3a5..eef5f949 100644 --- a/cuda.mk.in +++ b/cuda.mk.in @@ -1,21 +1,16 @@ # optional CUDA flags: # memory: # -DPINNED_HOST_MEMORY -# -DHAPI_MEMPOOL # # verbosity, debugging: # -DCUDA_DM_PRINT_TREES # -DCUDA_PRINT_TRANSFERRED_INTERACTIONS # -DCUDA_PRINT_TRANSFER_BACK_PARTICLES -# -DCUDA_NOTIFY_DATA_TRANSFER_DONE -# -DCUDA_VERBOSE_KERNEL_ENQUEUE +# -DCUDA_VERBOSE_OPS # -DCUDA_NO_KERNELS # -DCUDA_NO_ACC_UPDATES # -DCUDA_UNIT_TEST # -# print errors returned by CUDA calls: -# -DCUDA_PRINT_ERRORS -# # for performance monitoring via projections/stats # -DCUDA_STATS # -DHAPI_TRACE @@ -26,15 +21,13 @@ CUDA_DEBUG = $(DEBUG) CUDA_VERBOSE = $(VERBOSE) cuda_defines += -DCUDA -DSPCUDA -cuda_defines += -DHAPI_MEMPOOL -DPINNED_HOST_MEMORY -cuda_defines += -DCUDA_2D_TB_KERNEL -#cuda_defines += -DCUDA_2D_FLAT +cuda_defines += -DPINNED_HOST_MEMORY # Do the local treewalk on the GPU cuda_defines += @FLAG_GPU_LOCAL_TREE_WALK@ ifeq (1,$(CUDA_DEBUG)) - cuda_defines += -DCUDA_VERBOSE_KERNEL_ENQUEUE + cuda_defines += -DCUDA_VERBOSE_OPS endif ifeq (yes,@PROJECTIONS@) diff --git a/cuda_typedef.h b/cuda_typedef.h index c69b65a5..74654d07 100644 --- a/cuda_typedef.h +++ b/cuda_typedef.h @@ -30,14 +30,15 @@ typedef cosmoType hosttype; #define CUDA_XFER_LOCAL 9910 #define CUDA_XFER_REMOTE 9911 -#define CUDA_GRAV_LOCAL 9912 -#define CUDA_GRAV_REMOTE 9913 -#define CUDA_REMOTE_RESUME 9914 -#define CUDA_PART_GRAV_LOCAL 9915 -#define CUDA_PART_GRAV_LOCAL_SMALL 9916 -#define CUDA_PART_GRAV_REMOTE 9917 -#define CUDA_XFER_BACK 9918 -#define CUDA_EWALD 9919 +#define CUDA_GRAV_TREE_LOCAL 9912 +#define CUDA_GRAV_NODELIST_LOCAL 9913 +#define CUDA_GRAV_PARTLIST_LOCAL 9914 +#define CUDA_GRAV_NODELIST_REMOTE 9915 +#define CUDA_GRAV_PARTLIST_REMOTE 9916 +#define CUDA_GRAV_NODELIST_REMOTE_RESUME 9917 +#define CUDA_GRAV_PARTLIST_REMOTE_RESUME 9918 +#define CUDA_XFER_BACK 9919 +#define CUDA_EWALD 9920 #define SER_LOCAL_WALK 9922 #define SER_LOCAL_GATHER 9923 @@ -46,8 +47,6 @@ typedef cosmoType hosttype; #endif -// TODO: Fix small phase code -#define TP_LARGE_PHASE_THRESHOLD_DEFAULT 0.0 #define AVG_SOURCE_PARTICLES_PER_ACTIVE 10 /** @brief 3D vector of cudatype. diff --git a/memlog.h b/memlog.h new file mode 100644 index 00000000..8435f8e1 --- /dev/null +++ b/memlog.h @@ -0,0 +1,51 @@ +#ifndef MEMLOG_HINCLUDED +#define MEMLOG_HINCLUDED + +// Enum to represent memory operation types +enum MemLogOpType { + MEMLOG_ALLOC, + MEMLOG_FREE, + MEMLOG_ALLOC_FAIL, + MEMLOG_FREE_FAIL, + MEMLOG_FREE_SKIP // For attempting to free NULL +}; + +/** @brief Holds information about a single CUDA memory event */ +class MemLogEvent { +public: + int nodeId; // ID of the node where the event occurred (CkMyNode()) + MemLogOpType opType; // Type of memory operation + size_t size; // Size of memory block (for alloc/alloc_fail) + uintptr_t address; // Memory address (device pointer) + double timestamp; // Timestamp of the event (e.g., CkWallTimer()) + std::string location; // Code location ("file:line") + std::string pointerId; // User-provided tag for the operation (usually variable name) + std::string functionTag; // Tag indicating the calling function/context + + // Default constructor + MemLogEvent() : nodeId(-1), opType(MEMLOG_ALLOC), size(0), address(0), timestamp(0.0) {} + + // Constructor + MemLogEvent(int _nodeId, MemLogOpType _opType, size_t _size, uintptr_t _address, double _timestamp, const char* _file, int _line, const char* _pointerId, const char* _functionTag) + : nodeId(_nodeId), opType(_opType), size(_size), address(_address), timestamp(_timestamp), pointerId(_pointerId ? _pointerId : ""), functionTag(_functionTag ? _functionTag : "") { + // Combine file and line into location string + location = std::string(_file ? _file : "unknown") + ":" + std::to_string(_line); + } +}; + +/** @brief Manages buffering and flushing of memory log events */ +class MemLog { +public: + std::vector meTab; // Buffer for memory events + std::string fileName; // Output log file name + + MemLog() : fileName(".memlog") {} // Default filename + + // Method to flush buffered events to the file (implementation in .cpp) + void flush(); + + // Method to write metadata/header to the log file (implementation in .cpp) + static void logMetaData(std::ofstream &ofsLog); // Static as it doesn't depend on instance data +}; + +#endif diff --git a/parameters.h b/parameters.h index 955368cc..2d469a52 100644 --- a/parameters.h +++ b/parameters.h @@ -35,6 +35,11 @@ typedef struct parameters { int iMaxRung; #ifdef CUDA int nGpuMinParts; + int bGpuMemLogger; + int bCpuMemLogger; + double dHostPoolTargetCapacityGB; + int nHostPoolMinCapacityPerBucketMB; + int bHostPoolDebug; #endif int bCannonical; int bKDK; @@ -189,6 +194,11 @@ inline void operator|(PUP::er &p, Parameters ¶m) { p|param.nTruncateRung; #ifdef CUDA p|param.nGpuMinParts; + p|param.bGpuMemLogger; + p|param.bCpuMemLogger; + p|param.dHostPoolTargetCapacityGB; + p|param.nHostPoolMinCapacityPerBucketMB; + p|param.bHostPoolDebug; #endif p|param.iMaxRung; p|param.bCannonical;