function R = cholesky_inf(A, tolerance)
  %
  % Peform the Cholesky-Infinity factorization of ``A``. This differs
  % from the ordinary Cholesky factorization in that it handles
  % matrices which are not positive-definite.
  %
  % If a non-positive diagonal entry appears, it is replaced with
  % ``inf`` and the rest of that row is set to zero. This will force a
  % zero in the corresponding row of the solution vector.
  %
  % This function is based on the description of MATLAB's `cholinc(A,
  % 'inf')` function which has been deprecated and removed.
  %
  % INPUT:
  %
  %   - ``A`` -- The matrix to factor.
  %
  %   - ``tolerance`` -- (default: eps) The zero tolerance used to test
  %                      the diagonal entries. If `A(i,i) <= tolerance`,
  %                      then we consider `A(i,i) to be zero.
  %
  % OUTPUT:
  %
  %   - ``R`` -- The upper-triangular Cholesky factor of ``A``. ``R``
  %              will be sparse if ``A`` is.
  %
  if (nargin < 2)
    tolerance = eps;
  end

  [m,n] = size(A);

  if (issparse(A))
     % If `A` is sparse, `R` should be too.
    R = sparse(m,n);
  else
    R = zeros(m,n);
  end

  for i = [ 1 : n ]
    col_above = R(1:i-1, i);
    R(i,i) = A(i,i) - col_above'*col_above;

    if (R(i,i) <= tolerance)
      % A is not positive-definite. Rather than throw an error here,
      % we set this pivot to 'inf' and the rest of the row to zero.
      R(i,i) = inf;
      R(i, i+1 : end) = 0;
    else
      % Normal case with the pivot greater than zero.
      R(i,i) = sqrt(R(i,i));

      for j = [ i+1 : n ]
	row1 = R(1:i-1, i);
	row2 = R(1:i-1, j);
	R(i,j) = (A(i,j) - row1'*row2)/R(i,i);
      end
    end

  end

end
