The solution of vector-valued Cahn-Hilliard systems is of interest in many applications. We discuss strategies for the handling of smooth and nonsmooth potentials as well as for different types of constant mobilities. Whereas the use of smooth potentials leads to a system of parabolic partial differential equalities, the nonsmooth ones result in variational inequalities. Concerning the latter, we propose a Moreau-Yosida regularization technique that incorporates the necessary bound constraints. As a result, the variational inequalities are replaced by nonsmooth equations. Due to the use of fully implicit time discretizations, which are the most accurate, we have to solve in every time step nonlinear smooth or nonsmooth equations. This is done by standard Newton methods in the smooth case, and by semismooth Newton methods in the nonsmooth case. At the heart of both methods lies the solution of large and sparse linear systems for which we propose the use of preconditioned Krylov subspace solvers using an effective Schur complement approximation. Numerical results illustrate the efficiency of our approach. In particular, we numerically show mesh and phase independence of the developed preconditioner in the smooth case. The results in the nonsmooth case are also satisfying and the preconditioned version always outperforms the unpreconditioned one.