The paper studies an approach for solving a system of linear equations \(Ax=b\) numerically which consists of two steps. Firstly, a preconditioner which consists of some steps of the Gaussian elimination is used which produces zero entries in the matrix. Secondly, a classical iterative scheme is applied to the preconditioned system with the matrix \(\tilde A\) or to a reduced system with the system matrix \(\tilde A_1\) where the columns of \(\tilde A\) with the zero entries and the corresponding rows are eliminated. This algorithm has been previously studied for one Gaussian elimination step and pointwise iterative methods. In this paper, it is extended to more Gaussian elimination steps and blockwise iterative methods. The convergence of this algorithm is proven for nonsingular M-matrices \(A\), the block-Jacobi method, some block-Gauss-Seidel-type methods and block-successive overrelaxation-type methods applied to the systems with matrix \(\tilde A\) or \(\tilde A_1\), respectively, by bounding the spectral radii of the iteration matrices from above strictly by 1. The results are extended to nonsingular \(p\)-cyclic consistently ordered matrices. Their extensions to some types of singular matrices is discussed. For supporting the theoretical results, spectral radii of the iteration matrices for one particular M-matrix are presented. Finally, a numerical example for a finite difference discretization of the Poisson equation is presented, where two of the studied methods are compared to preconditioned conjugated conjugate gradient (PCG) methods with an incomplete Cholesky preconditioner. If the mesh width becomes finer, PCG becomes more and more superior.