We consider different parallel versions of the least squares methods in the Krylov subspaces which are based on computing various basis vectors. These algorithms are used for solving very large real, non-symmetric, in gerenal, sparse systems of linear algebraic equations (SLAEs) which arise in grid approximations of multi-dimensional boundary value problems. In particular, the Chebyshev acceleration approach, steepest descent and minimal residual, conjugate gradient and conjugate residual are applied as preliminary iterative processes. The resulting minimization of residuals is provided by the block, or implicit, orthogonalization procedures. The properties of the Krylov approaches proposed are analysed in the “pure form”, i.e. without preconditioning. The main criteria of parallelezation are estimated. The convergence rate and stability of the algorithms are demonstated on the results of numerical experiments for the model SLAEs which present the exponential fitting approximation of diffusion-convection equations on the meshes with various steps and with different coefficients.