This article is dedicated to the analysis of quadratic optimal control of linear complementarity systems (LCS), which are a class of strongly nonlinear and nonsmooth dynamical systems. Necessary first-order conditions are derived, that take the form of an LCS with inequality constraints, hence are numerically tractable. Then two numerical solvers are proposed, for the direct and the indirect approaches. They take advantage of MPEC solvers for computations. Numerical examples illustrate the theoretical developments and demonstrate the efficiency of the approach.