Well, here is a recent story with one of our customers that accomplished 98X performance speed up with Jacket in 16 days.  Of the 16 days, 15 days were spent sending emails back and forth about the problem and less than a day was spent getting the customer code in Jacket and running some performance tests!  Who would have imagined GPU programming with performance in 1 day.  Happy Reading.

<strong>Day 1:</strong>

Customer uses inverse radon (or iradon in MATLAB terms) extensively for their back projection algorithms.  They would like to know when the iradon function will be available/supported in Jacket.

AccelerEyes product management informs the customer that the inverse radon algorithm used in MATLAB is based on the filtered back project algorithm.  The key elements of that algorithm are interpolation, fft and frequency domain filtering which are all supported in Jacket and should provide decent speed-up for large enough data sizes.

<strong>Day 4:</strong>

The customer provides email response to product management indicating a couple of things:
<ol>
	<li>They would be interested in working with AccelerEyes and Jacket on an initial implementation of inverse radon given the information provided in the Day 1 response.</li>
	<li>If the initial implementation demonstrated a substantial savings in execution time, the customer could allocate the appropriate internal resource to leverage Jacket and GPUs for their needs.</li>
	<li>They had some (very limited) experience with programming GPUs directly (we assume CUDA), and had determined that the approach would require far too much development effort from a relatively small staff of software developers.</li>
	<li>The customer is ready to work on a prototype/pilot to get some results.</li>
</ol>
<strong>Day 17:</strong>

<span style="text-decoration: underline;">6:49am: </span> After 13 days, the customers sends a follow up email to product management and corporate email box asking whether or not AccelerEyes was interested in working with the customer on the inverse radon prototype/pilot.
<blockquote>“turns out the spam-filters are not always your best friend as the customers email on Day 4 got caught in the corporate spam filter and the product manager never got the message” – morale of the story is follow up email with a phone call in time sensitive situations!</blockquote>
<span style="text-decoration: underline;">11:31am:</span> After profuse apologies by product management regarding the delay, an exchange of information about the customer code begins.  A review of the use case for the function is undergone (use case = the MATLAB code developed by the customer for the application).   The MATLAB code for the iradon function itself , in addition to the MATLAB built-in functions are also reviewed with the customer to understand all aspects of the application/algorithm.

<span style="text-decoration: underline;">2:31pm: </span> It becomes clear what needs to be done with Jacket and how parts of the customer’s MATLAB code and other m-files need to be prepared for Jacket.
<blockquote><strong>For all Jacket users:</strong> many of you are interested in various MATLAB functions to be implemented in Jacket to ease the programming effort to leverage GPUs, increase performance on more applications and in general achieve more results.  Although the outline below of what needed to be done to get iradon working for this customer was not a design goal and is NOT for all Jacket users, this approach can be used by many Jacket customers to achieve results today without waiting for the specific function or functions to show up in Jacket documentation and fully supported.

The approach:

Step 1: Find the function in the MATLAB source tree that is bottlenecking your performance – in this case “iradon.m”

Step 2: Save the function or .m file as a GPU function like foo.m, in this case the function was saved as “giradon.m”

Step 3: In your foo.m function (this case giradon.m), pass the input variables from your application into the new function you are creating (example - function [img,H] = iradon(varargin) becomes function [img,H] = giradon(varargin))

Step 4: Add a “g” pre-fix to zeros when allocating memory which will allocate memory on the GPU for execution – if applicable

Step 5: Assign variables to be GPU variables p=gsingle(p); that will be used during execution of the most computational intensive calls within the m-file.  It is only necessary to make GPU variables of the data that goes to the GPU and where the Jacket function is supported – fft for example.  In this particalur case the following code was added to the iradon.m file to make it Jacket compatible for ‘linear’ case.
<p style="padding-left: 30px;"><span style="color: #008000;">%% Change 1: Add below line 158</span></p>
<p style="padding-left: 30px;"><span style="color: #008000;">%% Converting to gsingle values for GPU usage :</span></p>
<p style="padding-left: 30px;">costheta=gsingle(costheta);</p>
<p style="padding-left: 30px;">sintheta=gsingle(sintheta);</p>
<p style="padding-left: 30px;">p=gsingle(p);</p>
<p style="padding-left: 30px;">x=gsingle(x);</p>
<p style="padding-left: 30px;">y=gsingle(y);</p>
<p style="padding-left: 30px;">ctrIdx=gsingle(ctrIdx);</p>
Step 6: Vectorize code where possible. Especially where for-loops are involved.  In the iradon case the customer vectorized the for-loop in the function filterProjections() as shown below:
<p style="padding-left: 30px;"><span style="color: #008000;">%% Change 2: Lines 204-206</span></p>
<p style="padding-left: 30px;"><span style="color: #008000;">%% Following 'for' loop unrolled using vectorizations</span></p>
<p style="padding-left: 30px;"><span style="color: #008000;">% for i = 1:size(p,2)</span></p>
<p style="padding-left: 30px;"><span style="color: #008000;">%    p(:,i) = p(:,i).*H; % frequency domain filtering</span></p>
<p style="padding-left: 30px;"><span style="color: #008000;">% end</span></p>
<p style="padding-left: 30px;"><span style="color: #008000;">%% After Vectorizing the above for-loop :</span></p>
<p style="padding-left: 60px;">p(:,1:size(p,2))=p(:,1:size(p,2)).*repmat(H,1,size(p,2));</p>
<p style="padding-left: 30px;"><span style="color: #008000;">%%</span></p>
Step 6:  Open your MATLAB code that calls the bottleneck function.  You have two options here, you can program your MATLAB code to always use GPUs by changing the call to iradon.com (in this example) to foo.m or giradon.m or you can set a flag that will ask the user if they have GPUs or not.
<p style="padding-left: 30px;">If no, they follow the normal execution of the code using the standard m-files.</p>
<p style="padding-left: 30px;">If yes, you will need to follow an execution path that has the iradon.m file (from this example) changed to giradon.m or foo.m.</p>
Once you have modified your MATLAB application, save and run the code with the GPU option on.</blockquote>
<span style="text-decoration: underline;">3:54pm:</span> Customer needs to test and benchmark changes to code and chooses to use the following data sizes for initial testing while a determination is made of the most typical parameters used in reconstructions :
<p style="padding-left: 30px;">P = phantom(128);</p>
<p style="padding-left: 30px;">R = radon(P,0:179);</p>
<p style="padding-left: 30px;">I1 = iradon(R,0:179);</p>
<p style="padding-left: 30px;">I2 = iradon(R,0:179,'linear','none');</p>
<strong>Day 18:</strong>

<span style="text-decoration: underline;">10:17am:</span> Customer determines that the most typical parameters would be to increase the phantom size to <strong>1024 X 1024 ( P = phantom(1024) )</strong> from P = phantom(128).  It turned out the number of angles used is 180 so the initial use case of 0:179 was in line with the requirement.

<span style="text-decoration: underline;">5:55pm:</span> The following performance results are achieved with Jacket on GPUs using two different hardware configurations and various dimensions of data.

System 1:

CPU: Intel Xeon Dual Core 2.8 GHz, 1GB RAM, (2 Cores)

GPU: GeForce GTX 295, 1212 MHz, 895 MB VRAM, Compute 1.3

OS: Window XP (32-bit)

MATLAB Version: R2009b

Jacket Version: 1.2.2
<table border="0" cellspacing="0" cellpadding="0" width="606">
<tbody>
<tr>
<td width="96" valign="top">Phantom Dim</td>
<td width="108" valign="top">Projections Dim</td>
<td width="108" valign="top">'iradon' runtime (sec)</td>
<td width="114" valign="top">'giradon' runtime

(sec)</td>
<td width="72" valign="top">Speed GPU vs CPU</td>
<td width="108" valign="top">Error Norm (inf)</td>
</tr>
<tr>
<td width="96" valign="top">128x128</td>
<td width="108" valign="top">185x180</td>
<td width="108" valign="top">0.64</td>
<td width="114" valign="top">0.98</td>
<td width="72" valign="top">0.65</td>
<td width="108" valign="top">1.77E-006</td>
</tr>
<tr>
<td width="96" valign="top">256x256</td>
<td width="108" valign="top">367x180</td>
<td width="108" valign="top">2.45</td>
<td width="114" valign="top">0.91</td>
<td width="72" valign="top">2.69</td>
<td width="108" valign="top">3.35E-006</td>
</tr>
<tr>
<td width="96" valign="top">512x512</td>
<td width="108" valign="top">729x180</td>
<td width="108" valign="top">10.39</td>
<td width="114" valign="top">0.97</td>
<td width="72" valign="top">10.71</td>
<td width="108" valign="top">6.82E-006</td>
</tr>
<tr>
<td width="96" valign="top">1024x1024</td>
<td width="108" valign="top">1453x180</td>
<td width="108" valign="top">41.46</td>
<td width="114" valign="top">1.15</td>
<td width="72" valign="top">36.05</td>
<td width="108" valign="top">1.84E-005</td>
</tr>
<tr>
<td width="96" valign="top">2048x2048</td>
<td width="108" valign="top">2901x180</td>
<td width="108" valign="top">167.07</td>
<td width="114" valign="top">1.69</td>
<td width="72" valign="top">98.86</td>
<td width="108" valign="top">3.12E-005</td>
</tr>
<tr>
<td width="96" valign="top"></td>
<td width="108" valign="top"></td>
<td width="108" valign="top"></td>
<td width="114" valign="top"></td>
<td width="72" valign="top"></td>
<td width="108" valign="top"></td>
</tr>
</tbody>
</table>
System 2:

CPU: Intel(R) Dual Quad Core Xeon(R) CPUs W5580 @ 3.20GHz, 32GB RAM, (8 cores)

GPU: Tesla C1060, 1265 MHz, 4095 MB VRAM, Compute 1.3

OS: Linux Fedora 10 (32-bit)

MATLAB Version: R2009b

Jacket Version: 1.2.2
<table border="0" cellspacing="0" cellpadding="0" width="606">
<tbody>
<tr>
<td width="96" valign="top">Phantom Dim</td>
<td width="108" valign="top">Projections   Dim</td>
<td width="108" valign="top">'iradon'   runtime (sec)</td>
<td width="114" valign="top">'giradon' runtime

(sec)</td>
<td width="72" valign="top">Speed GPU vs   CPU</td>
<td width="108" valign="top">Error Norm   (inf)</td>
</tr>
<tr>
<td width="96" valign="top">128x128</td>
<td width="108" valign="top">185x180</td>
<td width="108" valign="top">0.1</td>
<td width="114" valign="top">0.4</td>
<td width="72" valign="top">0.25</td>
<td width="108" valign="top">1.77E-006</td>
</tr>
<tr>
<td width="96" valign="top">256x256</td>
<td width="108" valign="top">367x180</td>
<td width="108" valign="top">0.38</td>
<td width="114" valign="top">0.41</td>
<td width="72" valign="top">0.93</td>
<td width="108" valign="top">3.35E-006</td>
</tr>
<tr>
<td width="96" valign="top">512x512</td>
<td width="108" valign="top">729x180</td>
<td width="108" valign="top">1.82</td>
<td width="114" valign="top">0.43</td>
<td width="72" valign="top">4.23</td>
<td width="108" valign="top">6.82E-006</td>
</tr>
<tr>
<td width="96" valign="top">1024x1024</td>
<td width="108" valign="top">1453x180</td>
<td width="108" valign="top">8.35</td>
<td width="114" valign="top">0.51</td>
<td width="72" valign="top">16.37</td>
<td width="108" valign="top">1.84E-005</td>
</tr>
<tr>
<td width="96" valign="top">2048x2048</td>
<td width="108" valign="top">2901x180</td>
<td width="108" valign="top">38.14</td>
<td width="114" valign="top">1.93</td>
<td width="72" valign="top">19.76</td>
<td width="108" valign="top">3.12E-005</td>
</tr>
</tbody>
</table>
The “iradon” runtimes were run on the CPU where the “giradon” runtimes were leveraging the GPU.  Of course the title of this blog post was a little misleading as it is not true that 98X speed up was the case in every data size and every configuration as performance does vary in all applications based on code, data size, hardware and other factors.  In this case, an everyday laptop with a GeForce GPU can provide almost 100X performance improvement versus a garden variety Intel Core Duo 2.8 GHz processor for the right size data problem.  But even when you run on state of the art hardware, both CPU and GPU, and the data size is appropriate for the GPU architectures, it is not unusual to get 10 to 20X performance improvement.

There is actually a bigger story here than the performance improvements and speed up of this code!  It has to do with productivity.  This customer had already determined that they did not have the resources, knowledge and expertise to program CUDA to leverage the GPU resources.  Jacket provides a platform for engineers, scientists and analysts, that are familiar with MATLAB, to quickly (in this case, <strong>a couple of hours</strong>) make current applications and algorithms work on GPUs to achieve substantial performance gains, even when the latest and greatest hardware is not immediately at ones disposal.  Shortening the time to solution is one of the major value points for Jacket and should be considered even if an option to program in CUDA is available.

PS – one more point to the story, email is great but today’s spam filters can really disrupt communication if alternative methods are not used in conjunction with our new way of life.