Code covered by the BSD License  

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

5.0

5.0 | 8 ratings Rate this file 35 Downloads (last 30 days) File Size: 460 KB File ID: #37515
image thumbnail

MMX - Multithreaded matrix operations on N-D matrices

by

 

16 Jul 2012 (Updated )

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

| Watch this File

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
clear the threads from memory.

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
about compilation issues and options

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
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (24)
28 Jul 2014 Robbert

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

10 Mar 2014 Evan

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

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!

08 Oct 2013 Yuval

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])

08 Oct 2013 Syed

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!

08 Oct 2013 Syed  
02 Oct 2013 Johannes  
06 Aug 2013 Rafael Borges

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?

05 Aug 2013 Jaap

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

05 Aug 2013 Jaap  
08 Jul 2013 Yuval

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.

07 Jul 2013 Tamar Friedlander

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.

06 Jul 2013 Yuval

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.

05 Jul 2013 Tamar Friedlander

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.

04 Jun 2013 Alexandros Iliopoulos

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!

04 Jun 2013 Alexandros Iliopoulos  
04 Jun 2013 Yuval

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.

03 Jun 2013 Alexandros Iliopoulos

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.

03 Jun 2013 Alexandros Iliopoulos  
12 Apr 2013 Yuval

Thanks Adam!
Please spread the word.
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...

12 Apr 2013 Adam Cooman

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

Jan, look in the comments inside build_mmx.m

11 Sep 2012 Jan Simon

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

16 Jul 2012 Yuval

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

Updates
08 Jul 2013

Some more variable declaration fixes for lcc compiler (not tested).

Contact us