GPULib 1.8 released

GPULib 1.8 has been released with updates to the underlying libraries as well as many other features in many areas of the library. It has been updated to use the most recent versions of IDL and CUDA, IDL 8.4 and CUDA 6.5. The new features are:

• Support for integer data types. I have been wanting to support integer types in GPULib for awhile and now GPULib supports all the numeric types provided by IDL! We can finally do:

  dx = gpuIndgen(10)
  

• Added GPUREPMAT routine. This is a handy routine to create a new array by repeating a 2-dimensional array in a grid.

• Added GPUCREATEKERNEL routine to create the source code of a simple kernel. This is a code generation routine that can be loaded with GPULOADMODULE/GPULOADFUNCTION and executed with GPUEXECUTEFUNCTION.

• Added GPUFINITE routine similar to IDL's library routine.

• Added linear algebra routines GPULUDC, GPULUSOL, and GPULEAST_SQUARES. This fills out more of the GPU equivalent of the convenience routines provided by IDL so that the LAPACK interface of MAGMA is not required to perform linear algebra computations.

• Added support for RHO and THETA keywords in GPURADON.

• Added GPUMEAN routine. This routine has DIMENSION and NAN keywords with the same functionality as IDL's library routine.

GPULib 1.6.2 released

GPULib 1.6.2 has been released today! Updates include:

• support for CUDA 5.5
• updated to MAGMA 1.4
• improved syntax in GPUEXECUTEFUNCTION
• self-extracting executable installer on Linux and package installer on OS X
• small bug fixes

We intend to keep updating GPULib to use the newest versions of IDL and CUDA as those packages are released. If you need a different combination of IDL and CUDA, please let us know.

The update to MAGMA 1.4 should fix some of the issues that some were experiencing with GPUINVERT at particular matrix sizes.

For the most part, we updated the libraries that GPULib is built on in this release and not added new features, but there is one new convenience in using GPUEXECUTEFUNCTION in GPULb 1.6.2. After reading PTX code and loading a module and function, the CUDA kernel can be launched with the nicer syntax:

  gpuExecuteFunction, kernel, $
                      dx, $
                      dy, $
                      dz, $
                      n, $
                      GRID=[nBlocksPerGrid], $
                      BLOCKS=[nThreadsPerBlock], $
                      ERROR=err

instead of the old, cumbersome syntax of calling the _getHandle method for each GPU variable:

  gpuExecuteFunction, kernel, $
                      dx->_getHandle(), $
                      dy->_getHandle(), $
                      dz->_getHandle(), $
                      n, $
                      GRID=[nBlocksPerGrid], $
                      BLOCKS=[nThreadsPerBlock], $
                      ERROR=err

It should be possible to get the full functionality of GPULib without calling any of the private methods of the GPUVariable class now.

GPULib VISualize talk

My colleague Jon Rood presented on the recently released GPULib 1.6 at VISualize last week. Here are the slides for the talk.

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.

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.

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!

AGU 2012 poster

Here is the abstract for my poster at AGU this year, “GPU accelerated curve fitting with IDL”:

Curve fitting is a common mathematical calculation done in all scientific areas. The Interactive Data Language (IDL) is also widely used in this community for data analysis and visualization. We are creating a general-purpose, GPU accelerated curve fitting library for use from within IDL.

We have developed GPULib, a library of routines in IDL for accelerating common scientific operations including arithmetic, FFTs, interpolation, and others. These routines are accelerated using modern GPUs using NVIDIA’s CUDA architecture. We will add curve fitting routines to the GPULib library suite, making curve fitting much faster.

In addition, library routines required for efficient curve fitting will also be generally useful to other users of GPULib. In particular, a GPU accelerated LAPACK implementation such as MAGMA is required for the Levenberg- Marquardt curve fitting and is commonly used in many other scientific computations. Furthermore, the ability to evaluate custom expressions at runtime necessary for specifying a function model will be useful for users in all areas.

The poster is available to download on the AGU website.

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.

Multiplying a scalar by an array

Currently, you have to be a bit tricky to multiply a scalar by an array. (In the next version of GPULib, there will be a efficient, direct way to do scalar/array operations.)

Most GPULib operations have a short form and a long form. The short form is a simple version of the operation that takes two array arguments. For addition, the short form is equivalent to:

C = A + B

where A, B, and C are all arrays of the same size. In code, this is done with:

c = gpuadd(a, b)

The long form takes five arguments, three scalars and two arrays, in the form of:

C = s1*A + s2*B + s3

where A, B, and C are arrays of the same size and s1, s2, and s3 are scalars. In code, this is:

c = gpuadd(s1, a, s2, b, s3)

If you just want to do a scalar times an array, use the long form with A for both A and B, while also setting b and c to 0.0. For example, to do 2. * findgen(10), do:

IDL> gpuinit
Welcome to GPULib 1.5.0 (Revision: 2199M)
Graphics card: Tesla C2070, compute capability: 2.0, memory: 751 MB available, 1279 MB total
Checking GPU memory allocation...no errors
IDL> a = findgen(10)
IDL> da = gpuputarr(a)
IDL> dc = gpufltarr(10)
IDL> scalar = 2.0
IDL> ; dc = scalar * da + 0. * da + 0.
IDL> dc = gpuadd(scalar, da, 0., da, 0., lhs=dc) 
IDL> c = gpugetarr(dc)
IDL> print, c
      0.00000      2.00000      4.00000      6.00000      8.00000
      10.0000      12.0000      14.0000      16.0000      18.0000

Note, that most operations have a long form version that performs an affine transform along with the operation.

GPULib 1.4.4 released!

We are pleased to announce that GPULib 1.4.4 is available from the Tech-X website.


The following changes have been made for version 1.4.4:

• Support for CUDA 3.2 (rather than CUDA 3.2 rc2)
• Added more extensive checks for issues in GPUINIT.
• Added ability to transpose 1-dimensional arrays like TRANSPOSE does.
• Now handles ranges, index arrays, and * in any indexing position of an array in the same manner as IDL.
• Now can use single-dimensional indexing for multi-dimensional arrays.
• Trying to create a view from discontinuous memory is now disallowed.
• Fixed garbage collection issue with repeated assignment to the same variable name.
• Fixed an issue whereby multiple calls to GPUINIT were causing cudaErrorSetOnActiveProcess errors.
• Fixed memory allocation bug in GPUTRANSPOSE.
• Fixed typo in GPUMATRIX_MULTIPLY.
• Added unit tests for GPUTRANSPOSE, GPUMAXOP.
• Now prints an error message when trying to create a view from discontinuous memory.
• Fixed build issue on Mac OS X Leopard.
• Fixed several install issues for Windows when building from source.
• Improvements to the QuickStartGuide for Windows.
• Fixed getDimensions bug for certain views.
• Added error checking for changing dimensions.
• Fixed memory leak in GPUFFT.
• Fixed errors in GPUMAXOP.
• Fixed failing unit tests of GPUMAXOP.
• Added missing demo file to install list.
• Fixed digital signing of installer file. 
• Fixed disappearing graphic in transform3d demo.
• Corrected IDL_LIBDIR path.
• Fixed path of file gpulib_version.txt in CMakeLists.txt.
• Fixed full reduction in GPUTOTAL.
• Fixed -fPIC warnings when building from source on Windows.
• Fixed color table issue in GPU_SWIRL_DEMO.
• Added "press q to stop" option to GPU_FDTD_DEMO.
• Fixed flipped image in GPU_DECONHUBBLE_DEMO.

GPULib 1.4.2 released!

We are pleased to announce that GPULib 1.4.2 is available from the Tech-X website.

For Windows users, the installer contains the following pre-built versions of GPULib 1.4.2:

1. 32- or 64-bit single precision (compute capability 1.0) built with CUDA 3.2
2. 32- or 64-bit double precision (compute capability 1.3) built with CUDA 3.2
3. 32- or 64-bit Fermi (compute capability 2.0) built with CUDA 3.2

(As before, Windows users can opt to build GPULib from source if they so desire).

The following changes have been made for version 1.4.2:

• Fixed make install so that appropriate files are now installed.
• All demos cleaned up and numerous memory leaks fixed.
• Created a new comprehensive PDF for Windows users, covering general use as well as how to build from source.
• Windows installation now includes all source files in a zip archive, so that users who are not building from source do not have their install folder cluttered with source files.
• gpuinit() now displays amount of free memory as well as total memory available on the device.
• Fixed bug in GPUMATRIX_MULTIPLY. 
• Fixed several small bugs and memory leaks.

We are grateful to Mort Canty and David Grier for submitting bug reports which led to fixes for this release!

GPULib 1.4.0 released!

We are pleased to announce that GPULib 1.4.0 is available from the Tech-X website.

For Windows users, we created an installer which installs one of the following pre-built versions of GPULib 1.4.0:

1. 32- or 64-bit single precision (compute capability 1.0) built with CUDA 3.1
2. 32- or 64-bit double precision (compute capability 1.3) built with CUDA 3.1
3. 32- or 64-bit Fermi (compute capability 2.0) built with CUDA 3.1
4. 32-bit single precision (compute capability 1.0) built with CUDA 2.3
5. 32-bit double precision (compute capability 1.3) built with CUDA 2.3

(As before, Windows users can opt to build GPULib from source.)

In addition, the following changes have been made:

• Now builds with CMake for cross-platform compatibility.
• Now supports CUDA streams, enabling concurrent execution of multiple kernels.
• Now supports asynchronous data transfer. 
• Now leverages new features of IDL 8.0 enabling more seamless integration between the two products. 
• Includes a variety of new algorithms, such as functions for sorting and large histogramming. 
• gpuinit() now provides additional information, including GPULib version and model of your CUDA-enabled card, and additionally performs a small allocation to ensure the card is available.
• Demos updated to fail gracefully if you do not have enough memory on your card, or if your card does not have the compute capability to run the demo.
• Fixed a bug whereby assigning a regular IDL variable to a slice failed.
• gpushift() fixed for the case where a dimension is shifted by 0.
• Fixed problems with gpuInterpolate functions.
• Fixed memory leak in gpufix().
• Fixed memory leak in gpufft().
• MATLAB support has been dropped.

GPULib 1.2.2 Released

GPULib 1.2.2 is now available from the Tech-X website. 

We've made extensive changes to the build system, which is now cleaner and more robust. Full release notes follow.

Build system changes 
- Added --with-extra-nvcc-flags=... to configure which allows extra flags to be passed to nvcc.
- If --prefix is not set, make install will fail gracefully, instead of attempting to install in /.
- Fixed --with-matlab-dir=... configure option.
- Added IDL and MATLAB configuration info to config.summary to make it easier to troubleshoot problem.
- Added several missing Windows build files.
- Removed several obsolete files and directories. 
- Running 'make clean' will not affect documentation.
- Running 'make install' will properly build code if not already built.
- Install directory is now laid out properly.
- Fixed "No rule to make target `docs/GPULib_UsersGuide.pdf', needed by `all-am'" error.

IDL  changes  
- Fixed bug whereby FFT was only operating on a single row.
- Fixed bug whereby GPUPOW was not found.
- Corrected bwtest example.
- Fixed time reporting for FDTD demo.
- Fixed typo in FDTD demo README.

MATLAB changes 
- Added potential to specify the device number form gpuInit() and accInit().
- Fixed Bug in gpuSet function.
- Fixed Makefile which was incorrect for 32-bit Linux.

2nd edition of Mort Canty's book uses GPULib

From the CRC Press site for Image Analysis, Classification, and Change Detection in Remote Sensing: With Algorithms for ENVI/IDL, Second Edition:

This popular introduction to the processing of remote sensing imagery has been updated to include coverage of the latest versions of the ENVI software environment. This new edition covers support vector machines and other kernel-based methods. Illustrating many programming examples in the array-oriented language ID, the text includes coverage of basic Fourier, wavelet, principal components and minimum noise fraction transformations; convolution filters, topographic modeling, image-to-image registration and ortho-rectification; image fusion; supervised and unsupervised land cover classification with neural networks; hyperspectral analysis; multivariate change detection.

I was excited to hear that GPULib was used in this version of the book. Mort says:

In the text I discuss routines for nonlinear principal component analysis, supervised classification and nonlinear clustering, and explain that they can take advantage of GPULib/CUDA, if installed. (I use your routine GPU_DETECT() to check for GPULib).

GPULib docs from ENVI menu

A recent post by Mort Canty provides a handy program that adds an item to ENVI's help menu that will bring up the GPULib docs.

GPULib 1.2 released

GPULib 1.2 is available from the Tech-X website. This release focused on improved MATLAB bindings with a few important bug fixes for the IDL bindings along with a few new kernels. Full release notes follow:

Changes/new features in GPULib version 1.2

General

The main focus of this release is on the improved MATLAB bindings. Some new kernels were added since the release of version 1.0.8.

GPULib kernels

gpuAtan2, gpuFmod, gpuPow

IDL bindings

- Support for the new kernels. For the time being, these funtions only support float and double (so no complex types) and no affine transform arguments.
- Added example bwtest.pro showing the use of page-locked variables for fast CPU/GPU data transfer.
- Added finite-different time-domain example demonstrating the use of views for efficient array sub-selection
- Added spectral angle mapper example.
- Bug fixes for decon_hubble example
- improved documentation

MATLAB bindings

MATLAB GPULib version 1.2 has many major changes from the previous release. READ the README!

First and foremost, there are two distinct and completely separate interfaces to the library. They should NEVER be intermingled.

1.) The accArray class replaces the old gpuArray class from the previous release. This interface requires MATLAB R2008a or higher. This interface hass automatic garbage collection, overloaded operators, and overloaded versions of native MATLAB functions, ...
2.) "gpu" interface class can be used with older versions of MATLAB though it's not clear how far back one can go.

The interface was redesigned for speed. The accArray class is about 2.5X faster than the gpuArray interface for many functions tested. Some of the "gpu"-prefixed functions can be up to 10X faster than the gpuArray interface.

MATLAB GPULib has many new functions including
1.) fft, ifft, fft2, and ifft2
2.) Reduction operations, including sum, cumsum, prod, cumprod,... These functions support 1D vectors and 2D matrices currently.
3.) Single (Complex) and Double (Complex) precision versions of Matrix Multiplication, Transpose and Complex Conjugate Transpose.

The accArray class does not support the subsref.m (i.e. b=A(i)), subsasgn.m (i.e. A(i)=b), or array concatentation functions like A=[B; C; D], These will be supported in future releases.

The "gpu" interface supports subscripting through the gpuSubsref, gpuSubsasgn, and gpuSub2ind functions.

Both interfaces support page-locked host memory allocation via cudaMallocHost. This gives the possibility of much faster memory transfer from CPU memory to GPU memory and back.

Both interfaces include more comprehensive native MATLAB-like documentation.

New examples include: bench, bwtest, fdtd, and fftExample.

GPULib slides from VISualize 2009

Here are the slides from Peter Messmer's GPULib talk at VISualize 2009 a couple weeks ago. Peter and I also stayed around after the IDL presentations to do hands on training for GPULib. The training was well attended; it was good to see the interest in GPU computing with IDL. Most people did not bring a laptop, making it less “hands on” than we had originally intended, but it was good to be able to answer individual questions.

NVIDIA Releases OpenCL Driver

OpenCL is officially out -- have a look at the press release, or the OpenCL site.

Detecting presense of GPULib

Over the weekend, a user (Mort Canty) asked the GPULib user's mailing list for a way to detect the presence. One of our developers wrote a routine to do just that, and I thought it might be of use to other GPULib users as well.

Below is the routine GPU_DETECT that will return a boolean indicating if GPULib is present. It will also initialize the GPU if it finds GPUINIT. If you have any questions, feel free to ask!

; docformat = 'rst'

;+
; Initializes GPULib if it is present.
;
; :Returns:
; 1 if GPULib is present, 0 if not
;
; :Params:
; devId: in, optional, type=numtype
; id of the GPU device to be used for GPU computations
;
; :Keywords:
; _ref_extra : in, optional, type=keywords
; keywords to GPUINIT
;-
function gpu_detect, devId, _ref_extra=e
compile_opt strictarr

catch, error
if (error ne 0L) then begin
catch, /cancel
return, 0
endif

gpuinit, devId, _extra

return, 1
end

GPULib 1.0.8 update available

GPULib 1.0.8 is available at

http://gpulib.txcorp.com

This release is addressing several bugs pointed out by users; we're still working on the MATLAB refactoring for a future version. As always, if you have any questions, feel free to post in the comments, the GPULib mailing list, or contact us directly at support@txcorp.com.

Also, one note; during testing of this release, we found some problems with using it on Windows Vista. We're looking for solutions, but at the moment XP is the better choice for running GPULib.