GPU Computation in Mathematica

Ulises Cervantes / Abdul Dakkak

Wolfram Research, Inc.

Built on Mathematica 8 Technologies

Made Possible by Mathematica 8 Features

CUDALink and OpenCLLink are made possible by the performance and feature improvements introduced in Mathematica 8 technologies:

LibraryLink for loading C libraries into Mathematica

SymbolicC for generating C code in Mathematica

Compiler enhancements giving speed-ups in GPU computing in Mathematica

Getting Started with GPU Programming, CUDA and OpenCL

CUDA

Programming language developed by NVIDIA in late 2007.

C/C++ like language:

__global__ void addTwo_kernel(int * arry, len) {
	int index = threadIdx.x + blockIdx * blockDim.x;
	if (index >= len) return;
    arry[index] += 2;
}

Makes the GPU a general processing unit.

Very popular — thousands of developers.

CUDA Advantages

C++ features like templates are supported.

More mature.

Offline compilation allows one to run program without performance hit.

CUDALink’s function to detect if a CUDA-capable device is available.

Needs["CUDALink`"]

CUDA’s function to detect if a CUDA-capable device is available.

CUDAQ[]

OpenCL

Programming language developed by Khronos in late 2008.

Open standard, supported by Apple, NVIDIA, and AMD/ATI.

C like language:

__kernel void addTwo_kernel(__global int * arry, len) {
	int index = get_global_id(0);
	if (index >= len) return;
    arry[index] += 2;
}

Makes the GPU a general processing unit.

New. First stable implementation in mid 2009.

OpenCL Advantages

C features only — no templates.

Open standard.

Has a specification file.

First, load the OpenCLLink application.

Needs["OpenCLLink`"]

OpenCL’s function to detect if an OpenCL-capable device is available.

OpenCLQ[]

GPU Information

This gets more information on the OpenCL devices installed.

SystemInformation[]

CUDAInformation[]

OpenCLInformation[]

Compilers

CCompilerDriver allows one to invoke C compilers from within Mathematica.

No need to maintain Make on Linux, NMake on Windows, ...

Needs["CCompilerDriver`"]

CCompilers[]//TableForm

Basic GPU Programming

CUDA

Given the following input,

listSize = 1024;
data=Range[listSize];

suppose you want to compute the following:

For[index = 1, index ≤ listSize, index++,
data[[index]]=data[[index]]+2
]

If you translate the above into C, you get the following:

void addTwo(mint * in, mint * out, mint length) {

	mint index;
	for (index = 0; index < length; index++) {
		data[index] = data[index] + 2;
	}
}

In CUDA the above C code becomes:

__global__ void addTwo(mint * in, mint * out, mint length) {

	int index = threadIdx.x + blockIdx.x*blockDim.x;
	data[index] = data[index] + 2;
}

We will use and discuss the meaning of each line next.

Understanding a Basic CUDALink Program

On the Mathematica side, all you need to let CUDALink know is the program source, the function to load, the argument types, and the block dimensions. Consider the following:

The argument types are similar in syntax to that of Compile. Input lists (or C pointers) are denoted as a list and then the type, while scalars are denoted by just the type.

So, in our example we had the following CUDA function prototype:

__global__ void addTwo(mint * in, mint length)

Therefore, our CUDAFunctionLoad command is:

CUDAFunctionLoad[src,"addTwo",{{_Integer},_Integer},256]

OpenCL

This OpenCL program adds two to each element in a dataset.

src="
__kernel void addTwo(__global mint * data, mint length) {
	int index = get_global_id(0);
	if (index < length) {
		data[index] = data[index] + 2;
	}
}";

The program is loaded using:

openclAddTwo=OpenCLFunctionLoad[src,"addTwo",{{_Integer},_Integer},16,"Defines"->{"mint"->"int"}]

This defines the input elements.

listSize = 64;
data=Range[listSize];

This runs the function on the input, returning the result.

openclAddTwo[data, listSize]

Understanding a Basic OpenCL Program

On the OpenCL source side, the function head is:

__kernel void addTwo(__global mint * data, mint length)

__kernel — tells the compiler that this functions will run on the GPU and is callable from the CPU

void — a __kernel function cannot return a value, so the return type will always be void

addTwo — is the OpenCL “kernel” name

(__global mint * data, mint length) — are the function arguments. Here mint is the Mathematica integer

Getting the index:

int index = get_global_id(0);

int index — defines the index variable as an int

get_global_id(0) — is the offset of the OpenCL thread index in the x direction

The body of the OpenCL function is

data[index] = data[index] + 2;

stores data[index] + 2 into data[index] (i.e., increments data[index] by 2)

Terminology equivalence table

We can use the following table to translate between CUDA and OpenCL.

({{"CUDA", "OpenCL"}, {"__global__…", "__kernel…"}, {"…(int * f,…)", "…(__global int * f,…)"}, {"threadIdx", "get_local_id"}, {"blockIdx", "get_group_id"}, {"blockDim", "get_local_size"}, {"threadIdx+blockIdx*blockDim", "get_global_id"}, {"__shared__", "__local"}})

Why Program Using CUDALink and OpenCLLink?

Simplified Workflow

Traditional GPU programming requires the user to manage compilation, memory, and threads.

CUDALink and OpenCLLink simplify the workflow by allowing the user to concentrate only on the algorithmic part of the program (not the setup part).

This results in a more streamlined workflow.

Platform Independence

Mathematica is platform independent and as a result CUDALink and OpenCLLink are as well. The compiler will automatically detect and use the proper compiler and flags on the system.

Using CUDALink and OpenCLLink’s abstract types, like Real_t allows the user to use the maximum precision on the system.

We also give users the ability to find out whether a kernel is launched on NVIDIA or AMD/ATI hardware. This allows for lines like the following in the kernel.

#ifdef USING_DOUBLE_PRECISIONQ
#ifdef OPENCLLINK_USING_NVIDIA
#pragma OPENCL EXTENSION cl_khr_fp64 : enable
#else /* OPENCLLINK_USING_NVIDIA */
#pragma OPENCL EXTENSION cl_amd_fp64 : enable
#endif /* OPENCLLINK_USING_NVIDIA */
#endif /* USING_DOUBLE_PRECISIONQ */

Consistent Loading Syntax

The syntax of CUDAFunctionLoad and OpenCLFunctionLoad are the same. The user only needs to change the kernel and calls from CUDA* to OpenCL*

Coupled with Mathematica‘s symbolic code generation tools, this means that users can write code that executes code that runs on either CUDA or OpenCL depending on which is available.

Access to Mathematica’s Rich Set of Programming Constructs

Functions loaded via OpenCLFunctionLoad can be used seamlessly in any Mathematica command.

Example: Using Table and GraphicsGrid — ImageAdjust

This is the source code for a basic image adjust procedure.

src="
mint xclamp(mint val, mint low, mint high) {
	return val <= low ? low : (val >= high ? high : val);
}

mint adjust(mint pixel, float lowIn, float highIn, float lowOut, float highOut, float gamma) {

	float res, val;
	val = xclamp(pixel, lowIn, highIn);
	
	res = pow((val - lowIn) / (highIn - lowIn), gamma);
	res = res * (highOut - lowOut) - lowOut;
	
	return res + 0.5f;
}

__kernel void imageAdjust(__global mint * img, mint width, mint height, mint channels, float lowIn, float highIn, float lowOut, float highOut, float gamma) {
	
	int xIndex = get_global_id(0);
	int yIndex = get_global_id(1);
	if (xIndex >= width || yIndex >= height)
		return ;
	int pos = channels * (yIndex * width + xIndex);
	for (int ii = 0; ii < channels; ii++) {
		img[pos + ii] = adjust(img[pos + ii], 255*lowIn, 255*highIn, 255*lowOut, 255*highOut, gamma);
	}
}";

This loads the function into Mathematica, it also defines OpenCLImageAdjust which accepts arguments similar to ImageAdjust.

cOpenCLImageAdjust=OpenCLFunctionLoad[src,"imageAdjust",{{_Integer},_Integer,_Integer,_Integer,"Float","Float","Float","Float","Float"},{16,16},"Defines"->{"mint"->"int"}];
OpenCLImageAdjust[img_Image,{lowIn_Real,highIn_Real},gamma_:1.0]/;Head[gamma]==Real:=OpenCLImageAdjust[img,{lowIn,highIn},{0.0,1.0},gamma]
OpenCLImageAdjust[img_Image,{lowIn_Real,highIn_Real},{lowOut_Real,highOut_Real},gamma_:1.0]:=
Module[{width,height,channels},
{height,width,channels}=Flatten[{ImageDimensions[img],ImageChannels[img]}];
cOpenCLImageAdjust[img,width,height,channels,lowIn,highIn,lowOut,highOut,gamma,{width,height}]//First
]

We can test the functionality on an image.

OpenCLImageAdjust[,{0.3,0.8}]

Or we can use Map to apply it to a set of images and values.

GraphicsGrid@Table[OpenCLImageAdjust[img,{0.3,high}],{img,{, , , }},{high,0.5,1.0,0.1}]

Example: Using Manipulate with OpenCLFunctionLoad — ColorNegate

Define the OpenCL source which color negates an input image

src="
__kernel void imageColorNegate(__global mint * img, mint width, mint height, mint channels) {
    mint ii;
	int xIndex = get_global_id(0);
	int yIndex = get_global_id(1);
	int index = channels*(xIndex + yIndex*width);
	if (xIndex < width && yIndex < height) {
		for (ii = 0; ii < channels; ii++)
			img[index+ii] = 255 - img[index+ii];
    }
}";

Load the OpenCLFunction into Mathematica.

openclColorNegate=OpenCLFunctionLoad[src,"imageColorNegate",{{_Integer},_Integer,_Integer,_Integer},{16,16},"Defines"->{"mint"->"int"}]

Run the result by calling the OpenCLFunction with appropriate inputs (note you can pass images to an OpenCLFunction).

img=;
{width,height,channels}=Flatten[{ImageDimensions[img],ImageChannels[img]}];
openclColorNegate[img,width,height,channels]

Here are a set of images.

imgs={, , , };

It does not take much to convert the above into a Manipulate (essentially we are copying and pasting code, making img a Manipulate variable)

Manipulate[
{width,height,channels}=Flatten[{ImageDimensions[img],ImageChannels[img]}];
First[openclColorNegate[img,width,height,channels]],
{img,imgs}
]

Example: Creating animations with Animate — Game of Life

Users can use dynamics to create animations of data.

src="
__kernel void gol_kernel(__global mint * prev, __global mint * nxt, mint width, mint height) {
    int xIndex = get_global_id(0);
    int yIndex = get_global_id(1);
    int index = xIndex + yIndex*width;
    int ii, jj, curr, neighbrs;
    
    if (xIndex < width && yIndex < height) {
        curr = prev[index];
        neighbrs = 0;
        for (ii = -1; ii <= 1; ii++) {
            if (xIndex + ii >= 0 && xIndex+ii < width) {
                for (jj = -1; jj <= 1; jj++) {
                    if (yIndex+jj >= 0 && yIndex+jj < height) {
                        neighbrs += prev[xIndex + ii + (yIndex+jj)*width];
                    }
                }
            }
        }
        neighbrs -= 2*curr;
        if (curr == 1)
            nxt[index] = (neighbrs == 2 || neighbrs == 3) ? 1 : 0;
        else
            nxt[index] = (neighbrs == 3) ? 1 : 0;
    }
}
";

gol=OpenCLFunctionLoad[src,"gol_kernel",{{_Integer,"Input"},{_Integer,"Output"},_Integer,_Integer},{16,16},"Defines"->{"mint"->"int"}];

initialState=RandomChoice[{0.7,0.3}->{0,1},{512,512}];
outputState=ConstantArray[0,{512,512}];

Dynamic[
Refresh[
initialState=First[gol[initialState,outputState,512,512]];ArrayPlot[initialState,ImageSize->Medium],
UpdateInterval->1/60
]
]

Scalable to Multiple GPUs or Multiple Machines

CUDALink and OpenCLLink are scalable across multiple GPUs using Mathematica’s parallel tools. They are also scalable across machines using Wolfram Lightweight Grid Mathematica.

Example: using Parallelize with CUDALink

If[$CUDADeviceCount>0,
Parallelize[

   Table[CUDAErosion[img, 5, "Device" ->1 (* $KernelID *)], {img, 

     Flatten[Permutations[{, , , }, 2]]}]]
]

CUDALink and OpenCLLink can also be used in headless setting using remote kernels.

Examples

CUDA Threads and Blocks

CUDA launches functions in a thread and block configuration. Here is a 2D configuration:

Grids, blocks, and threads are:

CUDA Indexing (1D)

Up to this point, we have briefly mentioned how CUDA does indexing.

This code will allow us to experiment with the 1D CUDA indexing:

src="
__global__ void oned(mint * res) {
   int index = threadIdx.x + blockIdx.x*blockDim.x;
   res[index] = threadIdx.x;
}
";

Loads the function:

oned=CUDAFunctionLoad[src,"oned",{{_Integer}},{16}]

This invokes the function:

res=First[oned[ConstantArray[1,64]]];

Plot the result:

ListPlot[res]

CUDA Indexing (2D)

This code will allow us to experiment with the 2D CUDA indexing:

src="
__global__ void twod(mint * res, mint width) {
   int xIndex = threadIdx.x + blockIdx.x*blockDim.x;
   int yIndex = threadIdx.y + blockIdx.y*blockDim.y;
   int index = xIndex + yIndex*width;
   res[index] = blockIdx.y;
}
";

Loads the function:

twod=CUDAFunctionLoad[src,"twod",{{_Integer},_Integer},{16,16}]

This invokes the function:

res=First[twod[ConstantArray[1,{64,64}],64]];

Visualize the result:

ArrayPlot[res]

OpenCL Indexing (2D)

src="__kernel void myKernel( __global int * global0Id, __global int * global1Id, int width, int height) {
    int xIndex = get_global_id(0);
    int yIndex = get_global_id(1);
    int index = xIndex + yIndex*width;
	if (xIndex < width && yIndex < height) {
	   global0Id[index] = get_local_id(0);
       global1Id[index] = get_local_id(1);
    }
  }";

This loads the function:

fun=OpenCLFunctionLoad[src,"myKernel",{{_Integer},{_Integer},_Integer,_Integer},{16,16}]

Define the input parameters:

width=64;
height=64;
global0Id=ConstantArray[0,{width,height}];
global1Id=ConstantArray[0,{width,height}];

Invoke the function:

res=fun[global0Id, global1Id,width,height];

Visualize the result:

ArrayPlot/@res

Removing Noise artifacts from images

Given an image with noticeable artifacts, like the one bellow

img=;

It is known that the median filter is a good solution for removing the noise.

acc=MedianFilter[img,1]

CUDALink does not come with a median filter, but one can be easily implemented. This code implements the median filter in one direction.

medianFilterCode="
__device__ unsigned char median(unsigned char x, unsigned char y, unsigned char z) {
       if (x > y) {
            if (z > y) { return x > z ? z : x; } else { return y; }
       } else {
            if (y > z) { return x > z ? x : z; } else { return y; }
       }
}
__global__ void medianFilter(unsigned char * in, unsigned char * out, int width, int height) {
   	int xIndex = threadIdx.x + blockIdx.x*blockDim.x;
   	int yIndex = threadIdx.y + blockIdx.y*blockDim.y;
   	int index = xIndex + yIndex*width;
   	if (xIndex < width && yIndex < height) {
		   unsigned char x, y, z;
		   x = xIndex == 0 ? in[index] : in[index - 1];
		   y = in[index];
		   z = xIndex == (width-1) ? in[index] : in[index + 1];
  		 out[index] = median(x,y,z);
   	}
}";

This loads the code above as a Mathematica CUDAFunction.

medianFilter=CUDAFunctionLoad[medianFilterCode,"medianFilter",{{"UnsignedByte","Input"},{"UnsignedByte","Output"},_Integer,_Integer},{16,16}]

tmp=;

First compute the median in the horizontal dimension.

tmp=First@medianFilter[tmp,tmp,Sequence@@ImageDimensions[tmp]];

Rotate the image and filter in the vertical dimension

tmp=ImageRotate[tmp,Right];

tmp=First@medianFilter[tmp,tmp,Sequence@@ImageDimensions[tmp]];

Rotate the image to retrieve the filtered image.

tmp=ImageRotate[tmp,Left]

Here are the filter results.

GraphicsRow[{img,acc,tmp}]

Task : Implement a Mean Filter

Solution

meanFilterCode="
__device__ unsigned char mean(unsigned char x, unsigned char y, unsigned char z) {
	return (x+y+z)/3.0;
}

__global__ void meanFilter(unsigned char * in, unsigned char * out, int width, int height) {
   	int xIndex = threadIdx.x + blockIdx.x*blockDim.x;
   	int yIndex = threadIdx.y + blockIdx.y*blockDim.y;
   	int index = xIndex + yIndex*width;
   	if (xIndex < width && yIndex < height) {
		   unsigned char x, y, z;
		   x = xIndex == 0 ? in[index] : in[index - 1];
		   y = in[index];
		   z = xIndex == (width-1) ? in[index] : in[index + 1];
  		 out[index] = mean(x,y,z);
   	}
}";

meanFilter=CUDAFunctionLoad[meanFilterCode,"meanFilter",{{"UnsignedByte","Input"},{"UnsignedByte","Output"},_Integer,_Integer},{16,16}]

tmp=;

First compute the mean in the horizontal dimension.

tmp=First@meanFilter[tmp,tmp,Sequence@@ImageDimensions[tmp]];

Rotate the image and filter in the vertical dimension

tmp=ImageRotate[tmp,Right];

tmp=First@meanFilter[tmp,tmp,Sequence@@ImageDimensions[tmp]];

Rotate the image to retrieve the filtered image.

tmp=ImageRotate[tmp,Left]

Task : Implement a Median and Mean Filter in OpenCL

Finite Difference

A common operation in scientific computation is to fine the derivative of a function.

One can discritize space and approximate the derivative.

Recall from Calculus the definition of a derivative of a function f:

f'(x)= (lim_(h->0))(((f (x-h)-f (x)))/h)

Mathematica has a Derivative function, but let’s define our own using the above definition by first discretizing the function and calculating the differences. We choose the function Tan[x]

f[x_] := Tan[x]

We set h to a small value, 0.005 is small enough:

h=0.005

We discretize our function f

discretizedf=Map[f,Range[-2π,2π,h]];

We then compute our difference

diff=MapThread[(#2-#1)/h&,{discretizedf[[;;-2]],discretizedf[[2;;]]}];

To plot, we compute

fplot=Table[{-2π+ii*h,discretizedf[[ii]]},{ii,1,Length[discretizedf]}];
diffplot=Table[{-2π+ii*h,diff[[ii]]},{ii,1,Length[diff]}];

Plotting gives us

ListPlot[{fplot,diffplot}]

Being very similarly to Mathematica’s results with convergence if h is sufficiently small.

Plot[{f[x],f'[x]},{x,-2π,2π},PlotRangeClipping->False,ClippingStyle->None]

Let’s do this computation using CUDA.

Finite Difference Using Backward Difference

Backward difference:

f'(x)≈f(x)-f(x-h)

finiteDifferenceSrc="
__global__ void finiteDifference(float * diff, float h, int listSize) {
	int index = threadIdx.x + blockIdx.x * blockDim.x;
    float f_n = tanf(((float) index) / h);
    float f_n1 = tanf((index - 1.0) / h);
	if( index < listSize) {
		diff[index] = (f_n - f_n1) / h;
	}
}";

Load the function:

finiteDifference=CUDAFunctionLoad[finiteDifferenceSrc,"finiteDifference",{{"Float"}, "Float",_Integer},16]

Run the function:

res=First[finiteDifference[ConstantArray[1.0,1024],100,1024]];

Plot the results:

ListPlot[res]

Passing values to Finite Difference

Code is same as before, with the exception that we do not compute the function on the GPU:

finiteDifferenceSrc="
__global__ void finiteDifference(float * diffOut, float * diffIn, float h, int listSize) {
	int index = threadIdx.x + blockIdx.x * blockDim.x;
    float f_n = diffIn[index];
    float f_n1 = diffIn[index - 1];
	if( index < listSize) {
		diffOut[index] = (f_n - f_n1)/h;
	}
}";

Load the function:

finiteDifference=CUDAFunctionLoad[finiteDifferenceSrc,"finiteDifference",{{"Float"},{"Float","Input"},"Float",_Integer},16]

Run the function:

res=First[finiteDifference[ConstantArray[1.0,1024],N[Tan[Range[1024]/100]],1/100.,1024]];

ListPlot[res]

Task : Implement a Second Order Finite Difference

f''(x)≈((f(x + h) - 2 f(x) +f(x-h))/(h^2))

Solution :

finiteDifferenceSrc="
__global__ void finiteDifference(float * diffOut, float * diffIn, float h, int listSize) {
	int index = threadIdx.x + blockIdx.x * blockDim.x;
    float f_n0 = index == 0 ? diffIn[0] : diffIn[index-1];
	float f_n = index < listSize ? diffIn[index] : diffIn[listSize-1];
    float f_n1 = index+1 < listSize ? diffIn[index+1] : diffIn[listSize-1];
	if( index < listSize) {
		diffOut[index] = (f_n1 - 2 * f_n + f_n0)/h;
	}
}";

Load the function:

finiteDifference=CUDAFunctionLoad[finiteDifferenceSrc,"finiteDifference",{{"Float"},{"Float","Input"},"Float",_Integer},16]

Get discritize function:

vals=N[Tan[Range[1024]/100]];

Run the function:

res=First[finiteDifference[ConstantArray[1.0,1024],vals,1/100.^2,1024]];

Visualize the results:

ListPlot[res]

Complex Dynamical System

Julia Set of the following rational function

(Subscript[f, c])(z)=c+(z^2)

ComplexExpand[(x+I y)^2]

code="
__global__ void julia_kernel(Real_t * set, mint width, mint height, Real_t cx, Real_t cy) {
   int xIndex = threadIdx.x + blockIdx.x*blockDim.x;
   int yIndex = threadIdx.y + blockIdx.y*blockDim.y;
   int ii;

   Real_t x = ZOOM_LEVEL*(width/2 - xIndex);
   Real_t y = ZOOM_LEVEL*(height/2 - yIndex);
   Real_t tmp;
   Real_t c;

   if (xIndex < width && yIndex < height) {
       for (ii = 0; ii < MAX_ITERATIONS && x*x + y*y < BAILOUT; ii++) {
			tmp = x*x - y*y + cx;
			y = 2*x*y + cy;
			x = tmp;
		}
		c = logf(static_cast<Real_t>(0.1) + sqrtf(x*x + y*y));
		set[xIndex + yIndex*width] = c;
    }
}
";

{width,height}={512,512};
jset=CUDAMemoryAllocate[Real,{height,width}];

JuliaCalculate=CUDAFunctionLoad[code,"julia_kernel",{{_Real,_,"Output"},_Integer,_Integer, _Real, _Real},{16,16},"Defines"->{"MAX_ITERATIONS"->10,"ZOOM_LEVEL"->"0.0050","BAILOUT"->"4.0"}]

Manipulate[
JuliaCalculate[jset, width, height, c[[1]], c[[2]], {width, height}];
ReliefPlot[Reverse@CUDAMemoryGet[jset], ColorFunction->"Rainbow",DataRange->{{-2.0, 2.0}, {-2.0, 2.0}},ImageSize->256],
{{c, {0, 1}}, {-2, -2}, {2, 2}, Locator}]

Task : Compute the Julia set for the following iterative function

(Subscript[f, c])(z)=(z^3)+c

Hint

Note  the following

ComplexExpand[(x+ I y)^3]

Solution

code="
__global__ void julia_kernel(Real_t * set, mint width, mint height, Real_t cx, Real_t cy) {
   int xIndex = threadIdx.x + blockIdx.x*blockDim.x;
   int yIndex = threadIdx.y + blockIdx.y*blockDim.y;
   int ii;

   Real_t x = ZOOM_LEVEL*(width/2 - xIndex);
   Real_t y = ZOOM_LEVEL*(height/2 - yIndex);
   Real_t tmp;
   Real_t c;

   if (xIndex < width && yIndex < height) {
       for (ii = 0; ii < MAX_ITERATIONS && x*x + y*y < BAILOUT; ii++) {
			tmp = x*x*x - 3*x*y*y + cx;
			y = 3*x*x*y - y*y*y + cy;
			x = tmp;
		}
		c = logf(static_cast<Real_t>(0.1) + sqrtf(x*x + y*y));
		set[xIndex + yIndex*width] = c;
    }
}
";

{width,height}={512,512};
jset=CUDAMemoryAllocate[Real,{height,width}];

JuliaCalculate=CUDAFunctionLoad[code,"julia_kernel",{{_Real,_,"Output"},_Integer,_Integer, _Real, _Real},{16,16},"Defines"->{"MAX_ITERATIONS"->10,"ZOOM_LEVEL"->"0.0050","BAILOUT"->"4.0"}]

Manipulate[
JuliaCalculate[jset, width, height, c[[1]], c[[2]], {width, height}];
ReliefPlot[Reverse@CUDAMemoryGet[jset], ColorFunction->"Rainbow",DataRange->{{-2.0, 2.0}, {-2.0, 2.0}},ImageSize->256],
{{c, {0, 1}}, {-2, -2}, {2, 2}, Locator}]

Ray Tracing

Define Kernel

src="
#ifdef USING_DOUBLE_PRECISIONQ
#ifndef OPENCLLINK_USING_ATI
#pragma OPENCL EXTENSION cl_khr_fp64 : enable
#else /* OPENCLLINK_USING_ATI */
#pragma OPENCL EXTENSION cl_amd_fp64 : enable
#endif /* OPENCLLINK_USING_ATI */
#endif /* USING_DOUBLE_PRECISIONQ */

bool intersectsQ(__global Real_t * center, Real_t rad, Real_t * ray) {
	Real_t dx = ray[0] - center[0];
	Real_t dy = ray[1] - center[1];
	Real_t v = dx*dx + dy*dy;
	if (v < rad*rad) {
		return true;
	} else {
		return false;
	}
}
Real_t hitTest(__global Real_t * center, Real_t rad, Real_t * ray) {
	Real_t dx = ray[0] - center[0];
	Real_t dy = ray[1] - center[1];
	Real_t v = dx*dx + dy*dy;
	if (v < rad*rad) {
		return sqrt(rad*rad - v);
	} else {
		return 0;
	}
}
__kernel void rayTrace(__global Real_t * img, mint numSpheres, __global Real_t * centers, __global Real_t * radii, __global Real_t * colors, mint width, mint height, mint channels) {
	int xIndex = get_global_id(0);
	int yIndex = get_global_id(1);
	int index = xIndex + yIndex*width;
	if (xIndex < width && yIndex < height) {
		Real_t ray[2];
		int offset = index * channels;
		Real_t max = - INFINITY;
		Real_t t, k, r = 0, g = 0, b = 0;
		ray[0] = xIndex - width/2;
		ray[1] = yIndex - height/2;
		for (int ii = 0; ii < numSpheres; ii++) {
			if (intersectsQ(&centers[ii * 3], radii[ii], ray)) {
				t = hitTest(&centers[ii * 3], radii[ii], ray);
				k = t + centers[ii * 3 + 2];
				t /= radii[ii];
				if (k > max) {
					r = colors[ii * 3    ] * t;
					g = colors[ii * 3 + 1] * t;
					b = colors[ii * 3 + 2] * t;
					max = k;
				}
			}
		}
		img[offset] = r;
		img[offset+1] = g;
		img[offset+2] = b;
	}
}
";

Load Kernel

raytrace=OpenCLFunctionLoad[src,"rayTrace",{{_Real,"Output"},_Integer,{_Real,"Input"},{_Real,"Input"},{_Real,"Input"},_Integer,_Integer,_Integer},{16,16},"Defines"->{"mint"->"int"}];

Define Input Parameters

width=height=1024;
channels=3;
img=ConstantArray[0.0,{height,width,channels}];
numSpheres=1000;
spheres=Table[Sphere[{RandomReal[{-1000,1000}],RandomReal[{-2000,2000}],RandomReal[{-1000,33}]},RandomReal[{20,100}]],{numSpheres}];
centers=spheres[[All,1]];
radii=spheres[[All,2]];
colors=RandomReal[1,{numSpheres,3}];

Compute and Visualize

res=First@raytrace[img,numSpheres,centers,radii,colors,width,height,channels];

Image[res]

Make it Dynamic

Manipulate[
res=First[raytrace[img,numSpheres,(#+{x,y,0})&/@centers,radii,colors,width,height,channels]];
Image[res,ImageSize->256]
,
{{x,0},-1000,1000},
{{y,0},-1000,1000}
]

The following implements triplex Mandelbrot (or “Mandelbulb”) — an analog to the Mandelbrot set in 3D.

srcf=FileNameJoin[{$CUDALinkPath,"SupportFiles","mandelbulb.cu"}];

FilePrint[srcf]

height=480;
width=640;
iconfig={width,height,1,0,1,6};
config={0.001,0.0,0.0,0.0,8.0,15.0,10.0,5.0};
camera={{2.0,2.0,2.0},{0.0,0.0,0.0}};
AppendTo[camera,Normalize[camera[[2]]-camera[[1]]]];
AppendTo[camera,0.75*Normalize[Cross[camera[[3]],{0.0,1.0,0.0}]]];
AppendTo[camera,0.75*Normalize[Cross[camera[[4]],camera[[3]]]]];
config=Join[{config,Flatten[camera]}];
pixelsMem=CUDAMemoryAllocate["Float",{height,width,3}];

mandelbulb=CUDAFunctionLoad[{srcf},"MandelbulbGPU",{{"Float",_,"Output"},{"Float",_,"Input"},{_Integer,_,"Input"},_Integer,"Float","Float"},{16},"UnmangleCode"->False];
mandelbulb[pixelsMem,Flatten[config],iconfig,0,0.0,0.0,{width*height*3}];
pixels=CUDAMemoryGet[pixelsMem];
Image[pixels]

Recommended Next Steps

What is new in Mathematica 8

CUDALink white paper

Check CUDALink and OpenCLLink documentation pages

CUDALink functions, programming, and applications guides

Extensive real-world examples in all reference pages

Recommended Next Steps

Try Mathematica for free:
http://www.wolfram.com/mathematica/trial

Schedule a technical demo:
www.wolfram.com/mathematica-demo

Learn more about CUDA and OpenCL programming in Mathematica:
U.S. and Canada:						Europe:
1-800-WOLFRAM (965-3726)					+44-(0)1993-883400
info@wolfram.com						info@wolfram.co.uk

Outside U.S. and Canada (except Europe and Asia):		Asia:
+1-217-398-0700						+81-(0)3-3518-2880
info@wolfram.com						info@wolfram.co.jp