Modern CPUs are capable of parallel processing at multiple levels with the most obvious being the fact that a typical CPU contains multiple processor cores.\u00a0 My laptop, for example, contains a quad-core Intel Sandy Bridge i7 processor and so has 4 processor cores. You may be forgiven for thinking that, with 4 cores, my laptop can do up to 4 things simultaneously but life isn’t quite that simple.<\/p>\n
The first complication is hyper-threading<\/a> where each physical core appears to the operating system as two or more virtual cores.\u00a0 For example, the processor in my laptop is capable of using hyper-threading and so I have access to up to 8 virtual cores!\u00a0 I have heard stories where unscrupulous sales people have passed off a 4 core CPU with hyperthreading as being as good as an 8 core CPU…. after all, if you fire up the Windows Task Manager you can see 8 cores and so there you have it!\u00a0 However, this is very far from the truth since what you really have is 4 real cores with 4 brain damaged cousins<\/a>.\u00a0 Sometimes the brain damaged cousins can do something useful but they are no substitute for physical cores.\u00a0 There is a great explanation of this technology at makeuseof.com<\/a>.<\/p>\n The second complication is the fact that each physical processor core<\/strong> contains a SIMD<\/a> (Single Instruction Multiple Data) lane of a certain width. SIMD lanes, aka\u00a0 SIMD units or vector units, can process several numbers simultaneously with a single instruction rather than only one a time.\u00a0 The 256-bit wide SIMD lanes on my laptop’s processor, for example, can operate on up to 8 single (or 4 double) precision numbers per instruction.\u00a0 Since each physical core has its own SIMD lane this means that a 4 core processor could theoretically operate on up to 32 single precision (or 16 double precision) numbers per clock cycle!<\/p>\n So, all we need now is a way of programming for these SIMD lanes!<\/p>\n Intel’s SPMD Program Compiler<\/a>, ispc, is a free product that allows programmers to take direct advantage of the SIMD lanes in modern CPUS using a C-like syntax.\u00a0 The speed-ups compared to single-threaded code can be impressive with Intel reporting up to 32 times speed-up<\/a> (on an i7 quad-core) for a single precision<\/strong> Black-Scholes option pricing routine for example.<\/p>\n Using ispc on MATLAB Since ispc routines are callable from C, it stands to reason that we’ll be able to call them from MATLAB using mex.\u00a0 To demonstrate this, I thought that I’d write a sqrt function that works faster than MATLAB’s built-in version.\u00a0 This is a tall order since the sqrt function is pretty fast and is already multi-threaded<\/a>.\u00a0 Taking the square root of 200 million random numbers doesn’t take very long in MATLAB:<\/p>\n This might not be the most useful example in the world but I wanted to focus on how to get ispc to work from within MATLAB rather than worrying about the details of a more interesting example.<\/p>\n Step 1 – A reference single-threaded mex file<\/strong><\/p>\n Before getting all fancy, let’s write a nice, straightforward single-threaded mex file in C and see how fast that goes.<\/p>\n Save the above to a text file called sqrt_mex.c<\/strong> and compile using the following command in MATLAB<\/p>\n Let’s check out its speed:<\/p>\n Well, it works but it’s quite a but slower than the built-in MATLAB function so we still have some work to do.<\/p>\n Step 2 – Using the SIMD lane on one core via ispc<\/strong><\/p>\n Using ispc is a two step process.\u00a0 First of all you need the .ispc program<\/p>\n Save this to a file called ispc_sqrt.ispc<\/strong> and compile it at the Bash prompt<\/strong> using<\/p>\n This creates an object file, ispc_sqrt.o<\/strong>, and a header file, ispc_sqrt.h<\/strong>. Now create the mex file in MATLAB<\/p>\n Call this ispc_sqrt_mex.cpp <\/strong> and compile in MATLAB<\/strong> with the command<\/p>\n Let’s see how that does for speed:<\/p>\n So, we’ve improved on the single-threaded mex file a bit (1.37 instead of 2 seconds) but it’s still not enough to beat the MATLAB built-in.\u00a0 To do that, we are going to have to use the SIMD lanes on all 4 cores simultaneously.<\/p>\n Step 3 – A reference multi-threaded mex file using OpenMP<\/strong><\/p>\n Let’s step away from ispc for a while and see how we do with something we’ve seen before– a mex file using OpenMP (see here<\/a> and here<\/a> for previous articles on this topic).<\/p>\n Save this to a text file called openmp_sqrt_mex.c<\/strong> and compile in MATLAB<\/strong> by doing<\/p>\n Let’s see how that does (OMP_NUM_THREADS has been set to 4):<\/p>\n That’s very similar to the MATLAB built-in and I suspect that The Mathworks have implemented their sqrt function in a very similar manner. Theirs will have error checking, complex number handling and what-not but it probably comes down to a for-loop that’s been parallelized using Open-MP.<\/p>\n Step 4 – Using the SIMD lanes on all cores via ispc<\/strong><\/p>\n To get a ispc program to run on all of my processors cores simultaneously, I need to break the calculation down into a series of tasks. The .ispc file is as follows<\/p>\n Compile this by doing the following at the Bash prompt<\/strong><\/p>\n We’ll need to make use of a task scheduling system. The ispc documentation suggests that you could use the scheduler in Intel’s Threading Building Blocks<\/a> or Microsoft’s Concurrency Runtime<\/a> but a basic scheduler is provided with ispc in the form of tasksys.cpp (I’ve also included it in the .tar.gz file in the downloads section at the end of this post), We’ll need to compile this too so do the following at the Bash prompt<\/strong><\/p>\n Finally, we write the mex file<\/p>\n In the above, the input array is divided into tasks where each task takes care of 1 million elements. Our 200 million element test array will, therefore, be split into 200 tasks– many more than I have processor cores. I’ll let the task scheduler worry about how to schedule these tasks efficiently across the cores in my machine. Compile this in MATLAB<\/strong> by doing<\/p>\n Now for crunch time:<\/p>\n There we have it! A version of the sqrt function that works faster than MATLAB’s own by virtue of the fact that I am now making full use of the SIMD lanes in my laptop’s Sandy Bridge i7 processor thanks to ispc.<\/p>\n Although this example isn’t very useful as it stands, I hope that it shows that using the ispc compiler from within MATLAB isn’t as hard as you might think and is yet another tool in the arsenal of weaponry that can be used to make MATLAB faster.<\/p>\n Final Timings, downloads and links Finally, here’s some links and downloads<\/p>\n System Specs<\/strong><\/p>\n Modern CPUs are capable of parallel processing at multiple levels with the most obvious being the fact that a typical CPU contains multiple processor cores.\u00a0 My laptop, for example, contains a quad-core Intel Sandy Bridge i7 processor and so has 4 processor cores. You may be forgiven for thinking that, with 4 cores, my laptop […]<\/p>\n","protected":false},"author":3,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"jetpack_post_was_ever_published":false,"footnotes":"","jetpack_publicize_message":"","jetpack_publicize_feature_enabled":true,"jetpack_social_post_already_shared":false,"jetpack_social_options":{"image_generator_settings":{"template":"highway","default_image_id":0,"font":"","enabled":false},"version":2}},"categories":[53,4,11,41,7],"tags":[],"class_list":["post-3988","post","type-post","status-publish","format-standard","hentry","category-making-matlab-faster","category-math-software","category-matlab","category-parallel-programming","category-programming"],"jetpack_publicize_connections":[],"jetpack_featured_media_url":"","jetpack_shortlink":"https:\/\/wp.me\/p3swhs-12k","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/3988","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/users\/3"}],"replies":[{"embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=3988"}],"version-history":[{"count":27,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/3988\/revisions"}],"predecessor-version":[{"id":4035,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/3988\/revisions\/4035"}],"wp:attachment":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=3988"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=3988"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=3988"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}
\n<\/strong><\/p>\n>> x=rand(1,200000000)*10;\r\n>> tic;y=sqrt(x);toc\r\nElapsed time is 0.666847 seconds.<\/pre>\n
#include <math.h>\r\n#include \"mex.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 \/*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 {\r\n out[i] = sqrt(in[i]);\r\n }\r\n}<\/pre>\n mex sqrt_mex.c<\/pre>\n
>> x=rand(1,200000000)*10;\r\n>> tic;y=sqrt_mex(x);toc\r\nElapsed time is 1.993684 seconds.<\/pre>\n
export void ispc_sqrt(uniform double vin[], uniform double vout[],\r\n uniform int count) {\r\n foreach (index = 0 ... count) {\r\n vout[index] = sqrt(vin[index]);\r\n }\r\n}<\/pre>\nispc -O2 ispc_sqrt.ispc -o ispc_sqrt.o -h ispc_sqrt.h --pic<\/pre>\n
#include <math.h>\r\n#include \"mex.h\"\r\n#include \"ispc_sqrt.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 \/*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 ispc::ispc_sqrt(in,out,num_elements);\r\n}<\/pre>\n mex ispc_sqrt_mex.cpp ispc_sqrt.o<\/pre>\n
>> tic;y=ispc_sqrt_mex(x);toc\r\nElapsed time is 1.379214 seconds.<\/pre>\n
#include <math.h>\r\n#include \"mex.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;\r\n\r\n#pragma omp parallel for shared(in,out,num_elements)\r\n for(i=0; i<num_elements; i++){\r\n out[i] = sqrt(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,i; \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 mex openmp_sqrt_mex.c CFLAGS=\"\\$CFLAGS -fopenmp\" LDFLAGS=\"\\$LDFLAGS -fopenmp\"<\/pre>\n
>> tic;y=openmp_sqrt_mex(x);toc\r\nElapsed time is 0.641203 seconds.<\/pre>\n
task void\r\nispc_sqrt_block(uniform double vin[], uniform double vout[],\r\n uniform int block_size,uniform int num_elems){\r\n uniform int index_start = taskIndex * block_size;\r\n uniform int index_end = min((taskIndex+1) * block_size, (unsigned int)num_elems);\r\n\r\n foreach (yi = index_start ... index_end) {\r\n vout[yi] = sqrt(vin[yi]);\r\n }\r\n}\r\n\r\nexport void\r\nispc_sqrt_task(uniform double vin[], uniform double vout[],\r\n uniform int block_size,uniform int num_elems,uniform int num_tasks)\r\n{\r\n\r\n launch[num_tasks] < ispc_sqrt_block(vin, vout, block_size, num_elems) >;\r\n}\r\n<\/pre>\nispc -O2 ispc_sqrt_task.ispc -o ispc_sqrt_task.o -h ispc_sqrt_task.h --pic<\/pre>\n
g++ tasksys.cpp -O3 -Wall -m64 -c -o tasksys.o -fPIC<\/pre>\n
#include <math.h>\r\n#include \"mex.h\"\r\n#include \"ispc_sqrt_task.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,i;\r\n\r\n unsigned int num_elements;\r\n unsigned int block_size;\r\n unsigned int num_tasks; \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 block_size = 1000000;\r\n num_tasks = num_elements\/block_size;\r\n\r\n ispc::ispc_sqrt_task(in,out,block_size,num_elements,num_tasks);\r\n\r\n}<\/pre>\nmex ispc_sqrt_task_mex.cpp ispc_sqrt_task.o tasksys.o<\/pre>\n
>> x=rand(1,200000000)*10;\r\n>> tic;ys=sqrt(x);toc %MATLAB's built-in\r\nElapsed time is 0.670766 seconds.\r\n>> tic;y=ispc_sqrt_task_mex(x);toc %my version using ispc\r\nElapsed time is 0.393870 seconds.<\/pre>\n
\n<\/strong><\/p>\n\n
\n
\n