Code covered by the BSD License

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

5.0
5.0 | 9 ratings Rate this file 75 Downloads (last 30 days) File Size: 460 KB File ID: #37515 Version: 1.2

# MMX - Multithreaded matrix operations on N-D matrices

by

### Yuval (view profile)

16 Jul 2012 (Updated )

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

File Information
Description

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

Acknowledgements

Mtimesx Fast Matrix Multiply With Multi Dimensional Support inspired this file.

Required Products MATLAB
MATLAB release MATLAB 7.14 (R2012a)
Other requirements Works better if you can link it to a single-threaded BLAS library
11 Sep 2015 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.

20 Apr 2015 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

Comment only
28 Jul 2014 Robbert

### Robbert (view profile)

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

10 Mar 2014 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!

09 Mar 2014 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!

Comment only
08 Oct 2013 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])

Comment only
08 Oct 2013 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!

Comment only
08 Oct 2013 Syed

### Syed (view profile)

02 Oct 2013 Johannes

### Johannes (view profile)

06 Aug 2013 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?

Comment only
05 Aug 2013 Jaap

### Jaap (view profile)

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

Comment only
05 Aug 2013 Jaap

### Jaap (view profile)

08 Jul 2013 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.

Comment only
07 Jul 2013 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.

Comment only
06 Jul 2013 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.

Comment only
05 Jul 2013 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.

Comment only
04 Jun 2013 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!

Comment only
04 Jun 2013 Alexandros Iliopoulos

### Alexandros Iliopoulos (view profile)

04 Jun 2013 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.

Comment only
03 Jun 2013 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.

Comment only
03 Jun 2013 Alexandros Iliopoulos

### Alexandros Iliopoulos (view profile)

12 Apr 2013 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...

Comment only

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?

19 Sep 2012 Yuval

### Yuval (view profile)

Jan, look in the comments inside build_mmx.m

Comment only
11 Sep 2012 Jan Simon

### Jan Simon (view profile)

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

Comment only
16 Jul 2012 Yuval

### Yuval (view profile)

This package's homepage is here
http://www.cs.washington.edu/people/postdocs/tassa/code/#mmx

Comment only