{"id":3898,"date":"2011-10-24T16:58:55","date_gmt":"2011-10-24T15:58:55","guid":{"rendered":"http:\/\/www.walkingrandomly.com\/?p=3898"},"modified":"2011-10-24T16:58:55","modified_gmt":"2011-10-24T15:58:55","slug":"matlab-mex-functions-using-the-nag-c-library-in-windows","status":"publish","type":"post","link":"https:\/\/walkingrandomly.com\/?p=3898","title":{"rendered":"MATLAB mex functions using the NAG C Library in Windows"},"content":{"rendered":"

The NAG C Library<\/a> is one of the largest commercial collections of numerical software currently available and I often find it very useful when writing MATLAB mex files<\/a>.\u00a0 “Why is that?<\/em>” I hear you ask.<\/p>\n

One of the main reasons for writing a mex file is to gain more speed over native MATLAB.\u00a0 However, one of the main problems with writing mex files is that you have to do it in a low level language such as Fortran or C and so you lose much of the ease of use of MATLAB.<\/p>\n

In particular, you lose straightforward access to most of the massive collections of MATLAB routines that you take for granted.\u00a0 Technically speaking that’s a lie because you could use the mex function mexCallMATLAB<\/strong> to call a MATLAB routine from within your mex file but then you’ll be paying a time overhead every time you go in and out of the mex interface.\u00a0 Since you are going down the mex route in order to gain speed, this doesn’t seem like the best idea in the world.\u00a0 This is also the reason why you’d use the NAG C Library and not the NAG Toolbox for MATLAB<\/a> when writing mex functions.<\/p>\n

One way out that I use often is to take advantage of the NAG C library and it turns out that it is extremely easy to add the NAG C library to your mex projects on Windows.\u00a0 Let’s look at a trivial example.\u00a0 The following code, nag_normcdf.c<\/a>, uses the NAG function nag_cumul_normal <\/strong>to produce a simplified version of MATLAB’s normcdf function (laziness is all that prevented me from implementing a full replacement).<\/p>\n

\/* A simplified version of normcdf that uses the NAG C library\r\n * Written to demonstrate how to compile MATLAB mex files that use the NAG C Library\r\n * Only returns a normcdf where mu=0 and sigma=1\r\n * October 2011 Michael Croucher\r\n * www.walkingrandomly.com\r\n *\/\r\n\r\n#include <math.h>\r\n#include \"mex.h\"\r\n#include \"nag.h\"\r\n#include \"nags.h\"\r\n\r\nvoid mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )\r\n{\r\n    double *in,*out;\r\n    int rows,cols,num_elements,i;\r\n\r\n    if(nrhs>1)\r\n    {\r\n        mexErrMsgIdAndTxt(\"NAG:BadArgs\",\"This simplified version of normcdf only takes 1 input argument\");\r\n    } \r\n\r\n    \/*Get pointers to input matrix*\/\r\n    in = mxGetPr(prhs[0]);\r\n    \/*Get rows and columns of input matrix*\/\r\n    rows = mxGetM(prhs[0]);\r\n    cols = mxGetN(prhs[0]);\r\n    num_elements = rows*cols;\r\n\r\n    \/* Create output matrix *\/\r\n    plhs[0] = mxCreateDoubleMatrix(rows, cols, mxREAL);\r\n    \/* Assign pointer to the output *\/\r\n    out = mxGetPr(plhs[0]);\r\n\r\n    for(i=0; i<num_elements; i++){\r\n          out[i] = nag_cumul_normal(in[i]);\r\n         }\r\n}<\/pre>\n
To compile this in MATLAB, just use the following command<\/p>\n
 mex nag_normcdf.c CLW6I09DA_nag.lib<\/pre>\n
If your system is set up the same as mine then the above should ‘just work’ (see System Information<\/strong> at the bottom of this post).\u00a0 The new function works just as you would expect it to<\/p>\n
>> format long\r\n>> format compact\r\n>> nag_normcdf(1)\r\nans =\r\n   0.841344746068543<\/pre>\n
Compare the result to normcdf from the statistics toolbox<\/p>\n
>> normcdf(1)\r\nans =\r\n   0.841344746068543<\/pre>\n
So far so good.  I could stop the post here since all I really wanted to do was say ‘The NAG C library is useful for MATLAB mex functions and it’s a doddle<\/a> to use – here’s a toy example and here’s the mex command to compile it’<\/strong><\/p>\n
However, out of curiosity, I looked to see if my toy version of normcdf was any faster than The Mathworks’ version.\u00a0 Let there be 10 million numbers:<\/p>\n
>> x=rand(1,10000000);<\/pre>\n
Let’s time how long it takes MATLAB to take the normcdf of those numbers<\/p>\n
>> tic;y=normcdf(x);toc\r\nElapsed time is 0.445883 seconds.\r\n>> tic;y=normcdf(x);toc\r\nElapsed time is 0.405764 seconds.\r\n>> tic;y=normcdf(x);toc\r\nElapsed time is 0.366708 seconds.\r\n>> tic;y=normcdf(x);toc\r\nElapsed time is 0.409375 seconds.<\/pre>\n
Now let’s look at my toy-version that uses NAG.<\/p>\n
>> tic;y=nag_normcdf(x);toc\r\nElapsed time is 0.544642 seconds.\r\n>> tic;y=nag_normcdf(x);toc\r\nElapsed time is 0.556883 seconds.\r\n>> tic;y=nag_normcdf(x);toc\r\nElapsed time is 0.553920 seconds.\r\n>> tic;y=nag_normcdf(x);toc\r\nElapsed time is 0.540510 seconds.<\/pre>\n
So my version is slower!\u00a0 Never mind, I’ll just make my version parallel using OpenMP – Here is the code: nag_normcdf_openmp.c<\/a><\/p>\n
\/* A simplified version of normcdf that uses the NAG C library\r\n * Written to demonstrate how to compile MATLAB mex files that use the NAG C Library\r\n * Only returns a normcdf where mu=0 and sigma=1\r\n * October 2011 Michael Croucher\r\n * www.walkingrandomly.com\r\n *\/\r\n\r\n#include <math.h>\r\n#include \"mex.h\"\r\n#include \"nag.h\"\r\n#include \"nags.h\"\r\n#include <omp.h>\r\n\r\nvoid do_calculation(double in[],double out[],int num_elements)\r\n{\r\n    int i,tid;\r\n\r\n#pragma omp parallel for shared(in,out,num_elements) private(i,tid)\r\n    for(i=0; i<num_elements; i++){\r\n          out[i] = nag_cumul_normal(in[i]);\r\n         }\r\n}\r\n\r\nvoid mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )\r\n{\r\n    double *in,*out;\r\n    int rows,cols,num_elements;\r\n\r\n    if(nrhs>1)\r\n    {\r\n        mexErrMsgIdAndTxt(\"NAG_NORMCDF:BadArgs\",\"This simplified version of normcdf only takes 1 input argument\");\r\n    } \r\n\r\n    \/*Get pointers to input matrix*\/\r\n    in = mxGetPr(prhs[0]);\r\n    \/*Get rows and columns of input matrix*\/\r\n    rows = mxGetM(prhs[0]);\r\n    cols = mxGetN(prhs[0]);\r\n    num_elements = rows*cols;\r\n\r\n    \/* Create output matrix *\/\r\n    plhs[0] = mxCreateDoubleMatrix(rows, cols, mxREAL);\r\n    \/* Assign pointer to the output *\/\r\n    out = mxGetPr(plhs[0]);\r\n\r\n    do_calculation(in,out,num_elements);\r\n}<\/pre>\n
Compile that with<\/p>\n
 mex  COMPFLAGS=\"$COMPFLAGS \/openmp\"  nag_normcdf_openmp.c CLW6I09DA_nag.lib<\/pre>\n
and on my quad-core machine I get the following timings<\/p>\n
>> tic;y=nag_normcdf_openmp(x);toc\r\nElapsed time is 0.237925 seconds.\r\n>> tic;y=nag_normcdf_openmp(x);toc\r\nElapsed time is 0.197531 seconds.\r\n>> tic;y=nag_normcdf_openmp(x);toc\r\nElapsed time is 0.206511 seconds.\r\n>> tic;y=nag_normcdf_openmp(x);toc\r\nElapsed time is 0.211416 seconds.<\/pre>\n
This is faster than MATLAB and so normal service is resumed :)<\/p>\n
System Information<\/strong><\/p>\n