function [L,U] = blu1(A,r)
% Block LU factorization s.t. A = L*U
%   A: square matrix
%   r: block size

n = length(A);
k = 1;
while k <= n
    m = min(n,k+r-1);
    [L11,U11] = mylu1(A(k:m,k:m)); % A11 = L11*U11
    A(k:m,k:m) = tril(L11,-1) + U11; % overwrite A11
    U12 = L11\A(k:m,m+1:n); % solve L11*X = A12 for X
    L21 = A(m+1:n,k:m)/U11; % solve Y*U11 = A21 for Y
    A(k:m,m+1:n) = U12; % overwrite A12
    A(m+1:n,k:m) = L21; % overwrite A21
    A(m+1:n,m+1:n) = A(m+1:n,m+1:n) - L21*U12; % update A22
    k = k + r;
end
L = eye(n) + tril(A,-1);
U = triu(A);
