We present the results of an extensive lattice calculation of the renormalization constants of bilinear and four-quark operators for the non-perturbatively O(a)-improved Wilson action. The results are obtained in the quenched approximation at four values of the lattice coupling by using the non-perturbative RI/MOM renormalization method. Several sources of systematic uncertainties, including discretization errors and final volume effects, are examined. The contribution of the Goldstone pole, which in some cases may affect the extrapolation of the renormalization constants to the chiral limit, is non-perturbatively subtracted. The scale independent renormalization constants of bilinear quark operators have been also computed by using the lattice chiral Ward identities approach and compared with those obtained with the RI-MOM method. For those renormalization constants the non-perturbative estimates of which have been already presented in the literature we find an agreement which is typically at the level of 1%.