#define TILE_SIZE 32 /** * Compute a value based on pair interactions. */ __attribute__((reqd_work_group_size(64, 1, 1))) __kernel void computeN2Value(__global float4* posq, __local float4* local_posq, __global unsigned int* exclusions, __global unsigned int* exclusionIndices, __global float* global_value, __local float* local_value, __local float* tempBuffer, __global unsigned int* tiles, #ifdef USE_CUTOFF __global unsigned int* interactionFlags, __global unsigned int* interactionCount #else unsigned int numTiles #endif PARAMETER_ARGUMENTS) { #ifdef USE_CUTOFF unsigned int numTiles = interactionCount[0]; #endif unsigned int pos = get_group_id(0)*numTiles/get_num_groups(0); unsigned int end = (get_group_id(0)+1)*numTiles/get_num_groups(0); float energy = 0.0f; unsigned int lasty = 0xFFFFFFFF; while (pos < end) { // Extract the coordinates of this tile unsigned int x = tiles[pos]; unsigned int y = ((x >> 2) & 0x7fff)*TILE_SIZE; bool hasExclusions = (x & 0x1); x = (x>>17)*TILE_SIZE; unsigned int baseLocalAtom = (get_local_id(0) < TILE_SIZE ? 0 : TILE_SIZE/2); unsigned int tgx = get_local_id(0) & (TILE_SIZE-1); unsigned int valueBufferOffset = (tgx < TILE_SIZE/2 ? 0 : TILE_SIZE); unsigned int atom1 = x + tgx; float value = 0.0f; float4 posq1 = posq[atom1]; LOAD_ATOM1_PARAMETERS if (x == y) { // This tile is on the diagonal. local_posq[get_local_id(0)] = posq1; LOAD_LOCAL_PARAMETERS_FROM_1 barrier(CLK_LOCAL_MEM_FENCE); unsigned int xi = x/TILE_SIZE; unsigned int tile = xi+xi*PADDED_NUM_ATOMS/TILE_SIZE-xi*(xi+1)/2; #ifdef USE_EXCLUSIONS unsigned int excl = exclusions[exclusionIndices[tile]+tgx] >> baseLocalAtom; #endif for (unsigned int j = 0; j < TILE_SIZE/2; j++) { #ifdef USE_EXCLUSIONS bool isExcluded = !(excl & 0x1); #endif int atom2 = baseLocalAtom+j; float4 posq2 = local_posq[atom2]; float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f); #ifdef USE_PERIODIC delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X; delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y; delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z; #endif float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z; #ifdef USE_CUTOFF if (r2 < CUTOFF_SQUARED) { #endif float r = SQRT(r2); LOAD_ATOM2_PARAMETERS atom2 = y+baseLocalAtom+j; float tempValue1 = 0.0f; float tempValue2 = 0.0f; #ifdef USE_EXCLUSIONS if (!isExcluded && atom1 < NUM_ATOMS && atom2 < NUM_ATOMS && atom1 != atom2) { #else if (atom1 < NUM_ATOMS && atom2 < NUM_ATOMS && atom1 != atom2) { #endif COMPUTE_VALUE } value += tempValue1; #ifdef USE_CUTOFF } #endif #ifdef USE_EXCLUSIONS excl >>= 1; #endif } // Sum the values and write results. if (get_local_id(0) >= TILE_SIZE) tempBuffer[get_local_id(0)] = value; barrier(CLK_LOCAL_MEM_FENCE); if (get_local_id(0) < TILE_SIZE) { #ifdef USE_OUTPUT_BUFFER_PER_BLOCK unsigned int offset = x + tgx + (x/TILE_SIZE)*PADDED_NUM_ATOMS; #else unsigned int offset = x + tgx + get_group_id(0)*PADDED_NUM_ATOMS; #endif global_value[offset] += value+tempBuffer[get_local_id(0)+TILE_SIZE]; } } else { // This is an off-diagonal tile. if (lasty != y && get_local_id(0) < TILE_SIZE) { unsigned int j = y + tgx; local_posq[get_local_id(0)] = posq[j]; LOAD_LOCAL_PARAMETERS_FROM_GLOBAL } local_value[get_local_id(0)] = 0.0f; barrier(CLK_LOCAL_MEM_FENCE); // Compute the full set of interactions in this tile. unsigned int xi = x/TILE_SIZE; unsigned int yi = y/TILE_SIZE; unsigned int tile = xi+yi*PADDED_NUM_ATOMS/TILE_SIZE-yi*(yi+1)/2; #ifdef USE_EXCLUSIONS unsigned int excl = (hasExclusions ? exclusions[exclusionIndices[tile]+tgx] : 0xFFFFFFFF); excl = (excl >> tgx) | (excl << (TILE_SIZE - tgx)); excl >>= baseLocalAtom; #endif unsigned int tj = tgx%(TILE_SIZE/2); for (unsigned int j = 0; j < TILE_SIZE/2; j++) { #ifdef USE_EXCLUSIONS bool isExcluded = !(excl & 0x1); #endif int atom2 = baseLocalAtom+tj; float4 posq2 = local_posq[atom2]; float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f); #ifdef USE_PERIODIC delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X; delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y; delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z; #endif float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z; #ifdef USE_CUTOFF if (r2 < CUTOFF_SQUARED) { #endif float r = SQRT(r2); LOAD_ATOM2_PARAMETERS atom2 = y+baseLocalAtom+tj; float tempValue1 = 0.0f; float tempValue2 = 0.0f; #ifdef USE_EXCLUSIONS if (!isExcluded && atom1 < NUM_ATOMS && atom2 < NUM_ATOMS) { #else if (atom1 < NUM_ATOMS && atom2 < NUM_ATOMS) { #endif COMPUTE_VALUE } value += tempValue1; local_value[baseLocalAtom+tj+valueBufferOffset] += tempValue2; #ifdef USE_CUTOFF } #endif barrier(CLK_LOCAL_MEM_FENCE); #ifdef USE_EXCLUSIONS excl >>= 1; #endif tj = (tj+1)%(TILE_SIZE/2); } // Sum the values and write results. if (get_local_id(0) >= TILE_SIZE) tempBuffer[get_local_id(0)] = value; barrier(CLK_LOCAL_MEM_FENCE); if (get_local_id(0) < TILE_SIZE) { #ifdef USE_OUTPUT_BUFFER_PER_BLOCK unsigned int offset1 = x + tgx + (y/TILE_SIZE)*PADDED_NUM_ATOMS; unsigned int offset2 = y + tgx + (x/TILE_SIZE)*PADDED_NUM_ATOMS; #else unsigned int offset1 = x + tgx + get_group_id(0)*PADDED_NUM_ATOMS; unsigned int offset2 = y + tgx + get_group_id(0)*PADDED_NUM_ATOMS; #endif global_value[offset1] += value+tempBuffer[get_local_id(0)+TILE_SIZE]; global_value[offset2] += local_value[get_local_id(0)]+local_value[get_local_id(0)+TILE_SIZE]; } lasty = y; } pos++; } }