Visit additional Tabor Communication Publications
October 08, 2008
In the not-too-distant past, ENIAC was programmed with switches and a plugboard. Stored program computers soon followed that allowed one to write a program, load it into the computer memory, and run it. Initially, those programs had to be written in or manually translated into binary machine code, but soon assembly languages and assemblers were developed to simplify the process.
Soon followed operating systems, multiprogramming, and the concept of an application binary interface (ABI). The ABI defines the interface between an application and the operating system, libraries and other components. One aspect of an ABI is to define a calling convention, including how arguments are passed to a function and where the return value can be retrieved. For instance, the x64 ABI defines that the first six integer or pointer arguments are passed in registers (%rdi, %rsi, %rdx, %rcx, %r8, %r9); the first eight floating point arguments (single or double precision) are passed in SSE registers (%xmm0 to %xmm7); and any remaining arguments are pushed on the stack (in right-to-left order). This allows up to 14 arguments to be passed in registers, which surely captures most function calls.
But not all; WRF, the Weather Research and Forecasting Model, is used for both atmospheric research and numerical weather prediction. A version of WRF is included in the SPEC CPU2006 suite. One routine (copying from the WRF source code) is "a mixed phase ice microphysics scheme" (WSM5), with 49 arguments; it calls a subroutine (WSM52D) to handle the two-dimensional physics with 47 arguments (19 integers, 18 floating point scalars, 10 floating point arrays). Imagine writing the routine call by hand in assembly language; it takes over 100 instructions just to marshall the arguments and put them in the right place.
Instead, the computing community created higher level programming languages. While the first compiler (for the A-0 system) was more like what we would today call a loader, programming languages and compilers have progressed to where we use many higher level languages (C, Java, Fortran, others too many to enumerate) with a great increase in productivity. Much programming is done without using a textual language at all; for instance, a spreadsheet is a form of a program, and various visual programming interfaces exist. Now, the routine call in WRF with 47 arguments takes one Fortran statement, much easier to write and maintain than the corresponding assembly code:
CALL wsm52D(t, q(ims,kms,j), qci, qrs, &
w(ims,kms,j), den(ims,kms,j), &
p(ims,kms,j), delz(ims,kms,j), rain(ims,j), &
rainncv(ims,j),delt,g, cpd, cpv, rd, rv, t0c, &
ep1, ep2, qmin, &
XLS, XLV0, XLF0, den0, denr, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
If programming in binary is akin to using fingers and teeth, and assembly language is like using sticks and stone knives, think of higher level languages as the power tools of programming.
The earliest GPUs were hardware graphics accelerators to handle line drawing, area fill, image transfer, and so on, offloading the CPU. The adoption of standardized libraries such as OpenGL and Direct3D drove the development of hardware 3D graphics accelerators, particularly with programmable shading capability. Since 2000, the programmability of the graphics accelerator chips has improved to the point where they can be used for nongraphics applications. This has been called GPGPU (General Purpose computation on GPUs). Early GPGPU programming used the existing graphics libraries, such as OpenGL, mapping between computing concepts (array, loop, execute) and graphics concepts (texture, kernel, draw). This is truly heroic programming. It's like using a chain saw to carve blocks of ice: in the right hands, it can produce something beautiful, but one wrong mistake and all you have is ice cubes (or worse).
More recently, the programming research and development community have tried to come up with programming models that would map well onto GPUs and similar parallel computers, particularly stream programming, evidenced in several projects: StreamIt at MIT, Sh at Waterloo (which led to RapidMind), Brook at Stanford (which spun out briefly as PeakStream, and which AMD has adopted and extended as Brook+), and others.
The GPU programming model that has caught the most attention is NVIDIA's CUDA. The language is an extension to C; the software includes compiler, libraries, and many examples. There is a large user community, including a few dozen universities that use it in course work. Moreover, the software is free, though it (obviously) only targets NVIDIA GPUs.
Another programming model, very similar to CUDA and being sponsored by Apple and others, is OpenCL. The programming models are so similar that I'll only point out the differences here.
My last column discussed the GPU architectures and some of the problems facing programmers who want to compute on one. The first is to have an application with enough of the right type of parallelism to map onto the GPU. As the most parallelizable simple example, let's see what it would take to port a matrix multiplication to the GPU.
In Fortran, a matrix multiplication looks like a triply-nested loop:
do i = 1,n
do j = 1,m
do k = 1,p
a(i,j) = a(i,j) + b(i,k)*c(k,j)
In C, we have to decide whether to store the arrays linearized in a long vector, or whether to use a vector of pointers, or whether we have the degenerate case with fixed size arrays. Let's assume we use linearized arrays:
for( int i = 0; i < n; ++i )
for( int j = 0; j < m; ++j )
for( int k = 0; k < p; ++k )
a[i+n*j] += b[i+n*k] * c[k+p*j];
Matmul is a wonderful example to use when experimenting with loop optimizations, because it can be rewritten in so many ways. The three loops can be interchanged or reordered in six ways, strip-mined or tiled, parallelized and vectorized. To optimize for vector instructions, we want the i (stride-1) index innermost, to maximize the memory fetch/store bandwidth. For parallel multiprocessor or multicore execution, we want the j index outermost, so each processor or core is computing distinct columns of a. To optimize for cache memories, we want to tile all the loops, so the innermost loops compute a submatrix multiplication where the submatrices all fit in cache. An optimized, parallelized, vectorized matmul for a quad-core processor might look like:
jts = j tile size;
its = i tile size;
kts = k tile size;
parfor( int p = 0; p < 4; ++p ) /* parallel loop */
for( int jt = p; jt < m; jt += 4*jts )
for( int it = 0; it < n; it += its )
for( int kt = 0; kt < p; kt += kts )
for( int j = jt; j < min(m,jt+jts); ++j )
for( int k = 0; k < min(p,kt+kts); ++k )
for( int i = 0; i < min(n,it+its); ++i ) /* vector mode */
a[i+n*j] += b[i+n*k] * c[k+p*j];
So, even optimizing this for a modern parallel workstation or server takes significant work, knowledge of the memory hierarchy, and experimentation. In the past, programmers would have to do this all manually, though advanced compiler technology is now able to achieve this kind of optimization automatically.
But we want to compute the matmul on the GPU, using CUDA. Let's list the steps we must take in our program to get there.
cudaMalloc( &dev_a, n*m*sizeof(float) );
cudaMalloc( &dev_b, n*p*sizeof(float) );
cudaMalloc( &dev_c, p*m*sizeof(float) );
cudaMallocPitch( &dev_a, &pitch_a, n*sizeof(float), m );
cudaMallocPitch( &dev_b, &pitch_b, n*sizeof(float), p );
cudaMallocPitch( &dev_c, &pitch_c, p*sizeof(float), m );
cudaMemcpy2D( dev_b, pitch_b, b, n*sizeof(float), n*sizeof(float), p,
cudaMemcpy2D( dev_c, pitch_c, c, p*sizeof(float), p*sizeof(float), m,
__global__ void mmkernel( float* a, float* b, float* c,
int pitch_a, int pitch_b, int pitch_c,
int n, int m, int p )
int i = blockIdx.x*32+threadIdx.x;
int j = blockIdx.y;
float sum = 0.0;
for( int k = 0; k < p; ++k )
sum += b[i+pitch_b*k] * c[k+pitch_c*j];
a[i+pitch_a*j] = sum;
In CUDA, this is done with a few lines:
dim3 threads( 32 );
dim3 grid( n/32, m );
mmkernel<<< grid, threads >>>( dev_a, dev_b, dev_c,
pitch_a, pitch_b, pitch_c, n, m, p );
OpenCL is not so convenient. Since it is library-based, it can't depend on a compiler to simplify the steps. Instead, we will have to do each step explicitly, something approaching:
/* program is a prebuilt kernel program */
kernel = clCreateKernel( program, "mmkernel" );
grid = n;
grid = m;
threads = 32;
threads = 1;
/* context is the GPU compute context */
range = clCreateNDRangeContainer( context, 0, 2, grid, threads );
clSetKernelArg( kernel, 0, dev_a, sizeof(dev_a), NULL );
clSetKernelArg( kernel, 1, dev_b, sizeof(dev_b), NULL );
clSetKernelArg( kernel, 2, dev_c, sizeof(dev_c), NULL );
clSetKernelArg( kernel, 3, pitch_a, sizeof(pitch_a), NULL );
clSetKernelArg( kernel, 4, pitch_b, sizeof(pitch_b), NULL );
clSetKernelArg( kernel, 5, pitch_c, sizeof(pitch_c), NULL );
clSetKernelArg( kernel, 6, n, sizeof(n), NULL );
clSetKernelArg( kernel, 7, m, sizeof(m), NULL );
clSetKernelArg( kernel, 8, p, sizeof(p), NULL );
/* queue is a GPU work queue */
clExecuteKernel( queue, kernel, NULL, range, NULL, 0, NULL );
cudaMemcpy2D( a, n*sizeof(float), dev_a, pitch_a, n*sizeof(float), m,
cudaFree( dev_a );
cudaFree( dev_b );
cudaFree( dev_c );
CUDA (and, from reports, OpenCL) lets us program the GPU using a familiar language, C. However, a great deal of the programming is done using library calls, and the computation itself, the matrix multiplication, is divided into two parts: the kernel, on the GPU, and its invocation on the host. What was a multi-dimensional loop, which modern compilers are pretty good at optimizing, has turned into dozens of lines of code to manage memory, move data, and deal with the architectural specifics of the GPU.
I'm sure some (or most, or perhaps all) of the readers will at this point say "You shouldn't be programming matrix multiplication anyway; just call a library routine." That's true; in fact, NVIDIA provides versions of the BLAS routines, including SGEMM, which give very good performance. But, if it takes this much effort to get a matrix multiplication moved to the GPU, imagine how much effort it takes to move a real computation (say, 180,000 lines of WRF, or even just the 10,000 lines of microphysics). Here I had only three arrays with regular access patterns and full freedom to parallelize the loops. Can you begin to see the difficulties?
And I haven't even begun the real programming process, such as handling error returns from the runtime calls. Given the CUDA or OpenCL code, how portable will it be, how will I maintain it to keep it at the peak of efficiency?
As I mentioned above, the CUDA NVCC compiler simplifies some of the coding details. Compilers are good at bookkeeping, organizing details about memory addressing, alignment, and so on. OpenCL seems to be a step backwards, from a compiler-oriented solution aimed at raising the level of programming to a library-oriented solution aimed at giving low-level control to the programmer.
Low-level control isn't bad; that would be like saying assembly language is bad. It isn't bad, but it's only necessary for a very small bit of the programming we do today. If high level languages are the power tools of programming, we seem to be taking a step back to hand drills and saws. Many woodworkers prefer hand tools, and they can make beautiful furniture, but the cost is high and productivity is low. We need hand tools, but we'll be much more productive with better power tools.
In my next column, I'll look at the matrix multiplication kernel, exploring various ways to optimize for parallelism and memory bandwidth, and presenting performance. Whether you think programming the host side of a GPU program is hard, or just work, you'll be entertained, enlightened or frightened when you see what goes into the GPU side. I'm on record as saying that parallel programming isn't easy and never will be, but we can and should develop the tools and training to turn this Acceleration Nightmare into something about as scary as a Halloween Haunted House.
Jun 19, 2013 |
Supercomputer architectures have evolved considerably over the last 20 years, particularly in the number of processors that are linked together. One aspect of HPC architecture that hasn't changed is the MPI programming model.
Jun 18, 2013 |
The world's largest supercomputers, like Tianhe-2, are great at traditional, compute-intensive HPC workloads, such as simulating atomic decay or modeling tornados. But data-intensive applications--such as mining big data sets for connections--is a different sort of workload, and runs best on a different sort of computer.
Jun 18, 2013 |
Researchers are finding innovative uses for Gordon, the 285 teraflop supercomputer housed at the San Diego Supercomputer Center (SDSC) that has a unique Flash-based storage system. Since going online, researchers have put the incredibly fast I/O to use on a wide variety of workloads, ranging from chemistry to political science.
Jun 17, 2013 |
The advent of low-power mobile processors and cloud delivery models is changing the economics of computing. But just as an economy car is good at different things than a full size truck, an HPC workload still has certain computing demands that neither the fastest smartphone nor the most elastic cloud cluster can fulfill.
Jun 14, 2013 |
For all the progress we've made in IT over the last 50 years, there's one area of life that has steadfastly eluded the grasp of computers: understanding human language. Now, researchers at the Texas Advanced Computing Center (TACC) are utilizing a Hadoop cluster on its Longhorn supercomputer to move the state of the art of language processing a little bit further.
05/10/2013 | Cleversafe, Cray, DDN, NetApp, & Panasas | From Wall Street to Hollywood, drug discovery to homeland security, companies and organizations of all sizes and stripes are coming face to face with the challenges – and opportunities – afforded by Big Data. Before anyone can utilize these extraordinary data repositories, however, they must first harness and manage their data stores, and do so utilizing technologies that underscore affordability, security, and scalability.
04/15/2013 | Bull | “50% of HPC users say their largest jobs scale to 120 cores or less.” How about yours? Are your codes ready to take advantage of today’s and tomorrow’s ultra-parallel HPC systems? Download this White Paper by Analysts Intersect360 Research to see what Bull and Intel’s Center for Excellence in Parallel Programming can do for your codes.
Join HPCwire Editor Nicole Hemsoth and Dr. David Bader from Georgia Tech as they take center stage on opening night at Atlanta's first Big Data Kick Off Week, filmed in front of a live audience. Nicole and David look at the evolution of HPC, today's big data challenges, discuss real world solutions, and reveal their predictions. Exactly what does the future holds for HPC?
Join our webinar to learn how IT managers can migrate to a more resilient, flexible and scalable solution that grows with the data center. Mellanox VMS is future-proof, efficient and brings significant CAPEX and OPEX savings. The VMS is available today.