We present a numerical study of the Kelvin-Helmholtz instability occuring at the interface of a shear flow configuration in compressible magnetohydrodynamics. The calculations cover two and three dimensional configurations. The results are obtained using the Versatile Advection Code (VAC, see http://www.phys.uu.nl/~toth), a software package designed to solve general systems of conservation laws. This includes the shallow water model, the Euler equations and the magnetohydrodynamic equations (where magnetic fields control the plasma flow through the additional Lorentz force). The 2D results consider an initial magnetic field aligned with the shear flow, and analyze the differences between cases where the field is unidirectional everywhere and where the field changes sign at the interface. The 3D simulations consider shear flow in a cylindrical jet configuration, embedded in a uniform magnetic field. The growth of linear perturbations at specified poloidal and axial mode numbers demonstrate intricate non-linear coupling effects. We identify the physical mechanims leading to induced secondary Kelvin-Helmholtz instabilities at higher mode numbers. The initially weak magnetic field becomes locally dominant in the non-linear dynamics before and during saturation.