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

  % coded up directly from p. 216

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

    % 2(k-1) operations
    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);

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

      % 1 operation
      L(i,k) = entry/L(k,k);
    end
  end
end

