The paper considers preconditioned iterative methods in Krylov subspaces for solving large systems of linear algebraic equations with sparse coefficient matrices arising in solving multidimensional boundary-value problems by finite volume or finite element methods of different orders on unstructured grids. Block versions of the weighted Cimmino methods, based on various orthogonal and/or variational approaches and realizing preconditioning functions for two-level multi-preconditioned semi-conjugate residual algorithms with periodic restarts, are proposed. At the inner iterations between restarts, additional acceleration is achieved by applying deflation methods, providing low-rank approximations of the original matrix and playing the part of an additional preconditioner. At the outer level of the Krylov process, in order to compensate the convergence deceleration caused by restricting the number of the orthogonalized direction vectors, restarted approximations are corrected by using the least squares method. Scalable parallelization of the methods considered, based on domain decomposition, where the commonly used block Jacobi–Schwarz iterative processes is replaced by the block Cimmino–Schwarz algorithm, is discussed. Hybrid programming technologies for implementing different stages of the computational process on heterogeneous multi-processor systems with distributed and hierarchical shared memory are described.