In this paper, we present the block least squares method for solving nonsymmetric linear systems with multiple righthand sides. This method is based on the block bidiagonalization. We first derive two algorithms by using two different
convergence criteria. The first one is based on independently minimizing the 2-norm of each column of the residual matrix
and the second approach is based on minimizing the Frobenius norm of residual matrix. We then give some properties of
these new algorithms. Finally, some numerical experiments on test matrices from Harwell–Boeing collection are presented
to show the efficiency of the new method.