Domain-Specific Languages for
Heterogeneous GPU Computing
Pat Hanrahan
Pervasive Parallelism Laboratory- NVIDIA is a Founding Member Stanford University
NVIDIA Technology Conference
October 2, 2009

Jen-Hsun’s Argument for Co-Processing
Power efficiency
 CPUs more efficient for sequential workloads
 GPUs more efficient for data-parallel workloads

Amdahl’s Law
 Applications have a mixture of sequential andparallel code Parallelism may be limited by sequential

performance

Therefore, Optimal platform involve both fast sequential coresplus fast parallel cores

Page 1ZuneHD uses Tegra

Tegra APX2600: CPU+GPU+HD+…
CPU: ARM11 MPCore
GPU: ULP
HD: 720p decoder…

Page 2

Platforms even more Heterogeneous
Cluster (e.g. 100K node cluster) Distributed memory
 MPI for communication

Multi-core SMP (e.g. 32 core, 4-socket systems)
 Shared memory, transaction memory Threads/locks, OpenMP

Many-core GPU (e.g. Fermi)
 Separate GPU memory
 SIMT programming model (CUDA)+ combinations of these systems (accelerators)

Is is Possible to Write One ProgramandRun it on all these Machines?

Page 3

HYPOTHESIS: YES, BUT NEED
Domain SpecificProgramming Environments
andLanguages

Domain-Specific Languages
Widely used in many application areas
 matlab / R
 SQL / Microsoft’s LINQ / map-reduce OpenGL/D3D and Cg/HLSL …

DSLs are a hot topic now Programming language community (C#, Scala) Web programming environments (Ruby)

Page 4

Three Advantages
Productivity
 Separate domain expertise (computational science)from computer science expertisePortability
 Run on wide range of platforms

Performance Super-optimize using a combination of domain

knowledge and platform knowledge

RenderMan Shading Languagesurface corrode(float Ks=0.4, Ka=0.1, rough=0.25){
float i, freq=1, turb=0;// compute fractal texture
for( i=0; i<6; i++ ) {turb+=1/freq*noise(freq*P);
freq*=2;}// perturb surface
P -= turb * normalize(N);N = faceforward(normalize(calculatenormal(P)));// compute reflection and final color
Ci = Cs*(Ka*ambient()+Ks*specular(N,I,rough));}

Page 5Graphics Libraries
glMatrixMode(GL_PROJECTION);
glPerspective(45.0);
for( ;; ) {glBegin(TRIANGLES);
glVertex(…);
glVertex(…);…
glEnd();}
glSwapBuffers();

OpenGL “Grammar”<Scene> = <BeginFrame> <Camera> <World> <EndFrame><Camera> = glMatrixMode(GL_PROJECTION) <View><View> = glPerspective | glOrtho<World> = <Objects>*<Object> = <Transforms>* <Geometry><Transforms> = glTranslatef | glRotatef | …
<Geometry> = glBegin <Vertices> glEnd<Vertices> = [glColor] [glNormal] glVertex

Page 6

Productive, Performant, and Portable
Easy to write OpenGL programs (CG101)
Runs at extremely high speed
 Graphics drivers do amazing things Graphics chips do amazing things

Runs across very different architectures
 Wide range of GPUs on the market
 Backward and forward compatible

Three Four Advantages
Productivity
Portability
Performance
Encourage innovation
Allow NVIDIA to radically optimize implementation
Reduce cost of porting if the low-level programming
model changes

Page 7

DSLs in Other Parallel Applications?
Future graphics pipelines with rendering and physics
Statistics/machine learning and data analysis (beyond R)
Physical simulation
Computer vision and imaging
Brain simulation
Autonomous vehicles…

Liszt
Z. Devito, M. Medina, M. Barrientos,
J. Alonso, E. Darve, F. Ham
P. Hanrahan“…the most technically advanced and perhaps
greatest pianist of all time…”
1811-1886

Page 8

Stanford DOE PSAAP Center
Characterize the operability limits of a hypersonicpropulsion system using predictive computations.
Primary focus is the unstart phenomena triggered bythermal choking in a hydrogen-fueled scramjet.

Joe	  • State-­‐of-­‐the-­‐art	  unstructured	  RANS	  solver	  – Main	  tool	  for	  system-­‐level	  simula:on	  Combustion
Turbulence
ThermalFuel injection
Transition

Turbulence

Comparison	  of	  Joe	  simula0on	  result	  to	  experiment:	  DLR	  ground	  test	  based	  on	  atmospheric	  condi0ons	  at	  27	  km	  Page 9Typical Joe C Code Fragmentfor( int ifa = zone->ifa_f; ifa <= zone->ifa-1; ifa++ ) {}int naf_f = noofa_i[ifa];
int nof_l = noofa_i[ifa+1] – 1;
int nnof = nof_l – nof_f + 1;
double x_fa[3], u_fa_bc[3];for( int i = 0; i< 3; i++ ) {x_fa[i] = 0.0;
u_fa_bc[i] = 0.0;}
double phi_fa_bc = 0.0;
for( int nof = nof_f; nof <= nof_l; nof++ ) {int ino = noofa_v[nof];
int ibc = bc_flag[ino];
no_flag[ino] = nof – nof_f;for( int i = 0; i < 3; i++ ) {x_fa[i] += x_no[ino][i];u_fa_bc[i] += ubc->data[ibc][i];}phi_fa_bc += phibc->data[ibc];}

Joe in Liszt (Scala Version)
val pos = new Field[Vertex,double3]
val A = new SparseMatrix[Vertex,Vertex]for (c <- mesh.cells) {val center = avg pos of c.verticesfor (f <- c faces) {val face_dx = avg pos of f.vertices – centerfor (e <- f edgesWith c CounterClockwise) {val v0 = e.tailval v1 = e.headval v0_dx = pos(v0) – centerval v1_dx = pos(v1) – centerval face_normal = v0_dx cross v1_dx// calculate flux for face …
A(v0,v1) += …
A(v1,v0) -= …

Page 10

Built-in Features
Objects
 Short vectors, vectors, dense and sparse matrices
 Mesh, cells, faces, edges, vertices
 Fields
 Sets (unordered) and lists (ordered)

Solvers
 Sparse matrix solvers (e.g. hypre, PETSc, …)Parallel iteration
 Map, reduce, …

Ghost Cell

Owned Cell

Partitioning (ParMETIS)~100K Cells/Core

Page 11

Domain-Specific Optimizations
1.

Automatic partitioning

2.

Automate maintenance of distributed ghost cells

3.

Optimize mesh representation

4.

Optimize storage of fields

5.

…

Liszt uses program analysis
to perform these optimizations

Page 12

Program Analysis of Neighborhoods
Field<Cell, double> rho;
Field<Face, double> rhoOutside;forall (Face f in mesh.faces()) {rhoOutside(f) =+calc_flux( f,rho(f.outside() ))calc_flux( f,rho(f.inside()));}

Program Analysis of Neighborhoods
Field<Cell, double> rho;

Node 0Field<Face, double> rhoOutside;

Owned
Cells

forall (Face f in mesh.faces()) {rhoOutside(f) =+Ghost Cells

calc_flux( f,rho(f.outside() ))calc_flux( f,rho(f.inside()}

Page 13));

Program Analysis of Neighborhoods
Field<Cell,double> rho;
Field<Face,double> rhoOutside;forall (Face f in mesh.faces()) {Cell outside = f.outside();rhoOutside(f) = 0.5 * rho(c));
forall (Cell c in outside.cells()) {rhoOutside(f) += .25 * rho(c);}}Program Analysis of Neighborhoods
Field<Cell,double> rho;

Node 0Field<Face,double> rhoOutside;forall (Face f in mesh.faces()) {Cell outside = f.outside();rhoOutside(f) = 0.5 * rho(c));
forall (Cell c in outside.cells()) {rhoOutside(f) += .25 * rho(c);}}Ghost Cells

Page 14OwnedCells

Domain Decomposition / Ghost Cells
Node 0Liszt creates a graph ofmesh adjacencies needed
to run the algorithm
Graph is handed to Parmetis
to determine optimal partition

Owned
Cells

Communication ofinformation in ghost cells is
also automatically handed

Choose Mesh Representation
1.

Optimize based on what topological relationships areused by the simulation code

2.

Optimize based on special structure of the input mesh.
Input Meshes
Unstructured

TriangularFaces

Page 15

Regular GridChoose Mesh Representation (1)// Finite volume method

Meshforall(Cell c in mesh.cells()) {// calculate gradientsc0 c1 c2
f0 f1}
forall(Face f in mesh.faces()) {double flux = //flux for facedata(f.inside())-= flux;c0 c1

data(f.outside()) += flux;}
forall(Cell c in mesh.cells()) {// update cell state

f0c1

f1

c2c0}

Choose Mesh Representation (1)forall (Cell c in mesh.cells()) {double3 center = avg([pos(v) for v in c.vertices()]);
Meshc0 c1 c2
forall (Face f in c.faces()) {double3 face_dx = avg([pos(v) for v in f.vertices()]) – center;
forall (Edge e in f.edgesCCW(c)) {vertex v0 = e.tail();vertex v1 = e.head();v1 v2 f0 f1
double3 v0_dx = pos(v0) – center;v3 v5 f3 f5
double3 v1_dx = pos(v1) – center;double3 face_normal = v0_dx.cross(v1_dx);// calculate flux for face …
A[v0][v1] += …
f3
v3
v1
e0 e1 e2
A[v1][v0] -= …
e1v1

v5

v1 v2 v5e0v2 v1

f0f0
e2

c1c2c0v2

v2

Page 16f1Choose Mesh Representation (2)forall( Face f in mesh.faces() ) {forall( Face f2 in f.faces() ) {field(f)++;}}
Input Meshes
Unstructured

Triangular
Faces

Mesh 1 – Unstructured
forall( Face f in mesh.faces() ) {forall( Face f2 in f.faces() ) {field(f)++;}}

Page 17

Regular Grid

Unstructured
mixture ofquads andtriangles

Unstructured
mixture ofquads andtriangles

Mesh 1 – Unstructured
forall( Face f in mesh.faces() ) {forall( Face f2 in f.faces() ) {field(f)++;}}

e0e1v0f0f1

Winged Edge Data Structure
struct Face {Edge *e;int id;};

struct Edge {Vert *v0, *v1;
Face *f0, *f1;
Edge *e0, *e1,*e2, *e3;
};

v1
e2e3

Unstructured
mixture ofquads andtriangles

Mesh 1 – Unstructuredfor(int i = 0; i < nfaces; i++) {Face * f = faces[i];Edge * e = f->e;do {Face * f2 = faceAroundEdge(e,f);field[f2->id]++;} while( (e = edgeCCWAroundFace(e,f))!= f->e );}struct Face {Edge *e;
int id;};struct Edge {Vert *v0, *v1;
Face *f0, *f1;
Edge *e0, *e1,*e2, *e3;};

e0e1v0f0

f1
v1e2

Page 18e3

Mesh 2 – Triangular Faces

Triangular Faces

forall(Face f in mesh.faces() ) {forall(Face f2 in f.faces() ) {field(f)++;}}

Mesh 2 – Triangular Faces

Triangular Faces

forall(Face f in mesh.faces() ) {forall(Face f2 in f.faces() ) {field(f)++;}}

v1
Triangle mesh representation
struct Face {Edge *e0,*e1,*e2;
Face *f0,*f1,*f2;
Vert *v0,*v1,*v2;
int id;};f1
e1e0f0fv2e2
f2

Page 19v0Mesh 2 – Triangular Faces

Triangular Faces

for(int i = 0; i < nfaces; i++) {Face * f = faces[i];field[f->f0->id]++;field[f->f1->id]++;field[f->f2->id]++;}v1
Triangle mesh representation
struct Face {Edge *e0,*e1,*e2;
Face *f0,*f1,*f2;
Vert *v0,*v1,*v2;
int id;};f1
e1e0f0fv2e2v0f2

Mesh 3 – Regular Gridforall(Face f in mesh.faces() ) {forall(Face f2 in f.faces() ) {field(f)++;}}

Page 20

Regular Grid

Mesh 3 – Regular Grid

Regular Grid

forall(Face f in mesh.faces() ) {forall(Face f2 in f.faces() ) {field(f)++;}}

Simple mesh data structure
int nfaces = faces_x * faces_y;
int field[faces_x][faces_y];

Mesh 3 – Regular Gridfor(int x = 1; x < faces_x - 1; x++) {for(int y = 1; y < faces_y - 1; y++) {field[x+1][y]++;field[x-1][y]++;field[x][y+1]++;
field[x][y-1]++;}}• Can vectorize the loop body to handle 4 faces at once
• Can do intelligent blocking of the matrix operations
Simple mesh data structure
int nfaces = faces_x * faces_y;
int field[faces_x][faces_y];

Page 21

Regular Grid

Array of Structs vs. Struct of Arrays
K1	

K2	

K3	

K4	

K5	

enthalpy	

gamma	

gamma local_dt	

pressure	rho	

Lambda- local_dt	

rhoE	

muLam	OverCp	 pressure	 rhou	

temperature	

rhou	

rho	

rhoE	

muLam	rho	

RoM	

temperature	

3/3
gamma

3/3

pressure

3/3

1/3

K7	

rhou	

rho	u	enthalpy	

gamma	

LambdaOverCp	

muLam	

pressure	rho	u	- Fields accessed per kernel- Prob[B access|A access]- Determine optimal field layouts

rho
3/5

3/3

K6	

1/5

3/5
1/3
muLam
1/3

Target Application Areas
Future rendering pipelines (Gramps)
Data-parallel programming (Kore)
Statistics/machine learning - beyond R (OptiML)
Physical simulation (Liszt)…
Computer vision and imaging
Brain simulation
Autonomous vehicles…

Page 22

How to Create
Domain-Specific Languages“Little” Languages (UNIX)e.g. sed, awk, sh, Renderman, matlab, R, …
Typically built using compiler generators: yacc and lex
Disadvantages:- Continual requests for more generality- Proliferation of syntax: “$x”
- Cannot use multiple DSLs in the same application

Page 23

Domain-Specific Embedded Language
Examples: Lisp, C++ templates, C#, Haskell, Ruby, Scala
Advantages:+ Consistent base syntax+ Multiple DSLs may interoperate

Rake – Make in Ruby
SRC = FileList['*.c']
OBJ = SRC.ext('o')
task :default => ["hello”]
rule '.o' => '.c' do |t|
sh "cc -c -o #{t.name} #{t.source}"end
file "hello" => OBJ dosh "cc -o hello #{OBJ}"
end# File dependencies go here ...
file 'main.o' => ['main.c', 'greet.h']file 'greet.o' => ['greet.c']Page 24

Rake – Make in Ruby
Example Rakefile
Features of Ruby for DSELs
 Overloaded operators
 Simple (Smalltalk-like) syntax: a b c -> a.b(c) Blocks: Lambdas and closures Dynamic typing

LINQ (Language-Integrated Query)var personsNotInSeattle =from p in person
where p.Address != "Seattle”
orderby p.FirstName
select p;

Page 25

LINQ Driving Feature Set for C# 3.0
Language extensions
 Implicitly typed variables
 Lambdas
 Anonymous classes
 Extensions (implicit type wrapper)

Enables parallel data analysis
 SQL engines
 PLINQ: SMP DryadLINQ: Clusters

Domain-Specific Embedded Language
Approach 1: Objects w/ operators+ Access to most embedding language features+ Compiled code
- No general system for program analysis or programtransformations

Page 26

Domain-Specific Embedded Language
Approach 2: Build expression tree and interpret+ Allows program analysis and transformation+ Extend semantics of embedding language- Hides some features of the embedding language(e.g. types and type checking)- Slow because it is interpreted

Integrated DSEL Environment
1.

Flexible syntax and grammar
 The flexibility of lisp macros in a modern language
 Convert to simple abstract syntax tree2.

Program analysis
 Library supporting domain-specific analysis

3.

Program transformation
 Domain specific rewrite rules

4.

Code generation
 Generate code for different parallel platforms

5.

Language workbench / IDE

Page 27

DELITE Approach to DSELs
Support DSELs by embedding in Scala

Build dataflow graph by deferring method execution
Use static analysis to complement dynamicoptimizations

Apply a variety of generic and domain transformation
to dataflow graphIntelligently schedule work to maximize locality withdomain knowledge

DSELs in Scala
Scala combines both the strengths of FP and OOPEncourages use of best approach for a particular problem

Designed to embed DSLs
Concise syntax, implicit type conversions, …

Functional programming
Higher-order functions allow parallel control structures (fork,…)Declarative programming style ⇒ exposes parallelism

Object-oriented programming


Familiar programming model

Allows mutable data structures

Scala compiles to Java bytecode

JVM platform is mature and high performanceAllows the use of Java libraries

Page 28

PPL Overview
ScientificEngineering

Applications

DomainSpecific
Languages

Rendering

VirtualWorlds

Physics

Personal
Robotics

ScriptingData
informatics

Probabilistic
Resoning

MachineLearning

Parallel Object Language
DSL
Infrastructure

Common Parallel Runtime
Explicit / Static

Heterogeneous
HardwareOOO Cores
Programmable
Hierarchies

Implicit / Dynamic

SIMD Cores
Scalable
Coherence

Threaded Cores
Isolation &Atomicity

Pervasive
Monitoring

Acknowledgements
Bill Dally

Zach Devito

Alex Aiken

Montse Medina

Kunle Olukotun

Michael Barrientos

Eric Darve

Nathan Bronson

Juan Alonso

Hassan Chafi

Frank Ham
Vijay Pande
Arthur Whitney

Page 29

Thank you

Page 30