This paper develops and experimentally investigates macrogrid domain decomposition methods for solving large systems of linear algebraic equations (SLAEs) with sparse symmetric matrices obtained from grid approximations of multidimensional boundary value problems. The proposed algorithms are based on constructing two-layer macro-grids and special ordering of nodes according to their belonging to different topological primitives of the macro-grid: macro-nodes, macro-edges, macro-faces, and subdomains. With consistent numbering of vector components, the SLAE matrices in the three-dimensional case take a block-tridiagonal form of fourth order. For its solution, an incomplete factorization algorithm is used, based on block bisection of the original matrix and application of parallel direct or preconditioned iterative algorithms in subdomains. The justification of the proposed methods is given for symmetric positive definite (s.p.d.) matrices.