This paper introduces an efficient algorithm for computing the best approximation of a given matrix onto the intersection of linear equalities, inequalities and the doubly nonnegative cone (the cone of all positive semidefinite matrices whose elements are nonnegative). In contrast to directly applying the block coordinate descent type methods, we propose an inexact accelerated (two-)block coordinate descent algorithm to tackle the four-block unconstrained nonsmooth dual program. The proposed algorithm hinges on the efficient semismooth Newton method to solve the subproblems, which have no closed form solutions since the original four blocks are merged into two larger blocks. The $O(1/k^2)$ iteration complexity of the proposed algorithm is established. Extensive numerical results over various large scale semidefinite programming instances from relaxations of combinatorial problems demonstrate the effectiveness of the proposed algorithm.