Support vector machine training can be represented as a large quadratic program. We present an efficient and numerically stable algorithm for this problem using primal-dual interior point methods. Reformulating the problem to exploit separability of the Hessian eliminates the main source of computational complexity, resulting in an algorithm which requires only O(n) operations per iteration. Extensive use of L3 BLAS functions enables good parallel efficiency on shared-memory processors. As the algorithm works in primal and dual spaces simultaneously, our approach has the advantage of obtaining the hyperplane weights and bias directly from the solver.