function [L,U] = BlockLU(A,r)
% function [L,U] = BlockLU(A,r)
% Block LU without pivoting.
% A is nxn and has an LU factorization.
% L is nxn and unit lower triangular.
% U is nxn and upper triangular.
% r is the block size and n/r = N is an integer.
% A = LU
% GVL4: Algorithm 3.2.3
[n,n] = size(A); L = zeros(n,n); U = zeros(n,n);
N = n/r;
for k=1:N
i = (k-1)*r;
cols = i+1:i+r;
if k<N
[L(i+1:n,cols),U(cols,cols)] = RectLU(A(i+1:n,cols));
alfa = i+r+1:n;
U(cols,alfa) = L(cols,cols)\A(cols,alfa);
A(alfa,alfa) = A(alfa,alfa)-L(alfa,cols)*U(cols,alfa);
else
[L(cols,cols),U(cols,cols)] = RectLU(A(cols,cols));
end
end

function [L,U] = RectLU(A)
% function [L,U] = RectLU(A)
% Outer product LU for tall rectangular matrices without pivoting
% A is nxr and has an LU factorization, n>=r.
% L is nxr and unit lower triangular.
% U is rxr and upper triangular.
% A = LU
% GVL4: Section 3.2.10
[n,r] = size(A);
for k=1:r
rho = k+1:n;
% Compute L(k+1:n,k)...
A(rho,k) = A(rho,k)/A(k,k);
if k<r
% The outer product update...
mu = k+1:r;
A(rho,mu) = A(rho,mu) - A(rho,k)*A(k,mu);
end
end
L = tril(A);
for k=1:r
L(k,k) = 1;
end
U = triu(A(1:r,1:r));