Every iteration of an interior point method of large scale linear programming requires computing at least one orthogonal projection of the objective function gradient onto the null space of a linear operator defined by the problem constraint matrix A. The orthogonal projection itself is in turn dominated by the inversion of the symmetric matrix of form AθAT, where θ is a diagonal weighting matrix. In this paper several specific issues of implementation of the Cholesky factorization that can be applied for solving such equations are discussed. The code called CHFACT being the result of this work is shown to produce comparably sparse factors as the state-of-the-art implementation of the Cholesky decomposition of George and Liu. It has been used for computing projections in an efficient implementation of a higher order primal-dual interior point method of Altman and Gondzio. Although the primary aim of developing CHFACT was to include it into an LP optimizer, the code may equally well be used to solve general large sparse positive definite systems arising in different applications.