Efficient Multiplication of Submatrices whose results are to be safed in Submatrices

Hey Guys,

trying to find an efficient solution for quite some time now.

Given two matrices A,B, say both of size 1000x1000, compute A(1:100,1:100)*B(1:100,1:100) and safe that in (overwrite) A(1:100,1:100).

What's the most efficient way to do that? Since it is a mantra for Matlab to always preallocate, I precisely did that with worse results than before (before: appending new row and column in each loop). I've read about efficient BLAS-routines in Fortran and C for multiplying submatrices and saving the result in submatrices since there they don't actually "extract" a submatrix but merely point to which parts of the matrices you want to multiply and therefore you can still make use of the efficient routines available. Are there similar routines for Matlab or could you provide me with a mex-file since I can't program in C?

Original way without preallocation (exemplary, simplified code, say matrix A consists of my data):

The first part takes roughly 16 seconds, the second 14 seconds even though I preallocated. Having read some threads, I've heard that this is due to Matlab having to copy the submatrix at every extraction step and therefore having to find new memory.

In my real code, actually the latter variant is way slower than the first one. I can post it, if anyone is interested.

Thanks in advance!

tic
A=1;
for j=1:1000
    A=[A zeros(j,1);zeros(1,j) 1]*rand(j+1,j+1);
end
toc

tic
B=zeros(1000,1000);
for j=1:1000
    B(j+1,j+1)=1;
    B(1:j+1,1:j+1)=B(1:j+1,1:j+1)*rand(j+1,j+1);
end
toc

[James Tursa]

Yes, the BLAS library that MATLAB uses is already set up to work with sub-matrices without copying them first. That part can easily be done in a mex routine. The difficult part is stuffing the result back into one of the inputs. Because of the data sharing model that MATLAB uses, you might end up inadvertently overwriting other variables in the background. There are no official API functions available that can tell you when this data sharing is taking place, either. So you would be left with hacks for this, or you take special care to make sure the variables are not shared prior to entering the mex routine.

[Benjamin Kiefer]

Thanks for your answer!

Two reply questions:

1.

"That part can easily be done in a mex routine"

But it can be done with traditional Matlab as well, since it is already implemented there? Or what it is the advantage of doing it via a mex file?

2.

How can I take special care s.t. the variables are not shared using purely Matlab?

Best

Answer 1

[amin ya]

Why you don't simly use :

A(1:100,1:100)*B(1:100,1:100)

Matlab matrix multiplication already uses parallelization which is very fast

for preallocation use:

A=rand([1000,1000],'distributed')
B=rand([1000,1000],'distributed')

[Benjamin Kiefer]

Because in my actual code, extracting submatrices and multiplying them together just to be saved in submatrices again takes way longer than iteratively adding a column and row. Just as described in the little example method 2 takes equally much time compared to method 1 even though it's just a bunch of (sub)matrix multiplications. Is the only way to write in C?

My matrices are not too so big that distributing is worth doing. At least at my tests it performed much worse.