Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have to calculate billions of dot products. My current implementation uses Compile to generate C code and it's very fast (10s of seconds). But, I was wondering if I could make it faster by using my GPU (AMD 6670). Unfortunately I don't have an NVIDIA card since it seems that CUDADot is implemented.

I'm trying to use the OpenCL DotProduct code in FileNameJoin[{$OpenCLLinkPath, "SupportFiles", "DotProduct.cl"}]

code below:

src = "/*
 * Copyright 1993-2009 NVIDIA Corporation.  All rights reserved.
 * NVIDIA Corporation and its licensors retain all intellectual property and 
 * proprietary rights in and to this software and related documentation. 
 * Any use, reproduction, disclosure, or distribution of this software 
 * and related documentation without an express license agreement from
 * NVIDIA Corporation is strictly prohibited.

 __kernel void DotProduct (__global float* a, __global float* b, __global float* c, mint iNumElements)
    // find position in global arrays
    int iGID = get_global_id(0);

    // bound check (equivalent to the limit on a 'for' loop for standard/serial C code
    if (iGID >= iNumElements)

    // process 
    int iInOffset = iGID << 2;
    c[iGID] = a[iInOffset] * b[iInOffset] 
               + a[iInOffset + 1] * b[iInOffset + 1]
               + a[iInOffset + 2] * b[iInOffset + 2]
               + a[iInOffset + 3] * b[iInOffset + 3];

I have managed to load the OpenCL function and compute a single DotProdct but this is mostely useless. I've been going over the instructions but I can't seem to figure out how to run many dot products in parallel (the whole point of using the GPU).

srcf = FileNameJoin[{$OpenCLLinkPath, "SupportFiles", "DotProduct.cl"}]

OpenCLDotProduct = 
   "DotProduct", {{"Float"}, {"Float"}, {"Float"}, _Integer}, {16, 

output = {0};
a = {1., 2., 2.};
b = {1., 2., 2.};
OpenCLDotProduct[a, b, output, 3]


output -> {{1., 2., 2.}, {1., 2., 2.}, {9.}}


This is correct, great. However, I can't figure out how to do more than a.b:

output = {0, 0};
a = {1., 2., 2.};
b = {1., 2., 2.};
c = {3., 2., 8.};
OpenCLDotProduct[{a, c}, b, output, 6]


output -> {{{1., 2., 2.}, {3., 2., 8.}}, {1., 2., 2.}, {18., 68.}}


I can't figure out what this output means or what the proper syntax is to do multiple dotproducts. I'm also not sure if this implementation is even fast. I found a blog describing another implementation but I haven't been able to compile it (http://www.openclblog.com/2012/11/opencl-and-dot-product.html).

share|improve this question
If you are using a Mac, be careful: it may be possible to run the opencl code either on the CPU or the GPU and the CPU may be the default. –  Szabolcs Mar 8 '13 at 0:39
@Szabolcs I use both but at the moment I've been working on Win7 x64. My Macbook Air is slower than my PC and it only has the integrated Intel GPU. –  s0rce Mar 8 '13 at 0:40
The problem is that you can't do this directly using the out-of-the-box solution for a dot product, since GPU-s are directly helpful for massively parallel problems, while this one directly is not. What I would do is to glue all vectors you want to multiple into two large vectors, and pass them to GPU together with a list of lengths. Then, you will have to write your own custom thread scheduler, so that effectively a thread which is finished with one vector pair is rescheduled to another one. Should be doable, but some work. Besides, you need a really powerful GPU for all this to pay off. –  Leonid Shifrin Mar 8 '13 at 0:52
As an alternative, you may look for some Map-Reduce frameworks for GPUs, such as this. I never tried those, and it is probably also quite a bit of work. But the advantage here is that problems like your can be treated by Map-Reduce more directly and universally, while the first approach I suggested is not universal. –  Leonid Shifrin Mar 8 '13 at 0:56
Also, just to stick with convention, it's best to stick with lower case starting letters for functions. I was confused for a second, thinking, "I didn't know mathematica had an OpenCLDotProduct function already defined!" –  NeuroFuzzy Mar 13 '13 at 0:08
add comment

1 Answer

up vote 2 down vote accepted

The code is expecting vectors of length 4, and an equal number of a vectors and b vectors.

For example:

n = 10^7;
a = RandomReal[{0, 10}, {n, 4}];
b = RandomReal[{0, 10}, {n, 4}];
output = ConstantArray[0., n];

 rOpenCL = Last @ OpenCLDotProduct[a, b, output, n];]

(*  {1.2324022, Null}  *)

The OpenCL is faster than Mathematica

 rMMA = MapThread[Dot, {a, b}];]

 (*  {15.2568268, Null}  *)

But compiling to C is faster

cf = Compile[{{x, _Real, 2}, {y, _Real, 2}}, MapThread[Dot, {x, y}], 
   CompilationTarget -> "C"];

AbsoluteTiming[rC = cf[a, b];]

(*  {1.0140017, Null}  *)

Also the OpenCL code is less accurate as it uses single precision:

Max[Abs[rOpenCL - rMMA]]
(*  0.0000396708  *)

Max[Abs[rC - rMMA]]
(*  0.  *)

I recompiled the OpenCL function with the actual dot product calculation completely stripped out, and found that the timing barely changed - just copying the data to and from the GPU takes almost all the time here.

share|improve this answer
Thanks very much. I guess I'll stick with C. Also if you write the compiled function to be Listable and Parallelized its a bit faster Compile[{{x, _Real, 1}, {y, _Real, 1}}, Dot[x.y], CompilationTarget -> "C", RuntimeAttributes -> {Listable}, Parallelization -> True, RuntimeOptions -> {"Speed"}] –  s0rce Mar 22 '13 at 14:51
@s0rce, you're welcome. The real benefit of having GPU code for this sort of ultra-simple procedure is when you already have data on the GPU. A while ago I wrote a CUDA kernel to do an element-by-element 2D complex array multiply. As a standalone function it was slower than just doing C = A B in Mathematica, but the overall calculation was an iterative process including 2D FFTs (where the GPU really shines). Having the multiply function in CUDA meant I could load the GPU memory once, do everything on the GPU and then copy the result back to the CPU. –  Simon Woods Mar 22 '13 at 16:08
add comment

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.