Abstract <p>Some iterative processes in Krylov subspaces are considered for solving systems of linear algebraic equations (SLAE) with high-order sparse matrices that arise in grid approximations of multidimensional boundary value problems. The SLAE are preconditioned by a uniform combined method that includes domain decomposition and recursive application of a two-grid algorithm, which are implemented by forming block-tridiagonal algebraic and grid structures inverted by using incomplete factorization and diagonal compensation. Stability and convergence of iterations are studied for some Stieltjes systems. Parallelization and generalization of the methods to wider classes of relevant practical problems are discussed.</p>