function p7
  disp('-------------------------------')
  disp('Part a)')
  disp('-------------------------------')
  A = [16, -28, 0;-28, 53,10;0,10,29];
  b = [8, -2, 38]';
  solution = cholesky_solve(A, b)
  disp('-------------------------------')
  disp('Part b)')
  disp('-------------------------------')
  A = [9/4,3,3/2;3,25/4,7/2;3/2,7/2,17/4];
  b = [3,1,9]';
  solution = cholesky_solve(A, b)
  disp('-------------------------------')
  disp('Part c)')
  disp('-------------------------------')
  A = [4,-2,-2,0;-2,5,1,-2;-2,1,10,3;0,-2,3,18];
  b = [4,-4,4,-13]';
  solution = cholesky_solve(A, b)
  disp('-------------------------------')
  disp('Part d)')
  disp('-------------------------------')
  A = [1,-2,3,-2;-2,20,-2,8;3,-2,11,-5;-2,8,-5,9];
  b = [15,-12,56,-35]';
  solution = cholesky_solve(A, b)
end

function sol = cholesky_solve(A, b)
  A
  L = mycholesky(A)

  [n, dummy] = size(A);

  % backsubsitution for L (lower)
  tempsol = zeros(n,1);
  for row = 1:n
    tempsol(row) = b(row) - L(row,1:row-1)*tempsol(1:row-1);
    tempsol(row) = tempsol(row) / L(row,row);
  end

  % backsubsitution for L' (upper)
  sol = zeros(n,1);
  for step = 0:n-1
    row = n-step;
    sol(row) = tempsol(row) - L(row+1:n,row)'*sol(row+1:n);
    sol(row) = sol(row) / L(row,row);
  end

  residual = A*sol - b
end

function L = mycholesky(A)
  [n, dummy] = size(A);

  % coded up directly from p. 216

  for k=1:n
    diag_squared = A(k,k);

    for j=1:k-1
      diag_squared = diag_squared - L(k,j)^2;
    end

    L(k,k) = sqrt(diag_squared);

    for i=k+1:n
      entry = A(i,k);

      for j=1:k-1
        entry = entry - L(i,j)*L(k,j);
      end

      L(i,k) = entry/L(k,k);
    end
  end

  if (norm(A-L*L',2) > 1e-13)
    error('Cholesky did not work out.')
  end
end
