This paper deals with the solution of linear algebraic equations encountered in geostatistics. The kriging estimator in a global neighbourhood requires the inversion of a large symmetric positive definite matrix K or the search of the solution by any other direct or iterative method. Classical algorithms such as Gauss elimination, Jacobi methods or UtDU decomposition have been extensively used in the past (Journel and Huijbregts, 1981; Davis and Grivet, 1984). However, some of these methods require the knowledge of the second members and do not take into account the accuracy of the solution. An alternative method combining a direct and an iterative procedure is presented here. The basic idea is a UtUDU decomposition of the kriging matrix, then an iterative improvement of the original solution until a given accuracy is reached. The effective computer time depends on n^3 (more precisely n^3 /6 + n^2 (i+1) where n being the number of equations, i the number of iterations usually 1 to 2) which is similar to direct method (n^3 /6 + n^2), however the accuracy of the solution is guaranteed. This technique could be applied to solve direct or dual kriging systems. A similar technique based on the LU decomposition could be extended to the non stationary case.