clear;
k = 300000;
d = 3;
n = 8;

B = rand(k*d, n);
A = rand(k*d, d);

fprintf('Loop solution ... '); tic;
x1 = zeros(k*d, n);
for i = 1 : k
    idx = i*d-(d-1):i*d;
    x1(idx, :) = A(idx, :) \ B(idx, :);
end
fprintf('Ready in %.4f seconds.\n', toc);

fprintf('Block solution ... '); tic;
r = reshape(repmat((1:k*d).', 1, 3), [], 1);
c = bsxfun(@plus, kron(d*(0:k-1).', [1;1;1]), 1:d);
As = sparse(r, c, A(:));
x2 = As \ B;
fprintf('Ready in %.4f seconds.\n', toc);

x2 = multi_block_inverse(A, B);

err = norm(x1-x2)/norm(x1);
fprintf('Equal: %d\n', err < 10*eps);