function C = recgemm(A, B)
  
  [m, k] = size (A);
  [k2, n] = size (B);
  assert (k == k2);

  old_recdepth = max_recursion_depth (max (2*[m, k, n]));

  C = zeros (m, n);
  C = innergemm (m, n, k, A, B, C);

  max_recursion_depth (old_recdepth);

endfunction

function C = innergemm(m, n, k, A, B, C)
  DIMTHRESH = 16;

  if m < DIMTHRESH && n < DIMTHRESH && k < DIMTHRESH,
    C += A*B;
    return;
  endif

  lrgdim = max([m, n, k]);
  if lrgdim == m,
    mlow = floor(m/2);
    mhigh = m-mlow;
    C(1:mlow,:) = innergemm (mlow, n, k, A(1:mlow,:), B, C(1:mlow,:));
    C((mlow+1):m,:) = innergemm (mhigh, n, k, A((mlow+1):m,:), B, C((mlow+1):m,:));
  elseif lrgdim == n,
    nlow = floor(n/2);
    nhigh = n-nlow;
    C(:,1:nlow) = innergemm (m, nlow, k, A, B(:,1:nlow), C(:,1:nlow));
    C(:,(nlow+1):n) = innergemm (m, nhigh, k, A, B(:,(nlow+1):n), C(:,(nlow+1):n));
  else
    klow = floor(k/2);
    khigh = k-klow;
    C = innergemm (m, n, klow, A(:,1:klow), B(1:klow,:), C);
    C = innergemm (m, n, khigh, A(:,(klow+1):k), B((klow+1):k,:), C);
  endif

endfunction

%!test
%!  m = 17;
%!  n = 19;
%!  k = 22;
%!  A = rand(m, k);
%!  B = rand(k, n);
%!  C = A*B;
%!  Cr = recgemm(A, B);
%!  assert(Cr, C, -2*max([m,n,k])*eps)

%!test
%!  m = 746;
%!  n = 221;
%!  k = 957;
%!  A = rand(m, k);
%!  B = rand(k, n);
%!  C = A*B;
%!  Cr = recgemm(A, B);
%!  assert(Cr, C, -2*max([m,n,k])*eps)