DCS CUDA Block/Grid Decomposition (non-unrolled) Grid of thread blocks: Thread blocks: 0,0 0,1 … 64-256 threads 1,0 1,1 … … … … Threads compute 1 potential each Padding waste NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
DCS CUDA Block/Grid Decomposition (non-unrolled) • 16x16 CUDA thread blocks are a nice starting size with a satisfactory number of threads • Small enough that there’s not much waste due to padding at the edges NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
DCS Version 1: Const+Precalc 187 GFLOPS, 18.6 Billion Atom Evals/Sec (G80) • Pros: – Pre-compute dz^2 for entire slice – Inner loop over read-only atoms, const memory ideal – If all threads read the same const data at the same time, performance is similar to reading a register • Cons: – Const memory only holds ~4000 atom coordinates and charges – Potential summation must be done in multiple kernel invocations per slice, with const atom data updated for each invocation – Host must shuffle data in/out for each pass NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
DCS Version 1: Kernel Structure … float curenergy = energygrid[outaddr]; Start global memory reads float coorx = gridspacing * xindex; early. Kernel hides some of float coory = gridspacing * yindex; its own latency. int atomid; float energyval=0.0f; for (atomid=0; atomid<numatoms; atomid++) { float dx = coorx - atominfo[atomid].x; float dy = coory - atominfo[atomid].y; energyval += atominfo[atomid].w * rsqrtf(dx*dx + dy*dy + atominfo[atomid].z); } Only dependency on global memory read is at the end of energygrid[outaddr] = curenergy + energyval ; the kernel… NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
DCS CUDA Block/Grid Decomposition (unrolled, thread coarsening) • Reuse atom data and partial distance components multiple times • Use “unroll and jam” to unroll the outer loop into the inner loop • Uses more registers, but increases arithmetic intensity significantly • Kernels that unroll the inner loop calculate more than one lattice point per thread result in larger computational tiles: – Thread count per block must be decreased to reduce computational tile size as unrolling is increased – Otherwise, tile size gets bigger as threads do more than one lattice point evaluation, resulting on a significant increase in padding and wasted computations at edges NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
DCS CUDA Algorithm: Unrolling Loops • Add each atom’s contribution to several lattice points at a time, distances only differ in one component: potential[j ] += atom[i].charge / r ij potential[j+1] += atom[i].charge / r i(j+1) … Distances to Atom[i] Atom[i] NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
DCS CUDA Block/Grid Decomposition (unrolled) Unrolling increases Grid of thread blocks: computational tile size Thread blocks: 0,0 0,1 … 64-256 threads 1,0 1,1 … Threads compute … … … up to 8 potentials Padding waste NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
DCS Version 2: Inner Loop … for (atomid=0; atomid<numatoms; atomid++) { float dy = coory - atominfo[atomid].y; float dysqpdzsq = (dy * dy) + atominfo[atomid].z; float x = atominfo[atomid].x; Compared to non-unrolled float dx1 = coorx1 - x; kernel: memory loads are decreased by 4x, and FLOPS float dx2 = coorx2 - x; per evaluation are reduced, but float dx3 = coorx3 - x; register use is increased… float dx4 = coorx4 - x; float charge = atominfo[atomid].w; energyvalx1 += charge * rsqrtf(dx1*dx1 + dysqpdzsq); energyvalx2 += charge * rsqrtf(dx2*dx2 + dysqpdzsq); energyvalx3 += charge * rsqrtf(dx3*dx3 + dysqpdzsq); energyvalx4 += charge * rsqrtf(dx4*dx4 + dysqpdzsq); } NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
DCS Version 4: Const+Loop Unrolling+Coalescing 291.5 GFLOPS, 39.5 Billion Atom Evals/Sec (G80) • Pros: – Simplified structure compared to version 3, no use of shared memory, register pressure kept at bay by doing global memory operations only at the end of the kernel – Using fewer registers allows co-scheduling of more blocks, increasing GPU “occupancy” – Doesn’t have as strict of a thread block dimension requirement as version 3, computational tile size can be smaller • Cons: – The computation tile size is still large, so small potential maps don’t perform as well as large ones NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
DCS Version 4: Kernel Structure Processes 8 lattice points at a time in the inner • loop • Subsequent lattice points computed by each thread are offset by a half-warp to guarantee coalesced memory accesses • Loads and increments 8 potential map lattice points from global memory at completion of of the summation, avoiding register consumption Source code is available by request • NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
DCS Version 4: Inner Loop …float coory = gridspacing * yindex; float coorx = gridspacing * xindex; float gridspacing_coalesce = gridspacing * BLOCKSIZEX; Points spaced for int atomid; memory coalescing for (atomid=0; atomid<numatoms; atomid++) { float dy = coory - atominfo[atomid].y; float dyz2 = (dy * dy) + atominfo[atomid].z; Reuse partial distance components dy^2 + dz^2 float dx1 = coorx - atominfo[atomid].x; […] float dx8 = dx7 + gridspacing_coalesce; energyvalx1 += atominfo[atomid].w * rsqrtf(dx1*dx1 + dyz2); […] energyvalx8 += atominfo[atomid].w * rsqrtf(dx8*dx8 + dyz2); Global memory ops } occur only at the end energygrid[outaddr ] += energyvalx1; of the kernel, [...] decreases register use energygrid[outaddr+7*BLOCKSIZEX] += energyvalx7 ; NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
DCS CUDA Block/Grid Decomposition (unrolled, coalesced) Unrolling increases Grid of thread blocks: computational tile size Thread blocks: 0,0 0,1 … 64-256 threads 1,0 1,1 … … … … Threads compute up to 8 potentials, skipping by half-warps Padding waste NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Direct Coulomb Summation Performance Number of thread blocks modulo number of SMs results in significant performance variation for small workloads CUDA-Unroll8clx: fastest GPU kernel, 44x faster than CPU, 291 GFLOPS on GeForce 8800GTX CUDA-Simple: 14.8x faster, CPU 33% of fastest GPU kernel GPU computing. J. Owens, M. Houston, D. Luebke, S. Green, J. Stone, J. Phillips. Proceedings of the IEEE , 96:879-899, 2008. NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
DCS Version 4 Inner Loop, Scalar OpenCL … for (atomid=0; atomid<numatoms; atomid++) { float dy = coory - atominfo[atomid].y; Well-written CUDA code can float dyz2 = (dy * dy) + atominfo[atomid].z; often be easily ported to OpenCL if C++ features and pointer float dx1 = coorx – atominfo[atomid].x; arithmetic aren’t used in kernels. float dx2 = dx1 + gridspacing_coalesce; float dx3 = dx2 + gridspacing_coalesce; float dx4 = dx3 + gridspacing_coalesce; float charge = atominfo[atomid].w; energyvalx1 += charge * native_rsqrt(dx1*dx1 + dyz2); energyvalx2 += charge * native_rsqrt(dx2*dx2 + dyz2); energyvalx3 += charge * native_rsqrt(dx3*dx3 + dyz2); energyvalx4 += charge * native_rsqrt(dx4*dx4 + dyz2); } NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
DCS Version 4 Inner Loop (CUDA) (only 4-way unrolling for conciseness to compare OpenCL) … for (atomid=0; atomid<numatoms; atomid++) { float dy = coory - atominfo[atomid].y; float dyz2 = (dy * dy) + atominfo[atomid].z; float dx1 = coorx – atominfo[atomid].x; float dx2 = dx1 + gridspacing_coalesce; float dx3 = dx2 + gridspacing_coalesce; float dx4 = dx3 + gridspacing_coalesce; float charge = atominfo[atomid].w; energyvalx1 += charge * rsqrtf(dx1*dx1 + dyz2); energyvalx2 += charge * rsqrtf(dx2*dx2 + dyz2); energyvalx3 += charge * rsqrtf(dx3*dx3 + dyz2); energyvalx4 += charge * rsqrtf(dx4*dx4 + dyz2); } NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
DCS Version 4 Inner Loop, Vectorized OpenCL float4 gridspacing_u4 = { 0.f, 1.f, 2.f, 3.f }; CPUs, AMD GPUs, and Cell often perform gridspacing_u4 *= gridspacing_coalesce; better with vectorized kernels. float4 energyvalx=0.0f; Use of vector types may increase register pressure; sometimes a delicate balance… … for (atomid=0; atomid<numatoms; atomid++) { float dy = coory - atominfo[atomid].y; float dyz2 = (dy * dy) + atominfo[atomid].z; float4 dx = gridspacing_u4 + (coorx – atominfo[atomid].x); float charge = atominfo[atomid].w; energyvalx1 += charge * native_rsqrt(dx1*dx1 + dyz2); } NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Infinite vs. Cutoff Potentials • Infinite range potential: – All atoms contribute to all lattice points – Summation algorithm has quadratic complexity • Cutoff (range-limited) potential: – Atoms contribute within cutoff distance to lattice points – Summation algorithm has linear time complexity – Has many applications in molecular modeling: • Replace electrostatic potential with shifted form • Short-range part for fast methods of approximating full electrostatics • Used for fast decaying interactions (e.g. Lennard-Jones, Buckingham) NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Cutoff Summation • At each lattice point, sum potential contributions for atoms within cutoff radius: if (distance to atom[i] < cutoff) potential += (charge[i] / r) * s(r) • Smoothing function s(r) is algorithm dependent Cutoff radius r: distance to Atom[i] Lattice point being Atom[i] evaluated NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Cutoff Summation on the GPU Atoms spatially hashed into fixed- size “bins” in global memory Constant memory CPU handles overflowed bins Bin-Region Atoms neighborlist Global memory Potential Bins map regions of 8 Process atom bins for atoms current potential map region Shared memory Atom bin NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Using the CPU to Improve GPU Performance • GPU performs best when the work evenly divides into the number of threads/processing units • Optimization strategy: – Use the CPU to “regularize” the GPU workload – Handle exceptional or irregular work units on the CPU while the GPU processes the bulk of the work – On average, the GPU is kept highly occupied, attaining a much higher fraction of peak performance NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Cutoff Summation Runtime GPU cutoff with CPU overlap: 17x-21x faster than CPU core GPU acceleration of cutoff pair potentials for molecular modeling applications. C. Rodrigues, D. Hardy, J. Stone, K. Schulten, W. Hwu. Proceedings of the 2008 Conference On Computing Frontiers , pp. 273-282, 2008. NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Molecular Surface Visualization • Large biomolecular complexes are difficult to interpret with atomic detail graphical representations • Even secondary structure representations become cluttered • Surface representations are easier to use when greater abstraction is desired, but are computationally costly • Most surface display methods incapable of animating dynamics of large structures Poliovirus w/ millions of particles NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
VMD “QuickSurf” Representation • Displays continuum of structural detail: – All-atom models – Coarse-grained models – Cellular scale models – Multi-scale models: All-atom + CG, Brownian + Whole Cell – Smoothly variable between full detail, and reduced resolution representations of very large complexes Fast Visualization of Gaussian Density Surfaces for Molecular Dynamics and Particle System Trajectories. M. Krone, J. E. Stone, T. Ertl, K. Schulten. EuroVis Short Papers , pp. 67-71, 2012 NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
VMD “QuickSurf” Representation • Uses multi-core CPUs and GPU acceleration to enable smooth real-time animation of MD trajectories • Linear-time algorithm, scales to millions of particles, as limited by memory capacity Satellite Tobacco Mosaic Virus Lattice Cell Simulations NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
VMD “QuickSurf” Representation All-atom HIV capsid simulations NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
QuickSurf Representation of Lattice Cell Models Continuous particle Discretized lattice models derived based model – often 70 from continuous model shown in to 300 million particles VMD QuickSurf representation Lattice Microbes: High ‐ performance stochastic simulation method for the reaction ‐ diffusion master equation E. Roberts, J. E. Stone, and Z. Luthey ‐ Schulten. J. Computational Chemistry 34 (3), 245-255, 2013. NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
QuickSurf Algorithm Overview • Build spatial acceleration data structures, optimize data for GPU • Compute 3-D density map, 3-D volumetric texture map: 3-D density map lattice, spatial acceleration grid, and extracted surface • Extract isosurface for a user-defined density value NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
QuickSurf Particle Sorting, Bead Generation, Spatial Hashing • Particles sorted into spatial acceleration grid: – Selected atoms or residue “beads” converted lattice coordinate system – Each particle/bead assigned cell index, sorted w/NVIDIA Thrust template library • Complication: – Thrust allocates GPU mem. on-demand, no recourse if insufficient memory, have to re-gen QuickSurf data structures if caught by surprise! Coarse resolution spatial acceleration grid • Workaround: – Pre-allocate guesstimate workspace for Thrust – Free the Thrust workspace right before use – Newest Thrust allows user- defined allocator code… NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Spatial Hashing Algorithm Steps/Kernels 1) Compute bin index for each atom, store to memory w/ atom index 2) Sort list of bin and atom index tuples (1) by bin index (thrust kernel) 3) Count atoms in each bin (2) using a parallel prefix sum, aka scan , compute the destination index for each QuickSurf uniform atom, store per-bin starting index and grid spatial atom count (thrust kernel) subdivision data structure 4) Write atoms to the output indices computed in (3), and we have completed the data structure NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
QuickSurf and Limited GPU Global Memory • High resolution molecular surfaces require a fine lattice spacing • Memory use grows cubically with decreased lattice spacing • Not typically possible to compute a surface in a single pass, so we loop over sub- volume “chunks” until done… • Chunks pre-allocated and sized to GPU global mem capacity to prevent unexpected memory allocation failure while animating… • Complication: – Thrust allocates GPU mem. on-demand, no recourse if insufficient memory, have to re-gen QuickSurf data structures if caught by surprise! • Workaround: – Pre-allocate guesstimate workspace for Thrust – Free the Thrust workspace right before use – Newest Thrust allows user- defined allocator code… NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
QuickSurf Density Parallel Decomposition … QuickSurf 3-D density map Chunk 2 decomposes into thinner 3-D Chunk 1 slabs/slices (CUDA grids) Chunk 0 Large volume Small 8x8 thread computed in blocks afford large multiple passes, or per-thread register multiple GPUs count, shared memory … 0,0 0,1 Threads Each thread producing results that computes … 1,0 1,1 are used one or more density map … … … lattice points Inactive threads, Padding optimizes global region of memory performance, discarded output guaranteeing coalesced Grid of thread blocks global memory accesses NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
QuickSurf Density Map Algorithm • Spatial acceleration grid cells are sized to match the cutoff radius for the exponential, beyond which density contributions are negligible • Density map lattice points computed by summing density contributions from particles in 3x3x3 grid of neighboring spatial acceleration cells 3-D density map • Volumetric texture map is computed lattice point and by summing particle colors the neighboring normalized by their individual density spatial acceleration cells it references contribution NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
QuickSurf Density Map Kernel Optimizations • Compute reciprocals, prefactors, other math on the host CPU prior to kernel launch • Use of intN and floatN vector types in CUDA kernels for improved global memory bandwidth • Thread coarsening : one thread computes multiple output densities and colors • Input data and register tiling : share blocks of input, partial distances in regs shared among multiple outputs • Global memory (L1 cache) broadcasts : all threads in the block traverse the same atom/particle at the same time NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
QuickSurf Density Map Kernel Snippet… for (zab=zabmin; zab<=zabmax; zab++) { for (yab=yabmin; yab<=yabmax; yab++) { for (xab=xabmin; xab<=xabmax; xab++) { int abcellidx = zab * acplanesz + yab * acncells.x + xab; uint2 atomstartend = cellStartEnd[abcellidx]; if (atomstartend.x != GRID_CELL_EMPTY) { for (unsigned int atomid=atomstartend.x; atomid<atomstartend.y; atomid++) { float4 atom = sorted_xyzr[atomid]; float dx = coorx - atom.x; float dy = coory - atom.y; float dz = coorz - atom.z; float dxy2 = dx*dx + dy*dy; float r21 = (dxy2 + dz*dz) * atom.w; densityval1 += exp2f(r21); /// Loop unrolling and register tiling benefits begin here…… float dz2 = dz + gridspacing; float r22 = (dxy2 + dz2*dz2) * atom.w; densityval2 += exp2f(r22); /// More loop unrolling …. NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
QuickSurf Marching Cubes Isosurface Extraction • Isosurface is extracted from each density map “chunk”, and either copied back to the host, or rendered directly out of GPU global memory via CUDA/OpenGL interop • All MC memory buffers are pre-allocated to prevent significant overhead when animating a simulation trajectory … QuickSurf 3-D density map Chunk 2 decomposes into thinner 3-D Chunk 1 slabs/slices (CUDA grids) Chunk 0 Large volume computed in multiple passes NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Brief Marching Cubes Isosurface Extraction Overview • Given a 3-D volume of scalar density values and a requested surface density value, marching cubes computes vertices and triangles that compose the requested surface triangle mesh • Each MC “cell” (a cube with 8 density values at its vertices) produces a variable number of output vertices depending on how many edges of the cell contain the requested isovalue… • Use scan() to compute the output indices so that each worker thread has conflict-free output of vertices/triangles NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Brief Marching Cubes Isosurface Extraction Overview • Once the output vertices have been computed and stored, we compute surface normals and colors for each of the vertices • Although the separate normals+colors pass reads the density map again, molecular surfaces tend to generate a small percentage of MC cells containing triangles, we avoid wasting interpolation work • We use CUDA tex3D() hardware 3-D texture mapping: – Costs double the texture memory and a one copy from GPU global memory to the target texture map with cudaMemcpy3D() – Still roughly 2x faster than doing color interpolation without the texturing hardware, at least on GT200 and Fermi hardware – Kepler has new texture cache memory path that may make it feasible to do our own color interpolation and avoid the use of extra 3-D texture memory and associated copy, with acceptable performance NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
QuickSurf Marching Cubes Isosurface Extraction • Our optimized MC implementation computes per-vertex surface normals, colors, and outperforms the NVIDIA SDK sample by a fair margin on Fermi GPUs • Complications: – Even on a 6GB Quadro 7000, GPU global memory is under great strain when working with large molecular complexes, e.g. viruses – Marching cubes involves a parallel prefix sum (scan) to compute target indices for writing resulting vertices – We use Thrust for scan, has the same memory allocation issue mentioned earlier for the sort, so we use the same workaround – The number of output vertices can be huge, but we rarely have sufficient GPU memory for this – we use a fixed size vertex output buffer and hope our heuristics don’t fail us NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
QuickSurf Performance GeForce GTX 580 Molecular Atoms Resolution T sort T density T MC # vertices FPS system MscL 111,016 1.0Å 0.005 0.023 0.003 0.7 M 28 STMV capsid 147,976 1.0Å 0.007 0.048 0.009 2.4 M 13.2 Poliovirus 754,200 1.0Å 0.01 0.18 0.05 9.2 M 3.5 capsid STMV w/ water 955,225 1.0Å 0.008 0.189 0.012 2.3 M 4.2 Membrane 2.37 M 2.0Å 0.03 0.17 0.016 5.9 M 3.9 Chromatophore 9.62 M 2.0Å 0.16 0.023 0.06 11.5 M 3.4 Membrane w/ 22.77 M 4.0Å 4.4 0.68 0.01 1.9 M 0.18 water Fast Visualization of Gaussian Density Surfaces for Molecular Dynamics and Particle System Trajectories. M. Krone, J. E. Stone, T. Ertl, K. Schulten. EuroVis Short Papers , pp. 67-71, 2012 NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Extensions and Analysis Uses for QuickSurf Triangle Mesh • Curved PN triangles: – We have performed tests with post-processing the resulting triangle mesh and using curved PN triangles to generate smooth surfaces with a larger grid spacing, for increased performance – Initial results demonstrate some potential, but there can be pathological cases where MC generates long skinny triangles, causing unsightly surface creases • Analysis uses (beyond visualization): – Minor modifications to the density map algorithm allow rapid computation of solvent accessible surface area by summing the areas in the resulting triangle mesh – Modifications to the density map algorithm will allow it to be used for MDFF (molecular dynamics flexible fitting) – Surface triangle mesh can be used as the input for computing the electrostatic potential field for mesh-based algorithms NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Challenge: Support Interactive QuickSurf for Large Structures on Mid-Range GPUs • Structures such as HIV initially needed large (6GB) GPU memory to generate fully-detailed surface renderings • Goals and approach: – Avoid slow CPU-fallback! – Incrementally change algorithm phases to use more compact data types, while maintaining performance – Specialize code for different performance/memory capacity cases NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Improving QuickSurf Memory Efficiency • Both host and GPU memory capacity limitations are a significant concern when rendering surfaces for virus structures such as HIV or for large cellular models which can contain hundreds of millions of particles • The original QuickSurf implementation used single- precision floating point for output vertex arrays and textures • Judicious use of reduced-precision numerical representations, cut the overall memory footprint of the entire QuickSurf algorithm to half of the original – Data type changes made throughout the entire chain from density map computation through all stages of Marching Cubes NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Supporting Multiple Data Types for QuickSurf Density Maps and Marching Cubes Vertex Arrays • The major algorithm components of QuickSurf are now used for many other purposes: – Gaussian density map algorithm now used for MDFF Cryo EM density map fitting methods in addition to QuickSurf – Marching Cubes routines also used for Quantum Chemistry visualizations of molecular orbitals • Rather than simply changing QuickSurf to use a particular internal numerical representation, it is desirable to instead use CUDA C++ templates to make type-generic versions of the key objects, kernels, and output vertex arrays • Accuracy-sensitive algorithms use high-precision data types, performance and memory capacity sensitive cases use quantized or reduced precision approaches NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Minimizing the Impact of Generality on QuickSurf Code Complexity • A critical factor in the simplicity of supporting multiple QuickSurf data types arises from the so-called “gather” oriented algorithm we employ – Internally, all in-register arithmetic is single-precision – Data conversions to/from compressed or reduced precision data types are performed on-the-fly as needed • Small inlined type conversion routines are defined for each of the cases we want to support • Key QuickSurf kernels are genericized using C++ template syntax, and the compiler “connects the dots” to automatically generate type-specific kernels as needed NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Example Templated Density Map Kernel template<class DENSITY, class VOLTEX> __global__ static void gaussdensity_fast_tex_norm(int natoms, const float4 * RESTRICT sorted_xyzr, const float4 * RESTRICT sorted_color, int3 numvoxels, int3 acncells, float acgridspacing, float invacgridspacing, const uint2 * RESTRICT cellStartEnd, float gridspacing, unsigned int z, DENSITY * RESTRICT densitygrid, VOLTEX * RESTRICT voltexmap, float invisovalue) { NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Example Templated Density Map Kernel template<class DENSITY, class VOLTEX> __global__ static void gaussdensity_fast_tex_norm( … ) { … Triple - nested and unrolled inner loops here … DENSITY densityout; VOLTEX texout; convert_density(densityout, densityval1); densitygrid[outaddr ] = densityout; convert_color(texout, densitycol1); voltexmap[outaddr ] = texout; NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Net Result of QuickSurf Memory Efficiency Optimizations • Halved overall GPU memory use • Achieved 1.5x to 2x performance gain : – The “gather” density map algorithm keeps type conversion operations out of the innermost loop – Density map global memory writes reduced to half – Multiple stages of Marching Cubes operate on smaller input and output data types – Same code path supports multiple precisions • Users now get full GPU-accelerated QuickSurf in many cases that previously triggered CPU- fallback, all platforms (laptop/desk/super) benefit! NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
High Resolution HIV Surface NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Structural Route to the all-atom HIV-1 Capsid 1st TEM (1999) 1st tomography (2003) Crystal structures of separated hexamer and pentamer Pornillos et al. , Cell 2009 , Nature 2011 Ganser et al. Science , 1999 Briggs et al. EMBO J , 2003 High res. EM of hexameric tubule, tomography of capsid, all-atom model of capsid by MDFF w/ NAMD & VMD, Briggs et al. Structure , 2006 cryo-ET (2006) NSF/NCSA Blue Waters computer at Illinois hexameric tubule Zhao et al. , Nature 497: 643-646 (2013 ) Li et al., Nature , 2000 NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, Byeon et al., Cell 2009 U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Molecular Dynamics Flexible Fitting (MDFF) X-ray crystallography MDFF Electron microscopy APS at Argonne FEI microscope ORNL Titan Acetyl - CoA Synthase Flexible fitting of atomic structures into electron microscopy maps using molecular dynamics. NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/ L. Trabuco, E. Villa, K. Mitra, J. Frank, and K. Schulten. Structure, 16:673-683, 2008.
Evaluating Quality-of-Fit for Structures Solved by Hybrid Fitting Methods Compute Pearson correlation to evaluate the fit of a reference cryo-EM density map with a simulated density map produced from an all-atom structure. NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
GPUs Can Reduce Trajectory Analysis Runtimes from Hours to Minutes GPUs enable laptops and desktop workstations to handle tasks that would have previously required a cluster, or a very long wait… GPU-accelerated petascale supercomputers enable analyses were previously impossible, allowing detailed study of very large structures such as viruses GPU-accelerated MDFF Cross Correlation Timeline Regions with poor fit Regions with good fit NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Single-Pass MDFF GPU Cross-Correlation 3-D density map decomposes into 3-D grid of 8x8x8 tiles containing CC partial sums and local CC values Spatial CC map and overall CC value computed in a single pass Small 8x8x2 CUDA thread blocks afford large per-thread register count, shared memory Each thread computes … 0,0 0,1 Threads 4 z-axis density map producing lattice points and results that … 1,0 1,1 associated CC partial are used sums … … … Inactive threads, region of Padding optimizes global discarded memory performance, output guaranteeing coalesced global Grid of thread blocks memory accesses NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
VMD GPU Cross Correlation Performance RHDV Mm-cpn GroEL Aquaporin open Resolution (Å) 6.5 8 4 3 Atoms 702K 61K 54K 1.6K VMD-CUDA 0.458s 0.06s 0.034s 0.007s Quadro K6000 34.6x 25.7x 36.8x 55.7x VMD-CPU-SSE 0.779s 0.085s 0.159s 0.033s 32-threads, 2x Xeon E5-2687W 20.3x 18.1x 7.9x 11.8x Chimera 15.86s 1.54s 1.25s 0.39s 1-thread Xeon E5-2687W 1.0x 1.0x 1.0x 1.0x GPU-accelerated analysis and visualization of large structures solved by molecular dynamics flexible fitting. J. E. Stone, R. McGreevy, B. Isralewitz, and K. Schulten. Faraday Discussion 169, 2014. (In press). NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
VMD RHDV Cross Correlation Timeline on Cray XK7 RHDV Atoms 702K RHDV CC Timeline Component 720 Selections Single-node XK7 336 hours (14 days) (projected) 128-node XK7 3.2 hours 105x speedup Calculation would take 5 years using conventional non-GPU software on a workstation!! NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Molecular Orbitals • Visualization of MOs aids in understanding the chemistry of molecular system • MO spatial distribution is correlated with probability density for an electron(s) • Algorithms for computing other molecular properties are similar, and can share code High Performance Computation and Interactive Display of Molecular Orbitals on GPUs and Multi-core CPUs. J. Stone, J. Saam, D. Hardy, K. Vandivort, W. Hwu, K. Schulten, 2nd Workshop on General-Purpose Computation on Graphics Pricessing Units (GPGPU-2), ACM International Conference Proceeding Series , volume 383, pp. 9-18, 2009. NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Computing Molecular Orbitals • Calculation of high resolution MO grids can require tens to hundreds of seconds in existing tools • Existing tools cache MO grids as much as possible to avoid recomputation: – Doesn’t eliminate the wait for initial calculation, hampers interactivity – Cached grids consume C 60 100x-1000x more memory than MO coefficients NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Animating Molecular Orbitals • Animation of (classical mechanics) molecular dynamics trajectories provides insight into simulation results • To do the same for QM or QM/MM simulations one must compute MOs at ~ 10 FPS or more • >100x speedup (GPU) over existing tools now makes C 60 this possible! NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Molecular Orbital Computation and Display Process One-time Read QM simulation log file, trajectory initialization Initialize Pool of GPU Preprocess MO coefficient data Worker Threads eliminate duplicates, sort by type, etc… For current frame and MO index, retrieve MO wavefunction coefficients Compute 3-D grid of MO wavefunction amplitudes Most performance-demanding step, run on GPU… For each trj frame, for each MO shown Extract isosurface mesh from 3-D MO grid Apply user coloring/texturing and render the resulting surface NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
MO GPU Parallel Decomposition MO 3-D lattice … decomposes into 2-D GPU 2 slices (CUDA grids) GPU 1 GPU 0 Small 8x8 thread Lattice can be blocks afford large computed using per-thread register multiple GPUs count, shared memory … 0,0 0,1 Threads producing Each thread results that … 1,0 1,1 computes are used one MO lattice point. … … … Threads Padding optimizes global producing memory performance, results that are discarded guaranteeing coalesced Grid of thread blocks global memory accesses NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
MO GPU Kernel Snippet: Contracted GTO Loop, Use of Constant Memory [… outer loop over atoms …] float dist2 = xdist2 + ydist2 + zdist2; // Loop over the shells belonging to this atom (or basis function) for (shell=0; shell < maxshell; shell++) { float contracted_gto = 0.0f; // Loop over Gaussian primitives of this contracted basis function to build the atomic orbital int maxprim = const_num_prim_per_shell[shell_counter]; Constant memory: int shelltype = const_shell_types[shell_counter]; nearly register- for (prim=0; prim < maxprim; prim++) { float exponent = const_basis_array[prim_counter ]; speed when array float contract_coeff = const_basis_array[prim_counter + 1]; elements accessed contracted_gto += contract_coeff * __expf(-exponent*dist2); in unison by all prim_counter += 2; threads…. } [… continue on to angular momenta loop …] NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
MO GPU Kernel Snippet: Unrolled Angular Momenta Loop /* multiply with the appropriate wavefunction coefficient */ float tmpshell=0; switch (shelltype) { Loop unrolling: case S_SHELL: value += const_wave_f[ifunc++] * contracted_gto; • Saves registers break; [… P_SHELL case …] (important for GPUs!) case D_SHELL: tmpshell += const_wave_f[ifunc++] * xdist2; • Reduces loop control tmpshell += const_wave_f[ifunc++] * xdist * ydist; tmpshell += const_wave_f[ifunc++] * ydist2; overhead tmpshell += const_wave_f[ifunc++] * xdist * zdist; tmpshell += const_wave_f[ifunc++] * ydist * zdist; tmpshell += const_wave_f[ifunc++] * zdist2; • Increases arithmetic value += tmpshell * contracted_gto; intensity break; [... Other cases: F_SHELL, G_SHELL, etc …] } // end switch NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Preprocessing of Atoms, Basis Set, and Wavefunction Coefficients • Must make effective use of high bandwidth, low- latency GPU on-chip shared memory, or L1 cache: – Overall storage requirement reduced by eliminating duplicate basis set coefficients – Sorting atoms by element type allows re-use of basis set coefficients for subsequent atoms of identical type • Padding, alignment of arrays guarantees coalesced GPU global memory accesses NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
GPU Traversal of Atom Type, Basis Set, Shell Type, and Wavefunction Coefficients Monotonically increasing memory references Constant for all MOs, all timesteps Different at each timestep, and for Strictly sequential memory references each MO • Loop iterations always access same or consecutive array elements for all threads in a thread block: – Yields good constant memory and L1 cache performance – Increases shared memory tile reuse NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Use of GPU On-chip Memory • If total data less than 64 kB, use only const mem: – Broadcasts data to all threads, no global memory accesses! • For large data, shared memory used as a program-managed cache, coefficients loaded on-demand: – Tiles sized large enough to service entire inner loop runs, broadcast to all 64 threads in a block – Complications: nested loops, multiple arrays, varying length – Key to performance is to locate tile loading checks outside of the two performance-critical inner loops – Only 27% slower than hardware caching provided by constant memory (on GT200) • Fermi/Kepler GPUs have larger on-chip shared memory, L1/L2 caches, greatly reducing control overhead NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Array tile loaded in GPU shared memory. Tile size is a power-of-two, a multiple of coalescing size, and allows simple indexing in inner loops. Global memory array indices are merely offset to reference an MO coefficient within a tile loaded in fast on-chip shared memory. Surrounding data, unreferenced by next batch of loop iterations 64-byte memory coalescing block boundaries Full tile padding MO coefficient array in GPU global memory. Tiles are referenced in consecutive order. NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
VMD MO GPU Kernel Snippet: Loading Tiles Into Shared Memory On-Demand [… outer loop over atoms …] Shared memory tiles: if ((prim_counter + (maxprim<<1)) >= SHAREDSIZE) { prim_counter += sblock_prim_counter; • Tiles are checked sblock_prim_counter = prim_counter & MEMCOAMASK; s_basis_array[sidx ] = basis_array[sblock_prim_counter + sidx ]; and loaded, if s_basis_array[sidx + 64] = basis_array[sblock_prim_counter + sidx + 64]; necessary, s_basis_array[sidx + 128] = basis_array[sblock_prim_counter + sidx + 128]; immediately prior to s_basis_array[sidx + 192] = basis_array[sblock_prim_counter + sidx + 192]; prim_counter -= sblock_prim_counter; entering key __syncthreads(); arithmetic loops } for (prim=0; prim < maxprim; prim++) { float exponent = s_basis_array[prim_counter ]; • Adds additional float contract_coeff = s_basis_array[prim_counter + 1]; control overhead to contracted_gto += contract_coeff * __expf(-exponent*dist2); loops, even with prim_counter += 2; } optimized [… continue on to angular momenta loop …] implementation NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
New GPUs Bring Opportunities for Higher Performance and Easier Programming • NVIDIA’s Fermi, Kepler, Maxwell GPUs bring: – Greatly increased peak single- and double-precision arithmetic rates – Moderately increased global memory bandwidth – Increased capacity on-chip memory partitioned into shared memory and an L1 cache for global memory – Concurrent kernel execution – Bidirectional asynchronous host-device I/O – ECC memory, faster atomic ops, many others… NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
VMD MO GPU Kernel Snippet: Kernel based on L1 cache (Fermi) or Read-only Data Cache (Maxwell) [… outer loop over atoms …] L1 cache: // loop over the shells/basis funcs belonging to this atom for (shell=0; shell < maxshell; shell++) { • Simplifies code! float contracted_gto = 0.0f; • Reduces control int maxprim = shellinfo[(shell_counter<<4) ]; int shell_type = shellinfo[(shell_counter<<4) + 1]; overhead for (prim=0; prim < maxprim; prim++) { • Gracefully handles float exponent = basis_array[prim_counter ]; arbitrary-sized float contract_coeff = basis_array[prim_counter + 1]; problems contracted_gto += contract_coeff * __expf(- exponent*dist2); • Matches performance prim_counter += 2; } of constant memory on [… continue on to angular momenta loop …] Fermi and Maxwell NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
MO Kernel for One Grid Point (Naive C) … for (at=0; at<numatoms; at++) { Loop over atoms int prim_counter = atom_basis[at]; calc_distances_to_atom(&atompos[at], &xdist, &ydist, &zdist, &dist2, &xdiv); for (contracted_gto=0.0f, shell=0; shell < num_shells_per_atom[at]; shell++) { Loop over shells int shell_type = shell_symmetry[shell_counter]; for (prim=0; prim < num_prim_per_shell[shell_counter]; prim++) { Loop over primitives: float exponent = basis_array[prim_counter ]; largest component of float contract_coeff = basis_array[prim_counter + 1]; runtime, due to expf() contracted_gto += contract_coeff * expf(-exponent*dist2); prim_counter += 2; } Loop over angular for (tmpshell=0.0f, j=0, zdp=1.0f; j<=shell_type; j++, zdp*=zdist) { momenta int imax = shell_type - j; for (i=0, ydp=1.0f, xdp=pow(xdist, imax); i<=imax; i++, ydp*=ydist, xdp*=xdiv) (unrolled in real code) tmpshell += wave_f[ifunc++] * xdp * ydp * zdp; } value += tmpshell * contracted_gto; shell_counter++; } } ….. NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
VMD MO Performance Results for C 60 Sun Ultra 24: Intel Q6600, NVIDIA GTX 280 Kernel Cores/GPUs Runtime (s) Speedup CPU ICC-SSE 1 46.58 1.00 CPU ICC-SSE 4 11.74 3.97 CPU ICC-SSE-approx** 4 3.76 12.4 CUDA-tiled-shared 1 0.46 100. CUDA-const-cache 1 0.37 126. C 60 basis set 6-31Gd. We used an unusually-high resolution MO grid for accurate timings. A more typical calculation has 1/8 th the grid points. **Reduced-accuracy approximation of expf(), cannot be used for zero-valued MO isosurfaces NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
VMD Single-GPU Molecular Orbital Performance Results for C 60 on Fermi Intel X5550 CPU, GeForce GTX 480 GPU Kernel Cores/GPUs Runtime (s) Speedup Xeon 5550 ICC-SSE 1 30.64 1.0 Xeon 5550 ICC-SSE 8 4.13 7.4 CUDA shared mem 1 0.37 83 CUDA L1-cache (16KB) 1 0.27 113 CUDA const-cache 1 0.26 117 CUDA const-cache, zero-copy 1 0.25 122 Fermi GPUs have caches: match perf. of hand-coded shared memory kernels. Zero-copy memory transfers improve overlap of computation and host-GPU I/Os. NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Preliminary Single-GPU Molecular Orbital Performance Results for C 60 on Kepler Intel X5550 CPU, GeForce GTX 680 GPU Kernel Cores/GPUs Runtime (s) Speedup Xeon 5550 ICC-SSE 1 30.64 1.0 Xeon 5550 ICC-SSE 8 4.13 7.4 CUDA shared mem 1 0.264 116 CUDA RO-data-cache 1 0.228 134 CUDA const-cache 1 0.104 292 CUDA const-cache, zero-copy 1 0.0938 326 Kepler GK104 (GeForce 680) strongly prefers the constant cache kernels vs. the others. NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
VMD Orbital Dynamics Proof of Concept One GPU can compute and animate this movie on-the-fly! CUDA const-cache kernel, Sun Ultra 24, GeForce GTX 285 0.016 s GPU MO grid calc. CPU surface gen, 0.033 s volume gradient, and GPU rendering Total runtime 0.049 s threonine Frame rate 20 FPS With GPU speedups over 100x , previously insignificant CPU surface gen, gradient calc, and rendering are now 66% of runtime. Needed GPU- accelerated surface gen next… Wrote CUDA Marching Cubes to address surface gen perf gap. NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Multi-GPU Load Balance • Many early CUDA codes assumed all GPUs were identical • Host machines may contain a diversity of GPUs of varying capability (discrete, IGP, etc) • Different GPU on-chip and global memory capacities may need different problem “tile” sizes GPU 1 GPU N • Static decomposition works … poorly for non-uniform workload, 14 SMs 30 SMs or diverse GPUs NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
MO GPU Parallel Decomposition MO 3-D lattice … decomposes into 2-D GPU 2 slices (CUDA grids) GPU 1 GPU 0 Small 8x8 thread Lattice can be blocks afford large computed using per-thread register multiple GPUs count, shared memory … 0,0 0,1 Threads producing Each thread results that … 1,0 1,1 computes are used one MO lattice point. … … … Threads Padding optimizes global producing memory performance, results that are discarded guaranteeing coalesced Grid of thread blocks global memory accesses NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Multi-GPU Dynamic Work Distribution Dynamic work distribution // Each GPU worker thread loops over // subset of work items… while (!threadpool_next_tile(&parms, tilesize, &tile){ // Process one work item… // Launch one CUDA kernel for each // loop iteration taken… // Shared iterator automatically GPU 1 GPU N … // balances load on GPUs } NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Example Multi-GPU Latencies 4 C2050 GPUs, Intel Xeon 5550 6.3us CUDA empty kernel (immediate return) 9.0us Sleeping barrier primitive (non-spinning barrier that uses POSIX condition variables to prevent idle CPU consumption while workers wait at the barrier) 14.8us pool wake, host fctn exec, sleep cycle (no CUDA) 30.6us pool wake, 1x(tile fetch, simple CUDA kernel launch), sleep 1817.0us pool wake, 100x(tile fetch, simple CUDA kernel launch), sleep NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Multi-GPU Runtime Original Error/Exception Handling Workload • Competition for resources from other applications can cause runtime failures, e.g. Retry Stack GPU out of memory half way through an algorithm • Handle exceptions, e.g. convergence failure, NaN result, insufficient compute GPU 1 GPU N capability/features … SM 1.1 SM 2.0 • Handle and/or reschedule 128MB 3072MB failed tiles of work NIH BTRC for Macromolecular Modeling and Bioinformatics Beckman Institute, U. Illinois at Urbana-Champaign http://www.ks.uiuc.edu/
Recommend
More recommend
