function Feb8Lec()
n = 512;
m = 32;
nMin = 16;
A = randn(n,n);
[L,U] = luRecur(A,nMin,m);
err = norm(A - L*U,1)

  function [L,U] = luRecur(A,nMin,m)
% A = LU
% Recursive LU with block size m and base-case threshold nMin
[n,n] = size(A);
if n<=nMin
    [L,U] = lu(A);
else
%    A11   A12        L11   0    U11  U12
%                 =            *
%    A21   A22        L21  L22   0    U22
    A11 = A(1:m,1:m);
    A21 = A(m+1:n,1:m);
    A12 = A(1:m,m+1:n);
    A22 = A(m+1:n,m+1:n);
% A11 = L11*U11 A2 and so
    [L11,U11] = lu(A11);
% A12 = L11*U12 and so
    U12 = L11\A12;
% A21 = L21*U11 and so
    L21 = A21/U11;
% A22 = L21*U12 + L22*U22 and so
    [L22,U22] = lu(A22-L21*U12);
    L = [L11 zeros(m,n-m) ; L21 L22];
    U = [U11 U12 ; zeros(n-m,m) U22];
% The A22 flops are level-3 flops. That is, there are matrix-matrix multiplication flops
end