Mathematica has some predefined CUDA functions, like CUDAMap, CUDAFourier, etc. My question is how can one call functions which are part of CUDA library but not part of predefined Mathematica functions? For example, CUDA 6.0 has a large Math library; so if I am to call cuFFT from that library instead of using CUDAFourier, how would I go about doing this? Do I have to write a C code which calls the cuda library function and then call the c function from Mathematica. That seems too round about way of doing something. I could not quite find it in the documentation. Thanks very much!


1 Answer 1


CUDALink allows you to write custom kernels but unfortunately it doesn't allow to use its CUDAMemory in other functions. I found two methods to deal with it.

Documented method

You can use LibraryLink and call CUDA functions as in a regular CUDA program. For example, you can write something like fourier.cu. Then you can compile it and load "fourier_single" or "fourier" depending on your GPU architecture

CreateLibrary[{"fourier.cu"}, "libfourier", 
  "Compiler" -> NVCCCompiler, "Libraries" -> "cufft"];
AppendTo[$LibraryPath, NotebookDirectory[]];
fourier = LibraryFunctionLoad["libfourier", "fourier", {{_Complex, 1}}, {_Complex, 1}];

I don't consider this method in details because it has a big disadvantage: you have to copy your data to GPU and backward every time you call the function.

Undocumented method

Let us look through system files of Mathematica. There are brilliant examples in $InstallationDirectory/SystemFiles/Links/CUDALink/CSource/! They use undocumented MathematicaGPULibrary (WGL). This library allows to interact directly with CUDAMemory by its ID. I adapted these examples to demonstrate the usage of 1D single precision CUFFT (say, "fourier.cu" in the same directory as the notebook)

#include    <wgl.h>
#include    <wgl_cuda_runtime.h>

#define wglState                            (wglData->state)
#define wglErr                                (wglData->getError(wglData))
#define WGL_SuccessQ                        (wglErr->code == WGL_Success)
#define WGL_FailQ                            (!WGL_SuccessQ)
#define WGL_Type_RealQ(mem)                    ((mem)->type == WGL_Real_t)

#define WGL_SAFE_CALL(stmt, jmp)            stmt; if (WGL_FailQ) { goto jmp; }

WolframGPULibraryData wglData = NULL;

static int iFourier(WGL_Memory_t input, WGL_Memory_t output) {
    int err;


    if (input->type != output->type) {
        return LIBRARY_TYPE_ERROR;
    } else if (input->flattenedLength != output->flattenedLength) {

    WGL_SAFE_CALL(CUDA_Runtime_setMemoryAsInput(wglState, input, wglErr), cleanup);
    WGL_SAFE_CALL(CUDA_Runtime_setMemoryAsOutput(wglState, output, wglErr), cleanup);

    cufftHandle plan;
    cufftPlan1d(&plan, input->flattenedLength, CUFFT_C2C, 1);

        (cufftComplex *)CUDA_Runtime_getDeviceMemoryAsComplexFloat(input), 
        (cufftComplex *)CUDA_Runtime_getDeviceMemoryAsComplexFloat(output), 


    if (WGL_SuccessQ) {
        CUDA_Runtime_setMemoryAsValidOutput(wglState, output, wglErr);
    } else {
        CUDA_Runtime_setMemoryAsInvalidOutput(wglState, output, wglErr);
    CUDA_Runtime_unsetMemoryAsInput(wglState, input, wglErr);
    if (WGL_SuccessQ && err == LIBRARY_NO_ERROR) {
        return LIBRARY_NO_ERROR;
    } else {

EXTERN_C DLLEXPORT int oFourier(WolframLibraryData libData, mint Argc, MArgument * Args, MArgument Res) {
    WGL_Memory_t input, output;
    mint inputId, outputId;

    inputId                     = MArgument_getInteger(Args[0]);
    outputId                     = MArgument_getInteger(Args[1]);

    WGL_SAFE_CALL(wglData->setWolframLibraryData(wglData, libData), cleanup);

    WGL_SAFE_CALL(input = wglData->findMemory(wglData, inputId), cleanup);
    WGL_SAFE_CALL(output = wglData->findMemory(wglData, outputId), cleanup);

    err = iFourier(input, output);

    return err;

EXTERN_C DLLEXPORT int WolframGPULibrary_initialize(WolframGPULibraryData wglData0) {
    wglData = wglData0;
    return LIBRARY_NO_ERROR;

EXTERN_C DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData) {

EXTERN_C DLLEXPORT void WolframLibrary_uninitialize( ) {

Let's go back to Mathematica and compile the source and load the function

CreateLibrary[{"fourier.cu"}, "libfourier.so", 
  "TargetDirectory" -> ".", "Compiler" -> NVCCCompiler, 
  "Defines" -> "CONFIG_USE_CUDA", 
  "Libraries" -> FileNameJoin@{$CUDALinkLibraryPath, "libCUDALink_Single.so"}];

fourier = 
   "oFourier", {{_Integer, "Input"}, {_Integer, "Output"}}, 256, 
   TargetPrecision -> "Single"];


  • CONFIG_USE_CUDA flag is necessary for using CUDA in WGL. Another flag is CONFIG_USE_OPENCL.

  • Initialization of WGL fails without "Libraries" -> FileNameJoin@{$CUDALinkLibraryPath, "libCUDALink_Single.so"}. However everything will work if you evaluate any built-in CUDA* function before.

  • For double precision change Single to Double, Float to Double, C2C to Z2Z and cufftComplex to cufftDoubleComplex everywhere above and below.

  • CUFFT_INVERSE corresponds to forward Fourier with FourierParameters -> {1, 1}.

  • Input and output is _Integer because they transfer integer IDs, not the data itself.

  • Headers of WGL are located in GPUTools`Internal`$GPUToolsSystemResourcesPath

  • Built-in functions like CUDAFourier are loaded in $CUDALinkPath/CUDALink.m. It can be used as an example of loading WGL functions.

Tests (change 2^24 to an appropriate value if you don't have enough GPU memory):

x = RandomComplex[{-1 - I, 1 + I}, 2^24];
(* 268435600 *)

f1 = Fourier[x]; // AbsoluteTiming
(* {0.404416, Null} *)

gx = CUDAMemoryLoad[x, TargetPrecision -> "Single"]
gf = CUDAMemoryAllocate[Complex, Length[x], TargetPrecision -> "Single"]
(* CUDAMemory["<587248096>", "ComplexFloat"] *)
(* CUDAMemory["<1373924031>", "ComplexFloat"] *)

fourier[gx, gf]; // AbsoluteTiming
f2 = CUDAMemoryGet@gf/Sqrt@Length@x;
Norm[f2 - f1]/Norm[f1]
(* {0.010405, Null} *)
(* 2.64293*10^-7 *)

gf = CUDAFourier[gx]; // AbsoluteTiming
f3 = CUDAMemoryGet@gf;
Norm[f3 - f1]/Norm[f1]
(* {0.183712, Null} *)
(* 2.25279*10^-8 *)

CUDAMemoryUnload[gf, gx]

It works! I tested it with GTX 560 (2.0 architecture, showed above) and 9600M GT (1.1 architecture without double-precision) on Linux. Both tests show that our fourier is much faster then CPU and faster then built-in CUDAFourier. On GTX 560 CUDAFourier is relatively slow because it uses double precision and I don't know how to switch it to the single precision.

I will appreciate for any tests of WGL on other systems!

  • $\begingroup$ This looks great. CUDALink is broken (again) on my home PC but I'll test it at work on the Quadro 4000 tomorrow. $\endgroup$ Sep 21, 2014 at 10:41
  • $\begingroup$ I think this is amazing that we can use CUDA libraries directly in LibraryLink. But I don't quite understand how this work. For example, what's so special about CUDA functions that prevent them been called like ordinary functions in LibraryLink? And why do we need to interact with CUDAMemory object when using CUDA functions such as cufftExecC2C? Thanks. $\endgroup$ Sep 12, 2015 at 3:30
  • $\begingroup$ @xslittlegrass We can use LibraryLink like in my first example and call CUDA functions without CUDAMemory. However the bandwidth between CPU and GPU is very small so we have to avoid unnecessary copying to/from GPU. It is difficult to do without WGL because we cannot share GPU memory between different libraries and build-in routines (like CUDADot). $\endgroup$
    – ybeltukov
    Sep 14, 2015 at 10:35
  • $\begingroup$ @ybeltukov I guess I don't understand why the second approach avoids copying data from CPU to GPU. It looks like in the second example, gx is copied from CPU, is that right? As I understand, in both cases, we need to copy data from CPU to GPU, and after the computation on GPU finished, we need to copy data back. But why is the second approach faster? $\endgroup$ Sep 15, 2015 at 5:14
  • $\begingroup$ @xslittlegrass You are right, in this simple example the speed is the same. However, in real projects with multiple evaluation of function(s) you can permanently store a main data in the GPU memory. $\endgroup$
    – ybeltukov
    Sep 16, 2015 at 14:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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