  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:
\url{https://blog.nus.edu.sg/mattohkc/softwares/qppal/}.
