In this paper, we aim to solve high dimensional convex quadratic programming (QP) problems with a large number of quadratic terms, linear equality and inequality constraints. In order to solve the targeted {\bf QP} problems to a desired accuracy efficiently, we develop a two-phase {\bf P}roximal {\bf A}ugmented {\bf L}agrangian method {(QPPAL)}, with Phase I to generate a reasonably good initial point to warm start Phase II to obtain an accurate solution efficiently. More specifically, in Phase I, based on the recently developed symmetric Gauss-Seidel (sGS) decomposition technique, we design a novel sGS based semi-proximal augmented Lagrangian method for the purpose of finding a solution of low to medium accuracy. Then, in Phase II, a proximal augmented Lagrangian algorithm is proposed to obtain a more accurate solution efficiently. Extensive numerical results evaluating the performance of {QPPAL} against {existing state-of-the-art solvers Gurobi, OSQP and QPALM} are presented to demonstrate the high efficiency and robustness of our proposed algorithm for solving various classes of large-scale convex QP problems. {The MATLAB implementation of the software package QPPAL is available at: https://blog.nus.edu.sg/mattohkc/softwares/qppal/