Parallel methods for solving saddle-type algebraic systems that are relevant for modeling processes and phenomena in the problems of electromagnetism, hydro-gas dynamics, elastoplasticity, filtration and other applications are considered. Preconditioned iterative processes in the Krylov subspaces, including the efficient generalization of the Golub-Kahan-Arioli bidiagonalization method, are investigated as applied to large SLAEs with sparse matrices that arise when approximating multi-dimensional boundary value problems with a complex geometric configuration of computational domains and the contrasting material properties of various media on unstructured grids. It is supposed to store the matrices in compressed formats that require special technologies for working with big data. The parallelization of the proposed class of block algorithms is carried out by means of hybrid programming on supercomputers of a heterogeneous architecture with distributed and hierarchical shared memory, using the means of inter-node message transmission, multi-threaded computing, operation vectorization. A comparative analysis of various algorithmic approaches is carried out on the basis of the estimates of the performance and resource intensity of the corresponding software implementations.