We consider option pricing problems in the stochastic volatility jump diffusion model with correlated and contemporaneous jumps (SVCJ) in both the return and variance processes. The option value function solves a partial integro-differential equation (PIDE). We discretize this PIDE in space by the quadratic finite element (FE) method and integrate the resulting ordinary differential equation in time using an implicit–explicit Euler-based extrapolation scheme. The coefficient matrix of the resulting linear systems is block pentadiagonal with pentadiagonal blocks. The preconditioned biconjugate gradient stabilized (PBiCGSTAB) method is used to solve the linear systems. According to the structure of the coefficient matrix, several preconditioners are implemented and compared. The performance of preconditioning techniques for solving block tridiagonal systems resulting from the linear FE discretization of the PIDE is also investigated. The combination of the quadratic FE for spatial discretization, the extrapolation scheme for time discretization and the PBiCGSTAB method with an appropriate preconditioner is found to be very efficient for solving the option pricing problems in the SVCJ model. Compared with the standard second-order linear FE method combined with the popular successive over-relaxation (SOR) linear system solver, the proposed method reduces computational time by about a factor of 20 at the accuracy level of 1% and by more than fifty times at the accuracy level of 0.1% for the barrier option example tested in the paper.