NAMD: CudaNonbondedTables.C Source File

00001 #include "charm++.h"
00002 #include "NamdTypes.h"
00003 #include "ComputeNonbondedUtil.h"
00004 #include "LJTable.h"
00005 #include "CudaUtils.h"
00006 #include "CudaNonbondedTables.h"
00007 
00008 #ifdef NAMD_CUDA
00009 
00010 CudaNonbondedTables::CudaNonbondedTables(const int deviceID) : deviceID(deviceID) {
00011 
00012   vdwCoefTable = NULL;
00013   vdwCoefTableWidth = 0;
00014   vdwCoefTableTex = 0;  
00015 
00016   forceTableTex = 0;
00017   energyTableTex = 0;
00018 
00019   exclusionTable = NULL;
00020   r2_table = NULL;
00021   exclusionTableTex = 0;
00022   r2_table_tex = 0;
00023 
00024   modifiedExclusionForceTableTex = 0;
00025   modifiedExclusionEnergyTableTex = 0;
00026 
00027   cudaCheck(cudaSetDevice(deviceID));
00028   buildForceAndEnergyTables(4096);
00029   buildVdwCoefTable();
00030 }
00031 
00032 CudaNonbondedTables::~CudaNonbondedTables() {
00033   cudaCheck(cudaSetDevice(deviceID));
00034   if (vdwCoefTable != NULL) deallocate_device<float2>(&vdwCoefTable);
00035   if (exclusionTable != NULL) deallocate_device<float4>(&exclusionTable);
00036   if (r2_table != NULL) deallocate_device<float>(&r2_table);
00037 
00038   cudaCheck(cudaFreeArray(forceArray));
00039   cudaCheck(cudaFreeArray(energyArray));
00040 
00041   cudaCheck(cudaDestroyTextureObject(vdwCoefTableTex));
00042   cudaCheck(cudaDestroyTextureObject(forceTableTex));
00043   cudaCheck(cudaDestroyTextureObject(energyTableTex));
00044 
00045   cudaCheck(cudaDestroyTextureObject(exclusionTableTex));
00046   cudaCheck(cudaDestroyTextureObject(r2_table_tex));
00047 }
00048 
00049 template <typename T>
00050 void bindTextureObject(int size, T* h_table, T*& d_table, cudaTextureObject_t& tex, bool update=false) {
00051   // Copy to device
00052   if ( ! update) {
00053     allocate_device<T>(&d_table, size);
00054   }
00055   else {
00056     cudaCheck(cudaDestroyTextureObject(tex));
00057   }
00058   copy_HtoD_sync<T>(h_table, d_table, size);
00059 
00060   // Create texture object
00061   cudaResourceDesc resDesc;
00062   memset(&resDesc, 0, sizeof(resDesc));
00063   resDesc.resType = cudaResourceTypeLinear;
00064   resDesc.res.linear.devPtr = d_table;
00065   resDesc.res.linear.desc.f = cudaChannelFormatKindFloat;
00066   resDesc.res.linear.desc.x = sizeof(float)*8; // bits per channel
00067   if (sizeof(T) >= sizeof(float)*2) resDesc.res.linear.desc.y = sizeof(float)*8; // bits per channel
00068   if (sizeof(T) >= sizeof(float)*3) resDesc.res.linear.desc.z = sizeof(float)*8; // bits per channel
00069   if (sizeof(T) >= sizeof(float)*4) resDesc.res.linear.desc.w = sizeof(float)*8; // bits per channel
00070   resDesc.res.linear.sizeInBytes = size*sizeof(T);
00071 
00072   cudaTextureDesc texDesc;
00073   memset(&texDesc, 0, sizeof(texDesc));
00074   texDesc.readMode = cudaReadModeElementType;
00075 
00076   cudaCheck(cudaCreateTextureObject(&tex, &resDesc, &texDesc, NULL));
00077 }
00078 
00079 //
00080 // Builds the VdW Lennard-Jones coefficient table. Swiped from ComputeNonbondedCUDA.C
00081 // NOTE: Can only be called once
00082 //
00083 void CudaNonbondedTables::buildVdwCoefTable(bool update/*=false*/) {
00084 
00085   const LJTable* const ljTable = ComputeNonbondedUtil:: ljTable;
00086   const int dim = ljTable->get_table_dim();
00087 
00088   // round dim up to odd multiple of 16
00089   int tsize = (((dim+16+31)/32)*32)-16;
00090   if ( tsize < dim ) NAMD_bug("CudaNonbondedTables::buildVdwCoefTable bad tsize");
00091 
00092   float2 *h_vdwCoefTable = new float2[tsize*tsize];
00093   float2 *h_exclusionVdwCoefTable = new float2[tsize*tsize];
00094   float2 *row = h_vdwCoefTable;
00095   float2 *exclusionRow = h_exclusionVdwCoefTable;
00096   for ( int i=0; i<dim; ++i, row += tsize, exclusionRow += tsize ) {
00097     for ( int j=0; j<dim; ++j ) {
00098       const LJTable::TableEntry *e = ljTable->table_val(i,j);
00099       row[j].x = e->A * ComputeNonbondedUtil::scaling;
00100       row[j].y = e->B * ComputeNonbondedUtil::scaling;
00101       exclusionRow[j].x = ((e+1)->A - e->A)* ComputeNonbondedUtil::scaling;
00102       exclusionRow[j].y = ((e+1)->B - e->B)* ComputeNonbondedUtil::scaling;
00103     }
00104   }
00105 
00106   vdwCoefTableWidth = tsize;
00107 
00108   bindTextureObject<float2>(tsize*tsize, h_vdwCoefTable, vdwCoefTable, vdwCoefTableTex, update);
00109   bindTextureObject<float2>(tsize*tsize, h_exclusionVdwCoefTable, exclusionVdwCoefTable, exclusionVdwCoefTableTex, update);
00110 
00111   delete [] h_vdwCoefTable;
00112   delete [] h_exclusionVdwCoefTable;
00113 
00114   if ( ! CmiPhysicalNodeID(CkMyPe()) ) {
00115     CkPrintf("Info: Updated CUDA LJ table with %d x %d elements.\n", dim, dim);
00116   }
00117 }
00118 
00119 // Update tables from newer CPU values
00120 void CudaNonbondedTables::updateTables() {
00121   buildVdwCoefTable(true);
00122 }
00123 
00124 //
00125 // Builds force and energy tables. Swiped from ComputeNonbondedCUDA.C
00126 //
00127 template <typename T>
00128 void buildForceAndEnergyTable(const int tableSize, const double* r2list, const BigReal* src_table, const bool flip,
00129   const BigReal prefac, const int dst_stride, T* dst_force, T* dst_energy) {
00130 
00131   const BigReal r2_delta = ComputeNonbondedUtil:: r2_delta;
00132   const int r2_delta_exp = ComputeNonbondedUtil:: r2_delta_exp;
00133   const int r2_delta_expc = 64 * (r2_delta_exp - 1023);
00134 
00135   union { double f; int32 i[2]; } byte_order_test;
00136   byte_order_test.f = 1.0;  // should occupy high-order bits only
00137   int32 *r2iilist = (int32*)r2list + ( byte_order_test.i[0] ? 0 : 1 );
00138 
00139   for ( int i=1; i<tableSize; ++i ) {
00140     double r = ((double) tableSize) / ( (double) i + 0.5 );
00141     int table_i = (r2iilist[2*i] >> 14) + r2_delta_expc;  // table_i >= 0
00142 
00143     if ( r > ComputeNonbondedUtil::cutoff ) {
00144       dst_force[i*dst_stride] = 0.;
00145       dst_energy[i*dst_stride] = 0.;
00146       continue;
00147     }
00148 
00149     BigReal diffa = r2list[i] - ComputeNonbondedUtil::r2_table[table_i];
00150 
00151     BigReal table_a, table_b, table_c, table_d;
00152     if (flip) {
00153       table_a = src_table[4*table_i+3];
00154       table_b = src_table[4*table_i+2];
00155       table_c = src_table[4*table_i+1];
00156       table_d = src_table[4*table_i];
00157     } else {
00158       table_a = src_table[4*table_i];
00159       table_b = src_table[4*table_i+1];
00160       table_c = src_table[4*table_i+2];
00161       table_d = src_table[4*table_i+3];
00162     }
00163 
00164     BigReal grad = ( 3. * table_d * diffa + 2. * table_c ) * diffa + table_b;
00165     dst_force[i*dst_stride] = prefac * 2. * grad;
00166     BigReal ener = table_a + diffa * ( ( table_d * diffa + table_c ) * diffa + table_b);
00167     dst_energy[i*dst_stride] = prefac * ener;
00168   }
00169 
00170   dst_force[0] = 0.;
00171   dst_energy[0] = dst_energy[1*dst_stride];
00172 }
00173 
00174 template <typename T>
00175 void bindTextureObject(int tableSize, int tableWidth, T* h_table,
00176   cudaArray_t& array, cudaTextureObject_t& tableTex) {
00177 
00178   cudaChannelFormatDesc desc;
00179   memset(&desc, 0, sizeof(desc));
00180   desc.x = sizeof(T)*8;
00181   if (tableWidth >= 2) desc.y = sizeof(T)*8;
00182   if (tableWidth >= 3) desc.z = sizeof(T)*8;
00183   if (tableWidth >= 4) desc.w = sizeof(T)*8;
00184   desc.f = cudaChannelFormatKindFloat;
00185   cudaCheck(cudaMallocArray(&array, &desc, tableSize, 1));
00186   cudaCheck(cudaMemcpyToArray(array, 0, 0, h_table, tableSize*sizeof(T)*tableWidth, cudaMemcpyHostToDevice));
00187 
00188   cudaResourceDesc resDesc;
00189   memset(&resDesc, 0, sizeof(resDesc));
00190   resDesc.resType = cudaResourceTypeArray;
00191   resDesc.res.array.array = array;
00192 
00193   cudaTextureDesc texDesc;
00194   memset(&texDesc, 0, sizeof(texDesc));
00195   texDesc.addressMode[0] = cudaAddressModeClamp;
00196   texDesc.filterMode = cudaFilterModeLinear;
00197   texDesc.normalizedCoords = 1;
00198 
00199   cudaCheck(cudaCreateTextureObject(&tableTex, &resDesc, &texDesc, NULL));
00200 }
00201 
00202 void CudaNonbondedTables::buildForceAndEnergyTables(int tableSize) {
00203 
00204   // Build r2list
00205   const BigReal r2_delta = ComputeNonbondedUtil:: r2_delta;
00206   const int r2_delta_exp = ComputeNonbondedUtil:: r2_delta_exp;
00207   const int r2_delta_expc = 64 * (r2_delta_exp - 1023);
00208 
00209   double* r2list = new double[tableSize];  // double to match cpu code
00210   for ( int i=1; i<tableSize; ++i ) {
00211     double r = ((double) tableSize) / ( (double) i + 0.5 );
00212     r2list[i] = r*r + r2_delta;
00213   }
00214 
00215   // Non-bonded
00216   {
00217     float4* t = new float4[tableSize];
00218     float4* et = new float4[tableSize];
00219 
00220     buildForceAndEnergyTable<float>(tableSize, r2list, ComputeNonbondedUtil::fast_table, false, 1.0,
00221       4, &t[0].x, &et[0].x);
00222 
00223     buildForceAndEnergyTable<float>(tableSize, r2list, ComputeNonbondedUtil::vdwb_table, false, -1.0,
00224       4, &t[0].y, &et[0].y);
00225 
00226     buildForceAndEnergyTable<float>(tableSize, r2list, ComputeNonbondedUtil::vdwa_table, false, 1.0,
00227       4, &t[0].z, &et[0].z);
00228 
00229     buildForceAndEnergyTable<float>(tableSize, r2list, ComputeNonbondedUtil::scor_table, false, 1.0,
00230       4, &t[0].w, &et[0].w);
00231 
00232     bindTextureObject<float>(tableSize, 4, (float *)t, forceArray, forceTableTex);
00233     bindTextureObject<float>(tableSize, 4, (float *)et, energyArray, energyTableTex);
00234     delete [] t;
00235     delete [] et;
00236   }
00237 
00238   // Modified exclusions
00239   {
00240     float4* t = new float4[tableSize];
00241     float4* et = new float4[tableSize];
00242 
00243     buildForceAndEnergyTable<float>(tableSize, r2list, ComputeNonbondedUtil::fast_table, false, -1.0,
00244       4, &t[0].x, &et[0].x);
00245 
00246     buildForceAndEnergyTable<float>(tableSize, r2list, ComputeNonbondedUtil::vdwb_table, false, -1.0,
00247       4, &t[0].y, &et[0].y);
00248 
00249     buildForceAndEnergyTable<float>(tableSize, r2list, ComputeNonbondedUtil::vdwa_table, false, 1.0,
00250       4, &t[0].z, &et[0].z);
00251 
00252     buildForceAndEnergyTable<float>(tableSize, r2list, ComputeNonbondedUtil::slow_table, true, -1.0,
00253       4, &t[0].w, &et[0].w);
00254     // delete [] flip_slow_table;
00255 
00256     bindTextureObject<float>(tableSize, 4, (float *)t, modifiedExclusionForceArray, modifiedExclusionForceTableTex);
00257     bindTextureObject<float>(tableSize, 4, (float *)et, modifiedExclusionEnergyArray, modifiedExclusionEnergyTableTex);
00258     delete [] t;
00259     delete [] et;
00260   }
00261 
00262   // Exclusions
00263   {
00264     BigReal* corr_full_table = new BigReal[ComputeNonbondedUtil::table_length*4];
00265     for (int i=0;i < ComputeNonbondedUtil::table_length*4;i++) {
00266       corr_full_table[i] =  ComputeNonbondedUtil::corr_table[i] - ComputeNonbondedUtil::full_table[i];
00267     }
00268 
00269     float4* h_exclusionTable = new float4[ComputeNonbondedUtil::table_length];
00270     float* h_r2_table = new float[ComputeNonbondedUtil::table_length];
00271     for (int i=0;i < ComputeNonbondedUtil::table_length;i++) {
00272       h_exclusionTable[i].x = 6.0*corr_full_table[4*i + 3];
00273       h_exclusionTable[i].y = 4.0*corr_full_table[4*i + 2];
00274       h_exclusionTable[i].z = 2.0*corr_full_table[4*i + 1];
00275       h_exclusionTable[i].w = 1.0*corr_full_table[4*i + 0];
00276       h_r2_table[i] = ComputeNonbondedUtil::r2_table[i];
00277     }
00278 
00279     bindTextureObject<float>(ComputeNonbondedUtil::table_length, h_r2_table, r2_table, r2_table_tex);
00280     bindTextureObject<float4>(ComputeNonbondedUtil::table_length, h_exclusionTable, exclusionTable, exclusionTableTex);
00281 
00282     delete [] h_exclusionTable;
00283     delete [] h_r2_table;
00284 
00285   }
00286 
00287   if ( ! CmiPhysicalNodeID(CkMyPe()) ) {
00288     CkPrintf("Info: Updated CUDA force table with %d elements.\n", tableSize);
00289   }
00290 
00291   delete [] r2list;
00292 }
00293 
00294 #endif // NAMD_CUDA