function Feb8LecSol()
clc
n = 1024;
A = rand(n,n);
for m = [4 8 16 32 64]
   nMin = m;
   MMFlopsIn = 0;
   [L,U,MMFlopsOut] = luRecur(A,nMin,m,MMFlopsIn);
   AllFlops = 2*n^3/3;
   ratio = MMFlopsOut/AllFlops;
   fprintf('n = %4d     block size = %2d    (#MM Flops) / (Total Flops) = %5.3f\n',n,m,ratio)
end

  function [L,U,MMFlopsOut] = luRecur(A,nMin,m,MMFlopsIn)
% A = LU
% Recursive LU with block size m and base-case threshold nMin
% MMFlopsOut  = MMFlopsIn  + number of MM flops at this call
[n,n] = size(A);
if n<=nMin
    [L,U] = lu(A);
    MMFlopsOut = MMFlopsIn;
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,MMFlopsOut] = luRecur(A22-L21*U12,nMin,m,MMFlopsIn);
    MMFlopsOut = MMFlopsOut + 2*(n-m)*(n-m)*m;
    L = [L11 zeros(m,n-m) ; L21 L22];
    U = [U11 U12 ; zeros(n-m,m) U22];
end