function p5new
  a = random_pentadiag(6)
  [l,u] = crout_pentadiag(a)
  a_minus_lu = a - l*u
end
function [l, u] = crout_pentadiag(a)
  [n, dummy] = size(a);

  l = zeros(n,n);
  u = zeros(n,n);

  l(1:3,1) = a(1:3,1);
  u(1,1) = 1;

  % 1 + 3*2 = 7 op
  u(1,2) = a(1,2)/l(1,1);
  u(2,2) = 1;
  l(2,2) = a(2,2) - l(2,1)*u(1,2);
  l(3,2) = a(3,2) - l(3,1)*u(1,2);
  l(4,2) = a(4,2) - l(4,1)*u(1,2);

  for j = 3:n-2 % 10 op
    u(j-2,j) = a(j-2,j)  / l(j-2,j-2); % 1 op
    u(j-1,j) = (a(j-1,j) - l(j-1,j-2)*u(j-2,j))/l(j-1,j-1); % 3 op
    u(j,j) = 1;
    l(j,j)   =  a(j,j)   - l(j,j-2)*u(j-2,j) - l(j,j-1)*u(j-1,j); % 4 op
    l(j+1,j) =  a(j+1,j) - l(j+1,j-1)*u(j-1,j); % 2 op
    l(j+2,j) =  a(j+2,j);
  end

  j = n-1; % 10 op
  u(j-2,j) = a(j-2,j)  / l(j-2,j-2);
  u(j-1,j) = (a(j-1,j) - l(j-1,j-2)*u(j-2,j))/l(j-1,j-1);
  u(j,j) = 1;
  l(j,j)   =  a(j,j)   - l(j,j-2)*u(j-2,j) - l(j,j-1)*u(j-1,j);
  l(j+1,j) =  a(j+1,j) - l(j+1,j-1)*u(j-1,j);

  j = n; % 8 op
  u(j-2,j) = a(j-2,j)  / l(j-2,j-2);
  u(j-1,j) = (a(j-1,j) - l(j-1,j-2)*u(j-2,j))/l(j-1,j-1);
  u(j,j) = 1;
  l(j,j)   =  a(j,j)   - l(j,j-2)*u(j-2,j) - l(j,j-1)*u(j-1,j);

  lu = l*u

  assert(abs(a-l*u) < 1e-15)
end
function a = random_pentadiag(n)
  a = zeros(n,n);

  for i = 1:n
    for j = max(i-2,1):min(n,i+2)
      a(i,j) = rand;
    end
  end
end
function assert(v)
  if (v == 0)
    error('Assertion failed.')
  end
end
