In this paper, a polynomial-time algorithm is introduced for the solution of a generalized linear complementarity problem with upper bounds (BLCP) of the form w = q + Mz, a ⩽ z ⩽ b, { zi = ai ⇒ wi ⩾ 0, zi = bi ⇒ wi ⩽ 0, ai < zi < bi ⇒ wi = 0 } (i = 1 ,..., n), where q ∈ ℝn, M is a nonsingular M-matrix of order n, and −∞ ⩽ ai < bi ⩽ ∞ for all i = 1 ,..., n. The algorithm generalizes Pang's method for the solution of the same problem when all the bounds are finite. An implementation of the algorithm is discussed for the solution of large-scale BLCPs with a symmetric or unsymmetric nonsingular M-matrix. Computational experience indicates that the algorithm is quite efficient in solving large-scale BLCPs with sparse nonsingular M-matrices.