14Jul/10Off

## Using CULA in MATLAB, Part 1

by Kyle

It is possible to use CULA to accelerate a large number of popular linear algebra routines in MATLAB. With some relatively simple interface code, it's easy to seamlessly accelerate popular MATLAB routines such as: qr, lu, svd, eig, and mldivide.

When writing your own MATLAB code that uses an external library such as CULA, there are a few common hurdles to overcome. In the next few blog posts, we'll introduce these problems and discuss how to solve them. By the end of the series, we'll have enough code to implement a CULA accelerated Singular Value Decomposition (SVD), culasvd, that matches the functionality of MATLAB's svd routine.

Part 1 of this series will explain how to handle multiple CULA data types in MATLAB. Part 2 will detail some tips and tricks on interface wrapping. Finally, Part 3 will show how to compile, link, and benchmark CULA within MATLAB.

Excluding sparse, there are four main data types supported by MATLAB and CULA: single, double, single complex, and double complex. When writing a CULA accelerated MATLAB routine, your function should be able to support all of these types invisibly to the end user. Since both the MATLAB MEX compiler and CULA support C++ templates, we can easily implement all four datatypes through a single code path.

The following code snippet shows how to use C++ templates to handle all four MATLAB types. We use the MATLAB functions, mxGetClassID() and mxIsComplex(), to determine the precision and complexity of the input matrix and then switch into an explicitly typed CULA wrapper.

```// culaSvd.cpp
// Implements CULA accelerated SVD to be called within MATLAB

#include "mex.h"
#include "culapack.hpp"

// Templated wrapper
template <typename T>
void mexCulaGesvd(int nlhs,           /* number of expected outputs */
mxArray* plhs[],        /* output pointer array */
int nrhs,               /* number of inputs */
const mxArray* prhs[],  /* input pointer array */
mxClassID id,
mxComplexity complexity)
{
// Function core, details in "Using CULA in MATLAB, Part 2"
// culaGesvd(...)
}

// MATLAB Gateway Function
void mexFunction(int nlhs,              /* number of expected outputs */
mxArray* plhs[],       /* output pointer array */
int nrhs,              /* number of inputs */
const mxArray* prhs[]  /* input pointer array */ )
{
// We only support a full SVD in this example
if(nrhs != 1)
mexErrMsgTxt("culasvd: Must have 1 input argument [X]");
if(nlhs != 3)
mexErrMsgTxt("culasvd: Must have 3 output arguments [U,S,V]");

// Get precision (single or double)
mxClassID classID = mxGetClassID(prhs[0]);

// Get complexity (real or complex)
mxComplexity complexity = mxIsComplex(prhs[0]) ? mxCOMPLEX : mxREAL;

// Switch based on data type
if ( classID == mxSINGLE_CLASS && complexity == mxREAL)
mexCulaGesvd<culaFloat>(nlhs, plhs, nrhs, prhs, classID, complexity);
else if (classID == mxDOUBLE_CLASS && complexity == mxREAL )
mexCulaGesvd<culaDouble>(nlhs, plhs, nrhs, prhs, classID, complexity);
else if ( classID == mxSINGLE_CLASS && complexity == mxCOMPLEX )
mexCulaGesvd<culaFloatComplex>(nlhs, plhs, nrhs, prhs, classID, complexity);
else if ( classID == mxDOUBLE_CLASS && complexity == mxCOMPLEX )
mexCulaGesvd<culaDoubleComplex>(nlhs, plhs, nrhs, prhs, classID, complexity);
else
mexErrMsgTxt("culasvd: Unknown or unsupported data type");
}
```

As a quick overview, this file includes the 3 main components:

1) The CULA header for our SVD function and the MATLAB header for the MATLAB data types and API
2) A templated function that will support all four MATLAB data types
3) The MATLAB gateway function, 'mexFunction', required by all C/C++ programs called by MATLAB

These are the three main components needed by any MATLAB wrapper for CULA. As you can see, we really haven't gotten into any SVD specific code with the exception of the error checking.

Check back soon for Part 2 where we'll go into more details about how to wrap MATLAB's matrix data into CULA's SVD function. If you have any questions, feel free to ask in a comment here or on our forums.