I am looking for information from anyone that has tried to use OpenCL in their scientific code. Has anyone tried (recently) ViennaCL? If so, how does it compare to cusp?
What about OCLTools? Does it live up to the promise? If so, would it be a feasible way to start writing math kernels in OpenCL?
First of all I wish to thanks Aron Ahmadia for pointing me to this thread.
As for OpenCL in scientific code: OpenCL is meant to be a low-level API, thus it is crucial to wrap this functionality in some way in order to reach a reasonable productivity. Moreover, as soon as several compute kernels are involved, code can get VERY dirty if OpenCL kernel and memory handles need to be heavily passed around within an application. I don't know OCLTools, thus I can't say whether they are useful in this regard.
As for ViennaCL: I'm the head of ViennaCL, so I've worked recently with the library. :-) In the following I'll treat the request for a comparison with cusp in a slightly larger scope, namely ViennaCL versus the CUDA-based math libraries cusp and MAGMA. Only the present state is considered, even though there is a lot of ongoing development (at least on our side).
Functionality. MAGMA provides BLAS-functionality for dense matrices via the usual function interfaces. Most of this functionality is also provided with ViennaCL 1.2.0 using operator overloads and other syntactic sugar.
The same three iterative solvers (CG, BiCGStab, GMRES) are provided with cusp and ViennaCL. The set of preconditioners differs notably: Cusp provides diagonal, SA-AMG and various Bridson preconditioners. ViennaCL offers incomplete LU factorizations, diagonal preconditioners, and recently various AMG flavors and Sparse Approximate Inverse preconditioners. To my knowledge, all cusp preconditioners run entirely on the GPU, while ViennaCL relies particularly during the setup phase on CPU-based computations. Currently, the number of sparse matrix formats is larger in cusp: COO, CSR, DIA, ELL, HYB, while ViennaCL 1.2.0 provides COO and CSR.
There are a number of additional features provided with ViennaCL, which are not part of either MAGMA or cusp: Structured matrix types (Circulant, Hankel, etc.), fast Fourier transform, reordering algorithms (e.g. Cuthill-McKee) and wrappers for linear algebra types from other libraries.
Performance. The larger set of features and hardware support in ViennaCL typically comes at the cost of lower performance when compared to CUDA-based implementations. This is also partly due to the fact that CUDA is tailored to the architecture of NVIDIA products, while OpenCL represents in some sense a reasonable compromise between different many-core architectures.
Overall, ViennaCL is at present slower than MAGMA, particularly at BLAS level 3. The reasons is the different focus of ViennaCL (sparse instead of dense linear algebra) and thus the higher degree of optimization in MAGMA. Particularly BLAS level 3 operations are currently considerably faster in MAGMA.
Similarly, cusp provides slightly better overall performance in general. However, since sparse matrix operations are usually memory bandwidth limited, differences are considerably smaller and often negligible compared to data setup and the like. The choice of the preconditioner and its parameters usually has a higher impact on the overall execution time than any performance differences in sparse matrix-vector multiplications.
Portability. As for hardware portability, ViennaCL can use CPUs and GPUs from all major vendors thanks to OpenCL. In contrast, cusp and MAGMA rely on a suitable NVIDIA GPU.
ViennaCL is header-only, can be compiled on a wide range of C++ compilers and only needs to be linked with the OpenCL library if GPU-support is required. In principle, the generic algorithms in ViennaCL can also be used without any OpenCL linkage, while cusp and MAGMA require the NVIDIA compiler for compilation and the CUDA library on the target system for execution. MAGMA also requires a BLAS library, which can sometimes be a bit of a hassle to find or install on a new system.
API. MAGMA provides BLAS-style function interfaces for BLAS functionality. The C++ interface of cusp also uses some functions from BLAS, but no operator overloads. Since most interfaces in ViennaCL are intentionally similar to Boost.uBLAS and feature syntactic sugar such as operator overloads, ViennaCL is also intended to be used like Boost.uBLAS. Thus, in addition to just calling a predefined set of operations and algorithms, our intention is to make a transition from purely CPU-based execution to GPU code as simple as possible, even if non-standard algorithms are to be used. In the case that a dedicated OpenCL kernel is required, there is also a framework for integrating your own OpenCL kernels in ViennaCL. Thus, ViennaCL aims a lot more towards high productivity in the sense that the time required for implementing new algorithms on the GPU is minimized. These savings can significantly outweigh any performance penalty (if any) compared to cusp and MAGMA. (It has also be mentioned in the thread on unit testing that code development time is a precious resource in science.)
There are certainly a number of ideological issues (e.g. CUDA vs. OpenCL, BLAS-interface vs. operator overloads) throughout my comparison, but their discussion is beyond the scope of the initial question.
OpenCL can be used, however, there is a lack of infrastracture, e.g. important mature standard math libraries with tuned de facto standard linear algebra components and to some extent good profiling tools, albeit the latter problem has improved significantly for GPUs. This is available in CUDA as of today and can be contributed to a part of Nvidia's success with this programming model. However, OpenCL seems to be catching up (slowly).
Today, as a starting point for GPU programming CUDA is fine, and if it is needed, there is nothing that prevents using OpenCL as an alternative, e.g. to make the code more portable. Essentially, the same kernel code can be used both in CUDA and OpenCL so it should not be a major problem to go from CUDA to OpenCL. This is confirmed by own experiences testing this out. From a performance perspective, the same optimization techniques can be used and for trivial concurrent code compilers should do the job fine (e.g. loop unrolling, etc.).
