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 problems to a desired accuracy efficiently, we develop a two-phase proximal augmented Lagrangian method, 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 our proposed algorithm against the highly optimized commercial solver Gurobi and the open source solver OSQP are presented to demonstrate the high efficiency and robustness of our proposed algorithm for solving various classes of large-scale convex QP problems.