Least Squares Collocation (LSC) is used for the modeling of the gravity field, including prediction and error estimation of various quantities. The method requires that as many unknowns as number of data and parameters are solved for. Cholesky reduction must be used in a non-standard form due to missing positive-definiteness of the equation system. Furthermore the error estimation produces a rectangular or triangular matrix which must be Cholesky reduced in the non-standard manner. LSC have the possibility to add new sets of data without processing previously reduced parts of the equation system. Due to these factors standard Cholesky reduction programs using multi-processing cannot easily be applied. We has therefore implemented Fortran Open Multi-Processing (OpenMP) to the non-standard Cholesky reduction. In the computation of matrix elements (covariances) as well as the evaluation spherical harmonic series used in the remove/restore setting we also take advantage of multi-processing. We describe the implementation using quadratic blocks, which aids in reducing the data transport overhead. Timing results for different block sizes and number of equations are presented. OpenMP scales favorably so that e.g. the prediction and error estimation of grids from GOCE TRF vertical gradient data and ground gravity data can be done in the less than two hours for a 25° by 25° area with data selected close to 0.125° nodes.