# Presentation: 515105

## Leading Computational Methods on Modern Parallel Systems

- Leonid Oliker
- Lawrence Berkeley National Laboratory
- Joint work with:
- Julian Borrill, Andrew Canning,  Stephane Ethier, Art Mirin, David Parks, John Shalf, David Skinner, Michael Wehner, Patrick Worley

## Motivation

- Stagnating application performance is well-know problem in scientific computing
- By end of decade numerous mission critical applications expected to have 100X computational demands of current levels
- Many HEC platforms are poorly balanced for demands of leading applications
- Memory-CPU gap, deep memory hierarchies, poor network-processor integration, low-degree network topology
- Traditional superscalar trends slowing down
- Mined most benefits of ILP and pipelining,Clock frequency limited by power concerns
- Achieving next order of magnitude in computing will require O(100K) - O(1M) processors: interconnect that scale nonlinearly in cost (crossbar, fattree) become impractical

## Application Evaluation

- Microbenchmarks, algorithmic kernels, performance modeling and prediction, are important components of understanding and improving architectural efficiency
- However full-scale application performance is the final arbiter of system utility and necessary as baseline to support all complementary approaches
- Our evaluation work emphasizes full applications, with real input data, at the appropriate scale
- Requires coordination of computer scientists and application experts from highly diverse backgrounds
- Allows advanced optimizations for given architectural platforms
- In-depth studies reveal limitation of compilers, operating systems, and hardware, since all of these components must work together at scale
- Our recent efforts have focused on the potential of modern parallel vector systems
- Effective code vectorization is an integral part of the process
- First US team to conduct Earth Simulator performance study

## Vector Paradigm

- High memory bandwidth
- Allows systems to effectively feed ALUs (high byte to flop ratio)
- Flexible memory addressing modes
- Supports fine grained strided and irregular data access
- Vector Registers
- Hide memory latency via deep pipelining of memory load/stores
- Vector ISA
- Single instruction specifies large number of identical operations
- Vector architectures allow for:
- Reduced control complexity
- Efficiently utilize large number of computational resources
- Potential for automatic discovery of parallelism
- Can some of these features can be added to modern micros?
*However: most effective if sufficient regularity discoverable in program structure*
- Suffers even if small % of code non-vectorizable (Amdahl’s Law)

## Architectural Comparison

| Node Type | Where | Network | CPU/Node | ClockMHz | PeakGFlop | Stream BW GB/s/P | Peak byte/flop | MPIBWGB/s/P | MPI Latencysec | NetworkTopology |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Power3 | NERSC | Colony | 16 | 375 | 1.5 | 0.4 | 0.26 | 0.13 | 16.3 | Fat-tree |
| Itanium2 | LLNL | Quadrics | 4 | 1400 | 5.6 | 1.1 | 0.19 | 0.25 | 3.0 | Fat-tree |
| Opteron | NERSC | InfiniBand | 2 | 2200 | 4.4 | 2.3 | 0.51 | 0.59 | 6.0 | Fat-tree |
| X1 | ORNL | Custom | 4 | 800 | 12.8 | 14.9  | 1.16 | 6.3 | 7.1 | 4D-Hypercube |
| X1E | ORNL | Custom | 4 | 1130 | 18.0 | 9.7  | 0.54 | 2.9 | 5.0 | 4D-Hypercube |
| ES | ESC | IN | 8 | 1000 |  8.0 | 26.3 | 3.29 | 1.5 | 5.6 | Crossbar |
| SX-8 | HLRS | INX | 8 | 2000 | 16.0 | 41.0 | 2.56 | 2.0 | 5.0 | Crossbar |

-  Custom vectors : High Bandwidth, Flexible addressing, Many Register, vector ISA
-  Less control complexity, use many ALUs, auto parallelism
-  ES shows best balance between memory and peak performance
-  Data caches of superscalar systems and X1(E) potential reduce mem costs
-  X1E: 2 MSP’s per MCM - increases contention for memory and interconnect
-  A key ‘balance point’ for vector systems is the scalar:vector ratio (8:1 or even 32:1)
-  Opteron/IB shows best balance for superscalar, Itanium2/Quadrics lowest latency
-  Cost is a critical metric - however we are unable to provide such data
-  Proprietary, pricing varies based on customer and time frame
-  Poorly balanced systems cannot solve certain class of problems/resolutions regardless of processor count!

## Application Overview

| NAME | Discipline | Problem/Method | Structure |
| --- | --- | --- | --- |
| Cactus | Astrophysics | Theory of GRADM-BSSN, Method of Lines | Grid |
| LBMHD | Plasma Physics | Magneto-Hydrodyamics, Lattice-Boltzmann | Lattice/Grid |
| MADCAP | Cosmology | Cosmic Microwave Background,Dense Linear Algebra, high I/O | Dense Matrix |
| GTC | Magnetic Fusion | Particle in Cell,Vlasov-Poisson | Particle/Grid |
| PARATEC | Material Science | Density Functional Theory, Kohn Shan, 3DFFT | Fourier/Grid |
| FVCAM | Climate Modeling | AGCM,Finite Volume, Navier-Stokes, FFT | Grid |

- Examining candidate ultra-scale applications with abundant data parallelism
- Codes designed for superscalar architectures, required vectorization effort
- ES use requires minimum vectorization and parallelization hurdles_x0004_

## Astrophysics: CACTUS

- Numerical solution of Einstein’s equations from theory of general relativity
- Among most complex in physics: set coupled nonlinear hyperbolic & elliptic systems w/ thousands of terms
- CACTUS evolves these equations to simulate high gravitational fluxes, such as collision of two black holes
- Evolves PDE’s on regular grid using finite differences
- Gravitational Wave Astronomy: exciting new field about to be born, w/ new fundamental insight into universe

- Visualization of grazing collision of two black holes

- Developed at Max Planck Institute, vectorized by John Shalf LBNL

**Notes:**
One of the most challenging problems is astrophysics is the numerical solution of Einstein’s equation from the theory of general relatively..
CACTUS code developed to evolve this equations stably in 3D to simulate astrophysical phenomena with high gravitational fluxes, such as the collision of two black holes (shown) and the gravitational waves that radiate from that event.

## CACTUS: Performance

- SX8 attains highest per-processor performance ever attained for Cactus:  58X Power3!
- Vector performance related to x-dim (vector length)
- Excellent scaling on ES/SX8 using fixed data size per proc (weak scaling)
- Opens possibility of computations at unprecedented scale
- X1 surprisingly poor (~6x slower than  SX8) - low ratio scalar:vector
- Unvectorized boundary, required 15% of runtime on SX8 and 30+% on X1
- < 5% for the scalar version: unvectorized code can quickly dominate cost
- To date unable to successfully execute code on X1E
- Poor superscalar performance despite high computational intensity
- Register spilling due to large number of loop variables
- Prefetch engines inhibited due to multi-layer ghost zones calculations

| ProblemSize | P | Power 3
Seaborg |  | Itanium2Thunder |  | OpteronJaquard |  | X1Phoenix |  | SX6*ES |  | SX8HLRS |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  |  | GFs/P | %pk | GFs/P | %pk | GFs/P | %pk | GFs/P | %pk | GFs/P | %pk | GFs/P | %pk |
| 250x80x80perprocessor | 16 | 0.10 | 6% | 0.58 | 10% | 0.82 | 19% | 0.81 | 6% | 2.8 | 35% | 4.3 | 27% |
|  | 64 | 0.08 | 6% | 0.56 | 10% | 0.72 | 21% | 0.72 | 6% | 2.7 | 34% | 4.0 | 25% |
|  | 256 | 0.07 | 5% | 0.55 | 10% | 0.68 | 16% | 0.68 | 5% | 2.7 | 34% | 3.9 | 24% |

## Plasma Physics: LBMHD

- LBMHD uses a Lattice Boltzmann method to model magneto-hydrodynamics (MHD)
- Performs 2D/3D simulation of high temperature plasma
- Evolves from initial conditions and decaying to form current sheets
- Spatial grid coupled to octagonal streaming lattice
- Block distributed over processor grid
- Main computational components:
- Collision, Stream, Interpolation
- Vectorization: loop interchange, unrolling
- X1 compiler loop reordered automatically

- Ported by Jonathan Carter, developed by George Vahala’s group College of William & Mary

- Evolution of vorticity into turbulent structures

**Notes:**
Lattice boltzman methods popular approach for modeling the physics of fluid flows. Recently groups have applied the LBM to the problem of MHD.   
Mention cache oblivious work in progress

## LBMHD-3D: Performance

| Grid Size | P | Power3 Seaborg |  | Itanium2 Thunder |  | OpteronJacquard  |  | X1Phoenix |  | X1EPhoenix  |  | SX6 ES |  | SX8HLRS |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  |  | GFs/P | %pk | GFs/P | %pk | GFs/P | %pk | GFs/P | %pk | GFs/P | %pk | GFs/P | %pk | GFs/P | %pk |
| 2563 | 16 | 0.14 | 9% | 0.26 | 5% | 0.70 | 16% | 5.2 | 41% | 6.6 | 37% | 5.5 | 69% | 7.9 | 49% |
| 5123 | 64 | 0.15 | 9% | 0.35 | 6% | 0.68 | 15% | 5.2 | 41% | 5.8 | 32% | 5.3 | 66% | 8.1 | 51% |
| 10243 | 256 | 0.14 | 9% | 0.32 | 6% | 0.60 | 14% | 5.2 | 41% | 6.0 | 33% | 5.5 | 68% | 9.6 | 60% |
| 20483 | 512 | 0.14 | 9% | 0.35 | 6% | 0.59 | 13% |  |  | 5.8 | 32% | 5.2 | 65% |  |  |

- Not unusual to see vector achieve > 40% peak while superscalar architectures achieve ~ 10%
- There exists plenty of computation, however large working set causes register spilling scalars
- Examining cache-oblivious algorithms for cache-based systems
- Unlike superscalar approach, large vector register sets hide latency
- Opteron shows impressive superscalar performance, 2X speed vs. Itanium2
- Opteron has >2x STREAM BW, and Itanium2 cannot store FP in L1 cache
- ES sustains 68% of peak up to 4800 processors: 26TFlops - SC05 Gorden Bell Finalist	 The highest performance ever attained for this code by far
- SX8 shows highest raw performance, but ES shows highest efficiency
- SX8: Commodity DDR2-SDRAM vs. ES: high performance custom FPLRAM
- X1E achieved same performance as X1 using original code version
- By turning off caching resulted in about 10% improvement over X1

## Analysis of Bank Conflicts

- Comparison of bank conflicts shows that ES custom memory has clear advantages

| Code | P | Grid Size | Bank Conflict |  |
| --- | --- | --- | --- | --- |
|  |  |  | ES | SX8 |
| MHD2D | 64 | 81922 | 0.3% | 16.6% |
| MHD2D | 256 | 81922 | 0.3% | 10.7% |
| MHD3D | 64 | 2563 | 0.01% | 1.1% |
| MHD3D | 256 | 5123 | 0.01% | 1.2% |

## Magnetic Fusion: GTC

- Gyrokinetic Toroidal Code: transport of thermal energy (plasma microturbulence)
- Goal magnetic fusion is burning plasma power plant producing cleaner energy
- GTC solves 3D gyroaveraged kinetic system w/ _particle-in-cell approach_ (PIC)
- PIC scales *N *instead of *N**2* – particles interact w/ electromagnetic field on grid
- Allows solving equation of particle motion with ODEs (instead of nonlinear PDEs)
- Vectorization inhibited: multiple particles may attempt to concurrently update same point in charge deposition

- Whole volume and cross section of electrostatic potential field, showing elongated turbulence eddies

- Developed at PPPL, vectorized/optimized by Stephane Ethier and ORNL/Cray/NEC** **

**Notes:**
GTC is 3D PIC code used for studying turbulent transport in magnetic fusion plasma
The simulation geometry is that torus - the natural configuration of tokamak fusion devices
As the charged particles forming the plasma move within the magnetic field, they collectively create their own self-consistent electrostatic field, then quickly becomes turbulent under driving temperature and density gradients 
Performance of such device is determine by rate at which energy is transported out of hot core to the colder edge of plasma
PIC solves 5D gyro-averaged kinetic equation coupled to the Poisson equation.
PIC: Instead of interacting with one another, particles interact with electromagnetic field described on a grid

## GTC Particle Decomposition

**Vectorization**:
- Particle charge deposited amongst nearest grid points
- Several particles can contribute to same grid point preventing vectorization
- Solution: VLEN copies of charge deposition array with reduction after main loop
- Greatly increases memory footprint (8x)
- GTC originally optimized for superscalar SMPs using MPI/OpenMP
- However: OpenMP severely increase memory for vector implementation
- Vectorization and thread-level parallelism compete w/ each other
- Previous vector experiments limited to only 64-way MPI parallelism
- 64 is optimal domains for 1D toroidal (independent of # particles)
**Updated GTC algorithm**
- Introduces a third level of parallelism:
- Algorithm splits particles between several processors (within 1D domain)
- Allows increase concurrency and number of studied particles
- Larger particle simulations allow increase resolution studies
- Particles not subject to Courant condition (same timestep)
- Allows multiple species calculations

## GTC: Performance

- New decomposition algorithm efficiently utilizes high P (as opposed to 64 on ES)
- Breakthrough of Tflop barrier on ES for important SciDAC code
- 7.2 Tflop/s on 4096 ES processors
- SX8 highest raw performance (of any architecture)
- Opens possibility of new set of high-phase space-resolution simulations
- Scalar architectures suffer from low CI, irregular data access, and register spilling
- Opteron/IB is 50% faster than Itanium2/Quadrics and only 1/2 speed of X1
- Opteron: on-chip memory controller and caching of FP L1 data
- X1 suffers from overhead of scalar code portions
- Original (unmodified) X1 version performed 12% *slower* on X1E
- Additional optimizations increased performance by 50%!
- Recent ORNL work increases performance additional 75%
- SciDAC code, HPCS benchmark

| P | Part/Cell | Power3 Seaborg |  | Itanium2 Thunder |  | OpteronJacquard  |  | X1Phoenix |  | X1EPhoenix  |  | SX6 ES |  | SX8HLRS |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  |  | GFs/P | %pk | GFs/P | %pk | GFs/P | %pk | GFs/P | %pk | GFs/P | %pk | GFs/P | %pk | GFs/P | %pk |
| 128 | 200 | 0.14 | 9% | 0.39 | 7% | 0.59 | 13% | 1.2 | 9% | 1.7 | 10% | 1.9 | 23% | 2.3 | 14% |
| 256 | 400 | 0.14 | 9% | 0.39 | 7% | 0.57 | 13% | 1.2 | 9% | 1.7 | 10% | 1.8 | 22% | 2.3 | 15% |
| 512 | 800 | 0.14 | 9% | 0.38 | 7% | 0.51 | 12% |  |  | 1.7 | 9% | 1.8 | 22% |  |  |
| 1024 | 1600 | 0.14 | 9% | 0.37 | 7% |  |  |  |  |  |  | 1.8 | 22% |  |  |

**Notes:**
Recent X1E optimization at ORNL involving reversal order of main array dimensions, resulted in additional 75% improvement to 3 Gflops or 17% of peak (similar optimization has not been attempted on ES/SX8)

## Material Science: PARATEC

- First-principles quantum mechanical total energy calc using pseudopotentials & plane wave basis set
- Density Functional Theory to calc structure & electronic properties of new materials
*DFT calc are one of the largest consumers of supercomputer cycles in the world*
- 33% 3D FFT,  33% BLAS3,  33% Hand coded F90
- Part of calculation in real space other in Fourier space
- Uses specialized 3D FFT to transform wavefunction

- Conduction band minimum electron state forCadmium Selenide (CdSe) quantum dot

- Developed by Andrew Canning with Louie and Cohen’s groups (LBNL, UCB)

- Global transpose Fourier to real space Multiple
- FFTs: reduce communication latency
- Vectorize across (not within) FFTs
- Custom FFT: only nonzeros sent

**Notes:**
DFT has become standard technique in Mat Sci to calculate the structural and electronic properties of new materials with a full QUANTUM MECHANICAL treatment of the electrons.
CdSw quantum dots can be used as electronic dye tags by attaching them to organic molecules - they represent systems with important technological applications.
488 atom system, is the largest physical system ever run with Paratec

## PARATEC: Performance

- All architectures generally perform well due to computational intensity of code (BLAS3, FFT)
- ES achieves highest overall performance to date: 5.5Tflop/s on 2048 procs
- Main ES/SX8 advantage for this code is fast interconnect
- Allows never before possible, high resolution simulations
- CdSE Q-dot: Largest cell-size atomistic experiment ever run using PARATEC
- Important uses: can be attached to organic molecules as electronic dye tags
- SX8 achieves highest per-processor performance (>2x X1E)
- X1/X1E shows lowest % of peak
- Non-vectorizable code much more expensive on X1/X1E (32:1)
- Lower bisection bandwidth to computational ratio (4D-hypercube)
- X1 performance is comparable to Itanium2
- Itanium2 outperforms Opteron (unlike LBMHD/GTC) because
- PARATEC less sensitive to memory access issues (high computational intensity)
- Opteron lacks FMA unit
- Quadrics shows better scaling of all-to-all at large concurrencies

| Problem | P | Power3 Seaborg |  | Itanium2 Thunder |  | OpteronJacquard  |  | X1Phoenix |  | X1EPhoenix  |  | SX6 ES |  | SX8HLRS |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  |  | GFs/P | %pk | GFs/P | %pk | GFs/P | %pk | GFs/P | %pk | GF/s/P | %pk | GFs/P | %pk | GFs/P | %pk |
| 488 AtomCdSeQuantumDot | 128 | 0.93 | 62% | 2.8 | 51% |  |  | 3.2 | 25% | 3.8 | 21% | 5.1 | 64% | 7.5 | 49% |
|  | 256 | 0.85 | 67% | 2.6 | 47% | 2.0 | 45% | 3.0 | 24% | 3.3 | 18% | 5.0 | 62% | 6.8 | 43% |
|  | 512 | 0.73 | 49% | 2.4 | 44% | 1.0 | 22% |  |  | 2.2 | 12% | 4.4 | 55% |  |  |
|  | 1024 | 0.60 | 40% | 1.8 | 32% |  |  |  |  |  |  | 3.6 | 46% |  |  |

## Atmospheric Modeling: FVCAM

- CAM3.1: Atmospheric component of CCSM3.0
- Our focus is the Finite Volume (FV) approach
- AGCM: consists of physics (PS) and dynamical core (DC)
- DC approximates dynamics of atmosphere
- PS: calculates source terms to equations of motion:
- Turbulence, radiative transfer, clouds, boundary, etc
- DC default: Eulerian *spectral* transform - maps onto sphere
- Allows only 1D decomposition in latitude
- FVCAM grid is rectangular (longitude, latitude, level)
- Allows 2D decomp (latitude, level) in dynamics phase
- Singularity at pole prevents longitude decomposition
- Dynamical eqns bounded by Lagrangian surfaces
- Requires remapping to Eulerian reference frame
- Hybrid (MPI/OpenMP) programming
- MPI tasks limited by # of latitude lines
- Increase potential parallelism and improves S2V ratio
- Is not always effective on some platforms

- Experiments/vectorization Art Mirin, Dave Parks, Michael Wehner, Pat Worley

- Simulated Class IV hurricane at 0.5.  This storm was produced solely through the chaos of the atmospheric model.  It is one of the many events produced by FVCAM at resolution of 0.5.

## FVCAM Decomposition and Vectorization

**Processor communication topology** and volume for  1D and  2D FVCAM
- Generated by IPM profiling tool - used to understand interconnect requirements
- 1D approach straightforward nearest neighbor communication
- 2D communication bulk is nearest neighbor - however:
- Complex pattern due to vertical decomp and transposition during remapping
- Total volume in 2D remap is reduced due to improved surface/volume ratio
**Vectorization** - 5 routines, about 1000 lines
- Move latitude calculation to inner loops to maximize parallelism
- Reduce number of branches, performing logical tests in advance (indirect indexing)
- Vectorize across (not within) FFT’s for Polar filters
- However, higher concurrency for fixed size problem limit performance of vectorized FFTs

## FVCAM3.1: Performance

- FVCAM 2D decomposition allows effective use of >2X as many processors
- Increasing vertical discretizations (1,4,7) allows higher concurrencies
- Results showing high resolution vector performance 361x576x26 (0.5 x 0.625)
- X1E achieves speedup of over 4500 on P=672 - highest ever achieved
- To date only results from 1D decomp available on SX-8
- Power3 limited to speedup of 600 regardless of concurrency
- Factor of at least 1000x necessary for simulation to be tractable
- Raw speed X1E: 1.14X X1, 1.4X ES, 3.7X Thunder, 13X Seaborg
- At high concurrencies (P= 672) all platforms achieve low % peak (< 7%)
- ES achieves highest sustained performance (over 10% at P=256)
- Vectors suffer from short vector length of fixed problem size, esp FFTs
- Superscalars generally achieve lower efficiencies/performance than vectors
- Finer resolutions requires increased number of _more powerful_ processors

## Cosmology: MADCAP

- Anisotropy Dataset Computational Analysis Package
- Optimal general algorithm for extracting key cosmological data from Cosmic Microwave Background Radiation (CMB)
- Anisotropies in the CMB contains early history of the Universe
- Recasts problem in dense linear algebra: ScaLAPACK
- Out of core calculation: holds ~3 of the 50 matrices in memory

- Temperature anisotropies in CMB (Boomerang)

- Developed by Julian Borrill, LBNL

**Notes:**

Picture shows induced current and charge density in crystalized glcine.The simulations used to better undersand nuclear magnetic resonance.

## MADCAP: Performance

- Overall performance can be surprising low, for dense linear algebra code
- I/O takes a heavy toll on Phoenix and Columbia: I/O optimization currently in progress
- NERSC Power3 shows best system balance with respect to I/O
- ES lacks high-performance parallel I/O - code rewritten to utilize local disks
- Preliminary SX8 experiments unsuccessful due to poor I/O performance

| NumberPixels | P | NERSC (Power3) |  | Columbia (Itan2) |  | Phoenix (X1) |  | ES (SX6*) |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  |  | Gflops/P | %peak | Gflops/P | %peak | Gflops/P | %peak | Gflops/P | %peak |
| 10K | 64 | 0.73 | 49% | 1.2 | 20% | 2.2 | 17% | 2.9 | 37% |
| 20K | 256 | 0.76 | 51% | 1.1 | 19% | 0.6 | 5% | 4.0 | 50% |
| 40K | 1024 | 0.75 | 50% |  |  |  |  | 4.6 | 58% |

## Performance Overview

- Tremendous potential of vector systems: SX8 achieves unprecedented raw performance
- Demonstrated importance of vector platforms: allows solutions at unprecedented scale
- ES achieves highest efficiency on almost all test applications
- LBMHD-3D achieves 26 TF/s using 4800 procs, GTC achieves 7.2 TF/s on 4096 ES procs, PARATEC achieves 5.5 TF/s on 2048 procs
- Porting several of our apps to SX8 in progress:Madbench (slow I/O), FVCAM 2D decomp, Hyperclaw AMR (templated C++)
- Investigating how to incorporate vector-like facilities into modern micros
- Many methods unexamined (more difficult to vectorize): Implicit, Multigrid, AMR, Unstructured (Adaptive), Iterative and Direct Sparse Solvers, N-body, Fast Multiple
- One size does not fit all: need a range of architectural designs to best fit algorithms
- Heterogeneous supercomputing designs are on the horizon
- Moving on to latest generation of HEC platforms: Power5, BG/L, XT3

**Speedup vs. ES**

**Notes:**
With the emergence of the next generation of supercomputing systems - it is clear that no single architecture is best suited for the full scope of application domains - it is therefore critical to understand which computational areas are best suited for a given architectural paradigm, to maximize scientific productivity.