We propose an efficient computational method for linearly constrained quadratic optimization problems (QOPs) with complementarity constraints based on their Lagrangian and doubly nonnegative (DNN) relaxation and first-order algorithms. The simplified Lagrangian-CPP relaxation of such QOPs proposed by Arima, Kim, and Kojima in 2012 takes one of the simplest forms, an unconstrained conic linear optimization problem with a single Lagrangian parameter in a completely positive (CPP) matrix variable with its upper-left element fixed to 1. Replacing the CPP matrix variable by a DNN matrix variable, we derive the Lagrangian-DNN relaxation, and establish the equivalence between the optimal value of the DNN relaxation of the original QOP and that of the Lagrangian-DNN relaxation. We then propose an efficient numerical method for the Lagrangian-DNN relaxation using a bisection method combined with the proximal alternating direction multiplier and the accelerated proximal gradient methods. Numerical results on binary QOPs, quadratic multiple knapsack problems, maximum stable set problems, and quadratic assignment problems illustrate the superior performance of the proposed method for attaining tight lower bounds in shorter computational time.