We study iterative methods for solving large systems of linear algebraic equations (SLAEs) with sparse matrices in Krylov subspaces, applied to a three-dimensional Laplace operator problem discretized on cubic unstructured meshes. The preconditioner for the original SLAE is constructed using recursive algorithms and data structures, which generate operators for the multigrid incomplete factorization method. Here, the forward step corresponds to the traditional reduction step, while the backward step handles solution prolongation. We discuss the implementation of these recursive algorithms and data structures within the INMOST and PETSc software frameworks. Additionally, we investigate properties of node renumbering for initial unstructured meshes. Numerical results are presented for methodological applied problems with data typical of geophysical core modeling.