File Exchange

## MMX - Multithreaded matrix operations on N-D matrices

version 1.2 (460 KB) by

N-D matrix operations. Like Boettcher's ndfun or Tursa's mtimesx, only faster.

4.91667
14 Ratings

Updated

MMX treats an N-D matrix of double precision values as a set of pages
of 2D matrices, and performs various matrix operations on those pages.
MMX uses multithreading over the higher dimensions to achieve good
performance. Full singleton expansion is available for most operations.

C = MMX('mult', A, B) is equivalent to the matlab loop
for i=1:N,
C(:,:,i) = A(:,:,i) * B(:,:,i);
end
Singleton expansion is enabled on all dimensions so for example if
A = randn(5,4,3,10,1);
B = randn(4,6,3,1 ,6);
C = zeros(5,6,3,10,6);
then C = mmx('mult',A,B) equivalent to
for i = 1:3
for j = 1:10
for k = 1:6
C(:,:,i,j,k) = A(:,:,i,j,1) * B(:,:,i,1,k);
end
end
end

C = MMX('mult', A, B, mod) and where mod is a modifier string, will
transpose one or both of A and B. Possible values for mod are
'tn', 'nt' and 'tt' where 't' stands for 'transposed' and 'n' for
'not-transposed'. For example
>> size(mmx('mult',randn(4,2),randn(4,2),'tn'))
ans = 2 2

C = MMX('square', A, []) will perform C = A*A'
C = MMX('square', A, [],'t') will perform C = A'*A

C = MMX('square', A, B) will perform C = 0.5*(A*B'+B*A')
C = MMX('square', A, B, 't') will perform C = 0.5*(A'*B+B'*A)

C = MMX('chol', A, []) will perform C = chol(A)

C = MMX('backslash', A, B) will perform C = A\B
Unlike other commands, 'backslash' does not support singleton
expansion. If A is square, mmx will use LU factorization, otherwise it
will use QR factorization. In the underdetermined case, (i.e. when
size(A,1) < size(A,2)), mmx will give the least-norm solution which
is equivalent to C = pinv(A)*B, unlike matlab's mldivide.

C = MMX('backslash', A, B, 'U') or MMX('backslash', A, B, 'L') will
perform C = A\B assuming that A is upper or lower triangular,
respectively.

C = MMX('backslash', A, B, 'P') will perform C = A\B assuming that A
is symmetric-positive-definite.

MMX(n) does thread control: mmx will automatically start a number of
threads equal to the number of available processors, however the
number can be set manually to n using the command mmx(n). mmx(0) will

IMPORTANT NOTE: The functions which assume special types of square
matrices as input ('chol' and 'backslash' for 'U','L' or 'P'
modifiers) do not check that the inputs are indeed what you say they
are, and produce no error if they are not. Caveat computator.

COMPILATION: To compile run 'build_mmx'. Type 'help build_mmx' to read

Michal Kvasnicka

### Michal Kvasnicka (view profile)

@Israel Vaughn ... I can confirm the fact, that MMX produce very distorted results, as Israel Vaughn mentioned here! Something wrong with optional MKL and BLAS versions of MMX.

Any idea?

Michal Kvasnicka

### Michal Kvasnicka (view profile)

The MMX package is in bad shape. The crucial function build_mmx.m is not internally consistent:

1. When you start build_mmx.m on Linux PC you always got the following error message:
Compilation of 'mmx_mkl_single' failed with error:
ls: cannot access '*small_mkl_ilp64.*': No such file or directory

Of course, the MKL instalation and other addiotional steps described in build_mmx.m was already done!!!

2. You must fix the line 137
to

After this you still need solve other problem, see
Trying to compile 'mmx_mkl_single', using
CXXFLAGS="\\$CXXFLAGS -I/opt/intel/mkl/include ", LDFLAGS="\\$LDFLAGS -L/usr/local/MATLAB/R2017a/extern/lib/glnxa64 ", -llibsingle_mkl_ilp64, -lpthread, -DUNIX_SYSTEM, -DUSE_BLAS, -DMKL_ILP64,
Building with 'g++'.
Compilation of 'mmx_mkl_single' failed with error:
/usr/bin/ld: cannot find -llibsingle_mkl_ilp64
collect2: error: ld returned 1 exit status

Could be anybody here so kind and help me with this problem.

Israel Vaughn

### Israel Vaughn (view profile)

Note that this mmxc implementation is correct for 'mult' other functions may compute incorrect answers.

Israel Vaughn

### Israel Vaughn (view profile)

This comment seems to have disappeared, anyway, for those asking about a complex version of mmx, here is a straightforward way to do it :

function [ C ] = mmxc( actionStr,A,B, transString )
% complex mmx
% not very optimal, but hey, it works

if (strcmp(actionStr,'square') & isempty(B))
At = real(A);
Ai = imag(A);
clear A;

C = mmx('square',At,[]);

C = C - 1i.*mmx('mult',At,Ai,'nt');

C = C + 1i.*mmx('mult',Ai,At,'nt');

C = C + mmx('square',Ai,[]);
else if (exist('transString'))
At = real(A);
Ai = imag(A);
clear A;
Bt = real(B);
Bi = imag(B);
clear B;

C = mmx(actionStr,At,Bt,transString);

C = C + 1i.*mmx(actionStr,At,Bi,transString);

clear At;

C = C + 1i.*mmx(actionStr,Ai,Bt,transString);

clear Bt;

C = C - mmx(actionStr,Ai,Bi,transString);
else
At = real(A);
Ai = imag(A);
clear A;
Bt = real(B);
Bi = imag(B);
clear B;

C = mmx(actionStr,At,Bt);

C = C + 1i.*mmx(actionStr,At,Bi);

clear At;

C = C + 1i.*mmx(actionStr,Ai,Bt);

clear Bt;

C = C - mmx(actionStr,Ai,Bi);
end

end

Israel Vaughn

### Israel Vaughn (view profile)

I recently compiled and ran MMX on my moderate workstation, dual Xeon with 28 physical cores (v4 Xeons), 256GB ECC RAM, 2TB Samsung 960 nvme drive. I compiled with the Intel Math Kernel Libraries, and mmx does really poorly with the single or multi-threaded BLAS. It falls off at dimension=40 and is beaten by Matlab for loops.

Graph is here : https://dl.dropboxusercontent.com/u/19842835/mmxSlow.PNG

Any Ideas?

Yuval

### Yuval (view profile)

Hi everyone,

First of all the mention in Loren's blog [1] just wow! I guess the Matlab folks really are listening. Perhaps they will implement mmx's functionality natively one day.

[1] goo.gl/YH6JEM

Second, an apology. I love Matlab and used it constantly for a decade during my academic career, but those days are over. I will not be updating or fixing mmx.

That said, I feel like the code is pretty readable. This is simple C. Sure, LAPACK/BLAS are a pain, and you need to understand threading, but there's no magic here. The package comes with extensive documentation and tests.

Upgrading mmx to support complex numbers, which seems to be the biggest request, should be a nice project for an advanced CS undergraduate.

Convince your professor to give you credit for this, download the code, figure it out, and improve it in any way you like. I'll be happy to answer any questions to help you along.

Ronis Maximidis

### Ronis Maximidis (view profile)

Hi everybody,

Do we no for when to expect the complex matrix support?

It will be great help

yu xx

### yu xx (view profile)

HELLO, THE MMX SEEMS NOT CORRECTLY WHEN THE MATRIX CONTENS COMPLEX NUMBERS

Mark Leznik

### Mark Leznik (view profile)

Thank you for your work Yucal!

Here are two thing we have found when the first option does not work.

The check_dir function should look like this:
function check_dir(dir,files)
if ~isempty(dir)
here = cd(dir);
if nargin == 2
for i = 1:size(files)
try
if isempty(ls(['*' files{i} '.*']))
error('could not find file %s', files{i});
end
catch
cd(here)
error('could not find file %s', files{i});
end
end
end
cd(here);
end

Since, when an exception is thrown the cd(here) is not executed, and in the second iteration of the loop, mmx.cpp cannot be found.

Also, on Ubuntu with Matlab R2015b, we have changed the path from MATLAB_ROOT/extern/lib/glnxa64 to MATLAB_ROOT/bin/glnxa64.

Hope this might help anyone.

Tyler

### Tyler (view profile)

Hi Yuval, we have been using mmx for several years now and it has been very usefull. Thank you for putting in the time to develop it.

I have noticed strange behavior when using the positive definite option, 'P', when solving positive definite systems. I'm not sure if this is a bug or a configuration/use error though. Executing the following statements gives an incorrect solution to the system y = Rx:

%%%%%%%%%%%%%%%%% script 1 %%%%%%%%%%%%%%%%%%

clear all;
R = [2 1;1 2]
y = [3 3].'

%Compute using MATLAB.
'x_matlab = R\y'

%Now use mmx and pos.def. option
x_mmx_P = mmx('backslash',R,y,'P')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

My output is x_matlab = [1.0000 1.0000].' and x_mmx_P = [3 3].';

It appears that somehow x_mmx_P is somehow being set equal to y. However it is possible to obtain the correct solution if mmx is first called without the 'P' option:

%%%%%%%%%%%%%%%%% script 2 %%%%%%%%%%%%%%%%%%

clear all;
R = [2 1;1 2]
y = [3 3].'

%Fist use mmx without the 'P' option. Result
%is correct.

x_mmx_no_P = mmx('backslash',R,y)

%Now use mmx with the 'P' option. Now the result is also correct!

x_mmx_P2 = mmx('backslash',R,y,'P')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

My output for script 2 is x_mmx_no_P = [1 1].' and x_mmx_P2 = [1.0000 1.0000].'

It seems that mmx with the -P option requires some sort of initialization.

These results were obtained on R2015b with Windows 7 and built with the mmx_mkl_multi option. The same results were obtained with R2013a under Mac and linux.

Any idea what is going on here?

Thanks,

Tyler

Tyler

### Tyler (view profile)

I forgot the line 'x_matlab = R\y' in script 1. It should be

%%%%%%%%%%%%%%%%% script 1 %%%%%%%%%%%%%%%%%%

clear all;
R = [2 1;1 2]
y = [3 3].'

%Compute using MATLAB.
'x_matlab = R\y'

%Now use mmx and pos.def. option
x_mmx_P = mmx('backslash',R,y,'P')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Sorry for the confusion.

Guillaume

### Guillaume (view profile)

Hi, is it possible to modify the input to accept one matrix of logical ?

BR

guillaume

Kathleen G Kennedy

Yuval

### Yuval (view profile)

Emily, yes you can, but you'd need to permute the dimensions appropriately. In your example you'd want to A=permute(A,[2 3 1]) so that size(A)=[4 9 2] and then [4 9 2]x[9 7]=>[4 7 2] which you can permute back to [2 4 7].

Emily

### Emily (view profile)

Can I use this to multiply matrices with different dimensionality? For example if I want to multiply a 2x4x9 matrix with a 9x7 matrix to get a 2x4x7 matrix, with a sum occurring over the dimension with length 9?

Philip

### Philip (view profile)

This file is great. One request: can you enable singleton expansion. For example, if a(1,1,:)=[1 2] and b=[1 2; 3 4], I would like to compute c = mmx('mult',a,b), where c(:,:,1) = [1 2; 3 4] and c(:,:,2)= [2 4; 6 8]. Thanks for the submission.

Octavian

### Octavian (view profile)

Dear Yuval,

I want to use mmx for pagefun type multiplications on the cpu. I saw the speed benchmarking for the single and multithreaded BLAS, and I wanted to make sure that the hierarchy (single faster than multithread for large matrices) holds for page by page ND array multiplication. The reason is that both the Intel (for single thread) and Matlab (for multithread) libraries cost ~\$500 a piece, so chosing one not both is of interest. Please advise, thank you,

Octavian

Robbert

### Robbert (view profile)

Dear Yuval,
Thanks for the great tool, I would also very much appreciate complex number support.

Evan

### Evan (view profile)

I was able to compile on OSX 10.8 on 2013b using XCode 5.0.2 by changing link_dir from '/extern/lib/maci64' to 'bin/maci64.' The comments say that the 'mmx_mkl_muli' build uses the BLAS/LAPACK libraries that come with Matlab, but it seems as though it's still trying to use intel mkl libraries. I also had to stop the program from cd-ing out of the mmx package directory so that it wouldn't fail to find mmx.ccp. Works now though!

Evan

### Evan (view profile)

Hi Yuval,

Sorry if this is a naive question, but I'm trying to follow the instructions to download the Intel MKL, and it appears to be quite an expensive software package. Is there another configuration I can use? Thanks!

Yuval

### Yuval (view profile)

Thanks Syed. Unfortunately I don't have time to add that feature. Just use the identity

(B/A)' = (A'\B')

Remember that transposition is a very cheap operation. You can avoid loops by using
my_T = @(x) permute(x,[2 1 3 4])

Syed

### Syed (view profile)

Dear Yuval,

Thank you (and James!) for this wonderfully useful toolbox. I have a question, however. I am trying to run a "forward slash" so I can calculate the product
B \ A / B = (B \ A) / B = B \ (A / B) = B^{-1} A B^{-1}

Your code does great with B\A and I don't believe it's too hard to do A / B = (B' \ A')'. However, because backslash does not have any transpose flags, I can't force your code to do any tweaked backslash operation.

Is it possible for you to add a "forwardslash" option to mmx.cpp? That would be a small but tremendous addition to the Swiss army matrix-knife you've made. Thanks!

Syed

Johannes

Rafael Borges

### Rafael Borges (view profile)

Dear Yuval,

I'm running the Linux 64 bits edition of Matlab. I've followed the instructions you provided inside build_mmx.m for installing Intel MKL and creating libsingle_mkl_ilp64.so, which was then copied to the Matlab extern/lib/glnxa64 folder as instructed. But when I ran build_mmx afterwards, I got the error message

"Compilation of 'mmx_mkl_single' failed with error:
ls: cannot access *small_mkl_ilp64.*: No such file or directory"

and building with single-threaded MKL library failed. What has gone wrong?

Jaap

### Jaap (view profile)

works really great! much faster then using for loops. Indeed, support for complex numbers would be great.

Jaap

Yuval

### Yuval (view profile)

Tamar, I realized what the problem is. LCC is an ANSI C compiler. MMX uses some features of C++ and therefore cannot be compiled with LCC. Both Microsoft and Intel offer free versions of their compilers which matlab recognizes (mex -setup). Try using these.

Tamar Friedlander

### Tamar Friedlander (view profile)

Hi Yuval,

Thanks. It did solve the error message about line 151 (similarly I resolved the error message about lines 238-9), but then I still have the next error messages:

Trying to compile 'mmx_naive', using
-DWIN_SYSTEM,
Error mmx.cpp: 232 illegal statement termination
Error mmx.cpp: 232 skipping `char'
Error mmx.cpp: 232 undeclared identifier `commandStr'
Error mmx.cpp: 232 type error: pointer expected
Error mmx.cpp: 232 operands of = have illegal types `int' and `pointer to char'
Error mmx.cpp: 236 type error: pointer expected
Error mmx.cpp: 255 type error in argument 1 to `mxFree'; found `int' expected `pointer to void'
Error mmx.cpp: 280 illegal statement termination
Error mmx.cpp: 280 skipping `char'
Error mmx.cpp: 280 undeclared identifier `modifierStr'
Error mmx.cpp: 280 type error: pointer expected
Error mmx.cpp: 280 operands of = have illegal types `int' and `pointer to char'
Error mmx.cpp: 281 type error: pointer expected
Error mmx.cpp: 290 type error: pointer expected
Error mmx.cpp: 303 type error in argument 1 to `mxFree'; found `int' expected `pointer to void'
Warning mmx.cpp: 280 possible usage of modifierStr before definition
Error mmx.cpp: 517 illegal statement termination
Error mmx.cpp: 517 skipping `int'
Error mmx.cpp: 517 undeclared identifier `blksz'
Error mmx.cpp: 518 illegal statement termination
Error mmx.cpp: 518 skipping `int'
Error mmx.cpp: 518 too many errors

C:\PROGRA~1\MATLAB\R2011B\BIN\MEX.PL: Error: Compile of 'mmx.cpp' failed.

Compilation of 'mmx_naive' failed with error:
Unable to complete successfully.

What do you suggest to do?

Thanks,
--Tamar.

Yuval

### Yuval (view profile)

Tamar,

I believe you found a bug. Thank you. The declaration of the variable temp in line 151 is in the middle of the function code. This is standard in C++ but not in C, where declarations are supposed to come before statements. I suppose my compiler (VS) was a little too lenient.

Move line 151 to 129 and you should be fine.

Tamar Friedlander

### Tamar Friedlander (view profile)

Hi Yuval,

I tried to compile the mmx on Matlab R2011b (7.13) 32-bit and failed to do so.
My compiler is the default lcc compiler of Matlab (I chose it on mex -setup).

build_mmx failed to compile it on my computer.
I tested my compiler with another mex function that I wrote and it worked well.

For the naive compilation I received the following error messages:

Trying to compile 'mmx_naive', using
-DWIN_SYSTEM,
Error mmx.cpp: .\matrix_fun.c: 151 illegal statement termination
Error mmx.cpp: .\matrix_fun.c: 151 skipping `double'
Error mmx.cpp: .\matrix_fun.c: 151 undeclared identifier `temp'
Warning mmx.cpp: .\matrix_fun.c: 151 Statement has no effect
Error mmx.cpp: .\matrix_fun.c: 238 illegal statement termination
Error mmx.cpp: .\matrix_fun.c: 238 skipping `double'
Error mmx.cpp: .\matrix_fun.c: 238 undeclared identifier `temp'
Error mmx.cpp: .\matrix_fun.c: 239 illegal statement termination
Error mmx.cpp: .\matrix_fun.c: 239 skipping `int'
Error mmx.cpp: .\matrix_fun.c: 239 undeclared identifier `rank'
Error mmx.cpp: 232 illegal statement termination
Error mmx.cpp: 232 skipping `char'
Error mmx.cpp: 232 undeclared identifier `commandStr'
Error mmx.cpp: 232 type error: pointer expected
Error mmx.cpp: 232 operands of = have illegal types `int' and `pointer to char'
Error mmx.cpp: 236 type error: pointer expected
Error mmx.cpp: 255 type error in argument 1 to `mxFree'; found `int' expected `pointer to void'
Error mmx.cpp: 280 illegal statement termination
Error mmx.cpp: 280 skipping `char'
Error mmx.cpp: 280 undeclared identifier `modifierStr'
Error mmx.cpp: 280 type error: pointer expected
Error mmx.cpp: 280 too many errors

Any suggestions what to do?

Many thanks,
--Tamar.

Alexandros Iliopoulos

### Alexandros Iliopoulos (view profile)

Dear Yuval,

Indeed, I'm happy to report that the numerical errors were not due to your package, after all.

Interestingly, the errors I was having were not due to MMX giving me the pseudo-inverse solution (which is what I wanted), but rather they were due to MATLAB's pinv() function, whose output I passed to MMX. In particular, I found out that MATLAB does not use symmetric-matrix-specific algorithms for functions such as pinv, svd, or sqrt, which may result in output matrices having non-zero imaginary terms. Of course, this did mess up the MMX output (as noted earlier), but the input was wrong to begin with.

After going around that problem, MMX worked like a charm!

Alexandros Iliopoulos

Yuval

### Yuval (view profile)

Alexandros,

Thank you so much for your feedback. This is not a bug. As described in the documentation, for rank deficient matrices mmx gives the least-norm (pseudoinverse) solution rather than the maximally-sparse solution. You didn't specify how your matrices were generated but clearly they are not 'in general position'. Try using random A's and you will see no discrepancy (to machine precision). Alternatively, compare mmx's output to pinv(A)*B.

Alexandros Iliopoulos

### Alexandros Iliopoulos (view profile)

Dear Yuval,

It appears to me that your package suffers from numerical artifacts, at least with respect to 'BACKSLASH' computations.

To check, I ran something like this:
% A:3x3x1850, B:3x1x1850
C = mmx( 'backslash', A, B );
D = zeros( size( C ) );
for k = 1 : size(D,3)
D(:,:,k) = A(:,:,k) \ B(:,:,k);
end
E = sum( (C - D).^2, 1 );
stem( permute( E, [1 3 2] ) )

It is readily visible that E is zero everywhere except for 4-5 indices, where the error is pretty large.

Unfortunately, I haven't been able to consruct a MWE with rand() data or something similar. Please PM me if you want me to send you the data that gave me the errors.

I hope that you have the time to work on this issue (assuming that I haven't done something silly), as I think that this is a very handy package to have! I'd be willing to help if I can.

Alexandros Iliopoulos

Yuval

### Yuval (view profile)

Adding complex number support should be pretty easy.
If I see that mmx gets more traction in the community and that more people want it i'd be happy to do it myself.

As it stands it seems that it's buried in a dark corner of file exchange...

The package works great. It all compiled nicely and is a lot faster than everything else out there. Also, the backslash option is really usefull.

Watch out when passing complex numbers though, the function only works on the real part of the data. Would it be possible to extend the function to complex matrices too?

Yuval

### Yuval (view profile)

Jan, look in the comments inside build_mmx.m

Jan Simon

### Jan Simon (view profile)

Do you have explicite suggestions for "linking to a single-threaded BLAS library"?

Yuval