Showing posts with label IDL. Show all posts
Showing posts with label IDL. Show all posts

Tuesday, May 7, 2013

Custom kernels in GPULib 1.6

One of the new features of GPULib 1.6 is the ability to load and execute custom CUDA kernels at runtime. This means you can run CUDA code does not know about GPULib (or even IDL). Let's run through what is required to do this for a simple example.

Our example will simply add two n element arrays and store the result in a third array. The kernel for this will be written in a file `vectorAdd_kernel.cu` and should be straight-forward for anyone writing CUDA:

  extern "C" __global__ void VecAdd_kernel(const float* A,
                                           const float* B,
                                           float* C,
                                           int N)
  {
      int i = blockDim.x * blockIdx.x + threadIdx.x;
      if (i < N)
          C[i] = A[i] + B[i];
  }
This code is then compiled to PTX code using the --ptx flag of nvcc:
  $ nvcc --ptx vectorAdd_kernel.cu
This produces the file `vectorAdd_Kernel.ptx` which can be loaded and executed with a couple of new GPULib routines. First, load the source code of the kernel:
  ptx_source = gpu_read_ptxfile('vectorAdd_kernel.ptx')
This simply loads the source code from the PTX file into an IDL string. The source code could contain multiple kernel definitions; the module represents the entirety of the code, while the function represents one particular kernel. We call GPULib routines to load the module and function:
  module = gpuLoadModule(ptx_source, error=err)
  kernel = gpuLoadFunction(module, 'VecAdd_kernel', error=err)
Next, we create two arrays to add and a third array to contain the result:
  n = 20L
  dx = gpuFindgen(n)
  dy = gpuFindgen(n)
  dz = gpuFltarr(n)
Lastly, we are ready to call our kernel. The grid and block structure is set through the GRID and BLOCKS keywords, while otherwise we just pass the kernel and its arguments to GPUEXECUTEFUNCTION:
  nThreadsPerBlock = 256L
  nBlocksPerGrid = (n + nThreadsPerBlock - 1L) / nThreadsPerBlock

  gpuExecuteFunction, kernel, $
                      dx->_getHandle(), $
                      dy->_getHandle(), $
                      dz->_getHandle(), $
                      n, $
                      GRID=[nBlocksPerGrid], $
                      BLOCKS=[nThreadsPerBlock], $
                      ERROR=err
There is a SHARED_MEMORY keyword to pass the amount of shared memory that your kernel requires, but that is not needed for our simple example.

Friday, May 3, 2013

Linear algebra in GPULib 1.6

GPULib 1.6 adds 100+ linear algebra routines available through MAGMA, a GPU accelerated LAPACK implementation.

The IDL bindings for the low-level MAGMA routines are automatically generated wrappers for the standard LAPACK interface, but GPULib also offers a higher-level routine mimicking the LA_ routines present in the IDL library. GPUINVERT is comparable to LA_INVERT, inverting a non-singular matrix in a single call:

  dinverse = gpuInvert(da, lhs=dinverse)
While these routines document their technique for performing their operations, they do not require much user knowledge of the underlying techniques used. More high-level routines similar to the LA_ are planned for future GPULib versions.

Calling the low-level routines directly in MAGMA requires more knowledge of the mathematical algorithms and a bit more cumbersome notation. For example, the following example shows how to call the MAGMA version of SGELS:

  status = gpusgels((byte('N'))[0], $
                    m, n, nrhs, $
                    da->_getHandle(), lda, $
                    db->_getHandle(), ldb, $
                    work, lwork, $
                    info)
SGELS solves overdetermined or underdetermined real linear systems using a QR or LQ factorization of a matrix of full rank. This low-level operation has options for performing multiple solves in a call and operating on the transpose.

The calling syntax for the MAGMA routines is in the API documentation for GPULib as well as in lib/magma_routines.h in the GPULib distribution.

Wednesday, May 1, 2013

GPULib 1.6 released

Now available for download from the Tech-X website, GPULib 1.6 offers some exciting new features and a lot more stability. Here’s the brief rundown of what's new:

  • All platforms, Windows, Linux, and OS X, are now distributed as binaries. No building from source required!
  • Added wrappers for MAGMA, a GPU accelerated LAPACK library. This makes over 100 low-level linear algebra routines available. In addition to these low-level routines, we have one high-level routine, GPUINVERT, that is similar to the IDL routine LA_INVERT. We intend to write analogs to the rest of the LA_ IDL routines as well. This feature requires the Basic membership plan.
  • GPULib can now load and execute custom CUDA kernels without having to link to GPULib; you just compile your kernel to a .ptx file. GPULib provides routines to load and execute that kernel at runtime. This feature requires the Basic membership plan.
  • Support for CUDA 5.0.
  • Added support for up to 8-dimensional arrays. This is basic support for higher-dimensional arrays in GPULib, not every routine understands them yet.
  • Added optimized scalar/array operations. No more tricks just to add a scalar to an array!
  • Fixed bugs in GPUHISTOGRAM and GPUHIST_2D.
  • Fixed bugs in GPUATAN2 for complex and double complex.
  • Fixed bugs in GPUREAL.

I have planned a few more posts with details and examples about some of the big features above, like the MAGMA wrappers and the custom kernels.

A lot of work was done on infrastructure to make releasing an easier process, hopefully resulting in more frequent updates in the future. We have plans for some very exciting features in the coming year!

Monday, August 27, 2012

Experiments with OpenCL

I've been looking at what it would take to convert GPULib from CUDA to OpenCL. In the effort, I have created a simple DLM that provides some basic arithmetic operations and some more advanced features. For example, we can create variables on the GPU in standard ways:

IDL> dx = cl_findgen(10)
IDL> dresult = cl_fltarr(10)
We can also transfer data from the CPU to the GPU:
IDL> hy = 2.0 * findgen(10)
IDL> dy = cl_putvar(hy)
Performing a simple operation is easy (though I will eventually change the form of this call to return the value of the result):
IDL> status = cl_add(dx, dy, dresult)
Retrieve the result:
IDL> result = cl_getvar(dresult)
I have implemented memory operations for allocating, transferring data, and freeing memory. I also have routines for obtaining metadata on variables, as well as the available OpenCL devices and platforms.

OpenCL kernels are always constructed from strings, so it is easy to make custom kernels on-the-fly to evaluate user entered expressions. In the next example we will evaluate the expression:

dz = 2. * dx + sin(dy)
First, allocate and initialize the inputs/outputs:
dx = cl_findgen(10)
dy = cl_putvar(2. * findgen(n))
dz = cl_fltarr(10, /nozero)
Next, create the kernel using a string representing the expression, a string array of the variable names in the expression, and an array of type codes of the variables in the expression:
kernel = cl_compile('z[i] = 2. * x[i] + sin(y[i])', $
                    ['x', 'y', 'z'], $
                    lonarr(3) + 4L, $
                    /simple)
The SIMPLE keyword indicates that we are just filling in an expression in a simple element-wise kernel. To run the kernel on our data, call CL_EXECUTE with the kernel identifier and a structure containing the values for the variables in the expression:
status = cl_execute(kernel, { x: dx, y: dy, z: dz })
Finally, retrieve the results as before:
z = cl_getvar(dz)
Without the SIMPLE keyword set, the user can specify the entire full code for the kernel.

There are some great advantages to using OpenCL:

  1. OpenCL is an open standard.
  2. OpenCL is not limited to NVIDIA GPUs (or even GPUs at all, OpenCL can run on various platforms including CPUs).
  3. It is easy to create on-the-fly kernels.
  4. OpenCL uses the same concepts as CUDA, so many kernels can be translated fairly easily.
Unfortunately, there is one big disadvantage right now:
  1. OpenCL libraries are not as far along as CUDA libraries, both those provided by NVIDIA and third-parties.
This would mean that anything beyond basic arithmetic in the current GPULib would be eliminated right now. I don't think it will be long before OpenCL catches up in this regard, but this is a serious issue currently.

Thursday, January 22, 2009

32- vs. 64-bit IDL on Mac OS X

IDL 7.0.4 was recently released, providing full 64-bit support for Mac OS X, and some users who have upgraded have asked how this impacts GPULib, as using the library when running IDL in 64-bit mode doesn't work.

Currently, in order to build 64-bit code using GCC on OS X, one needs to set an additional flag, depending on the architecture targeted. We don't do this at the moment for GPULib, so it defaults to building the more compatible 32-bit libraries. Full details at

http://developer.apple.com/documentation/Darwin/Conceptual/64bitPorting/building/chapter_5_section_2.html


Apple goes so far as to say, "You should transition your application to a 64-bit executable format only when the 64-bit environment offers a compelling advantage for your specific application." We haven't switched to using 64-bit libraries with GPULib on Mac OS X, since thus far there hasn't been a need. If you do need the extra address space using GPULib, please let us know! Otherwise, you can add the "-32" flag to IDL 7.0.4 and run IDL in 32-bit mode, which will be compatible with the current release of GPULib.

Tuesday, January 13, 2009

Performing a shift in IDL with GPULib

A couple of users recently asked about doing a shift of a matrix using GPULib. While there is no interface that mimics 'shift' at the moment, all the functionality for shifting arrays is there. For example, you can use the following to shift an array along the x-direction:

IDL> x = findgen(5, 5)
IDL> print, x
0.00000 1.00000 2.00000 3.00000 4.00000
5.00000 6.00000 7.00000 8.00000 9.00000
10.0000 11.0000 12.0000 13.0000 14.0000
15.0000 16.0000 17.0000 18.0000 19.0000
20.0000 21.0000 22.0000 23.0000 24.0000

IDL> a = gpuPutArr(x)
IDL> b = gpuFltarr(5, 5)

; here comes the actual shift:
; we first perform b[0:3, *] = a[1:*, *]
; see the documentation for 'gpuSubArr' for an
; explanation of the arguments.
IDL> gpuSubArr, a, [1, -1], -1, b, [0, 3], -1

; and then b[4, * ] = a[0,*]
IDL> gpuSubArr, a, 0, -1, b, 4, -1

IDL> res = gpuGetArr(b)
IDL> print, res

1.00000 2.00000 3.00000 4.00000 0.00000
6.00000 7.00000 8.00000 9.00000 5.00000
11.0000 12.0000 13.0000 14.0000 10.0000
16.0000 17.0000 18.0000 19.0000 15.0000
21.0000 22.0000 23.0000 24.0000 20.0000

Shifts in other directions can be implemented similarly. This is not the fastest possible shift, but at least it should allow you to perform the shift on the GPU, rather than transferring it back to the CPU.

As always, if you have questions about this or any other use of GPULib, feel free to leave comments here, or email us at support@txcorp.com.

Wednesday, November 19, 2008

Error: "Attempt to call undefined procedure/function: 'GPUADD'."

Some users have seen the following error when using GPULib's IDL bindings:

IDL> gpuinit
% Attempt to call undefined procedure/function: 'GPUADD'.
% Error occurred at: GPUINIT 137 C:\Documents and
Settings\Sizemore\Desktop\gpulib-1.0.0\IDL\lib\gpuinit.pro
% $MAIN$
% Execution halted at: $MAIN$

If this happens to you, you may need to copy the file IDL\lib\gpuinit.pro.orig to IDL\lib\gpuinit.pro, overwriting the original gpuinit.pro file.

UPDATE: We've released a new version of GPULib, v1.0.2, that fixes this error. See

http://gpulib.txcorp.com

to get the new version.

Wednesday, August 13, 2008

IDL tip: compile gpuinit before running your own procedures

Another "gotcha" came up today from a user; when using GPULib in IDL, you may need to compile gpuinit before running your code. For example, to run a procedure you've written called foo, you should do something like the following:
IDL> .compile gpuinit
IDL> .compile foo
IDL> foo
As always, if you have questions or suggestion regarding GPULib, feel free to leave comments or email support@txcorp.com.