title: block matrix multiplication with matlab cells
body: <p>I want to do block matrix multiplication of A and B using cells in Matlab. More specifically suppose </p>
  
  <pre><code>a=
  
   1     1     2     2
   1     1     2     2
   3     3     4     4
   3     3     4     4
  
  b=
  
   2     2     4     4
   2     2     4     4
   6     6     8     8
   6     6     8     8
  </code></pre>
  
  <p>We can now convert a and b to cell arrays containing their blocks. </p>
  
  <pre><code>A = mat2cell(a,[2,2],[2,2])
  
  ans = 
  
  [2x2 double]    [2x2 double]
  [2x2 double]    [2x2 double]
  
  B = mat2cell(b,[2,2],[2,2])
  
  ans = 
  
  [2x2 double]    [2x2 double]
  [2x2 double]    [2x2 double]
  </code></pre>
  
  <p>I need a function C=foo(A,B) that will return a cell array C such that the blocks of C are the blocks of matrix product A*B, for example in this case:</p>
  
  <pre><code>C{1,1} = A{1,1}*B{1,1} + A{1,2}*B{2,1}
  C{1,2} = A{1,1}*B{1,2} + A{1,2}*B{2,2}
  ...
  </code></pre>
  
  <p>and <code>cell2mat(C)</code> should return:</p>
  
  <pre><code>ans =
  
  28    28    40    40
  28    28    40    40
  60    60    88    88
  60    60    88    88
  </code></pre>
  
  <ul>
  <li><p>The reason I can't just do <code>cell2mat(A)*cell2mat(B)</code> is because in my application most blocks are zero and that would be inefficient. </p></li>
  <li><p>Even though most blocks are zero, I also can't do <code>sparse(cell2mat(A))*sparse(cell2mat(B))</code> because the blocks that are non-zero are fully dense so that would also be inefficient. </p></li>
  </ul>
  
  <p>Is there a better way to do this problem without using clunky and slow for loops? Thanks!</p>
  
  <p><strong>Edit</strong>: I wrote a small code that illustrates what I want to do. It is slow, however, and I'm wondering if there's a better way.</p>
  
  <pre><code>function C = celltimes(A,B,nn,blocksize)
  
  C = cell(nn);
  [C{:}] = deal(sparse(blocksize,blocksize));
  
  for ii = 1:nn
     for jj = 1:nn
        row = A(ii,:);
        col = B(:,jj);
        for kk = 1:nn
           if ( nnz(row{kk}~=0) &amp;&amp; nnz(col{kk}~=0) )
              C{ii,jj} = C{ii,jj}+row{kk}*col{kk};
           end
        end
     end
  end
  </code></pre>
  
  <p>and a test code: </p>
  
  <pre><code>%test
  
  nn = 3; %number of blocks
  blocksize = 3; %block size
  
  a = randi([0,10],nn*blocksize)
  b = randi([0,10],nn*blocksize)
  
  A = mat2cell(a,repmat(blocksize,[1,nn]),repmat(blocksize,[1,nn]));
  B = mat2cell(b,repmat(blocksize,[1,nn]),repmat(blocksize,[1,nn]));
  
  C = celltimes(A,B,nn,blocksize);
  
  %verify result
  c = a*b;
  max(max(abs(cell2mat(C)-c)))
  </code></pre>