I've also stucked with such problem. Since this is old question, but I didn't find the full answer on it, I'll write down my attemption. It doesn't represent the full solution; however, I think that it almost gives the answer.
The density of states $\rho (E)$ is the imaginary part of the self-energy $\Sigma (\mathbf r , \mathbf r, E+i\epsilon)$, where $\epsilon \to 0^{+}$:
$$
\rho (E) = -\frac{1}{\pi}\text{Im}\left(\lim_{\epsilon \to 0^{+}}\Sigma (\mathbf r , \mathbf r , E + i\epsilon) \right)
$$
How can we determine $\Sigma (\mathbf r , \mathbf r , E + i\epsilon)$?
By the definition, the Green operator is
$$
\hat{G} (T) = \frac{1}{T\hat{I} - \hat{H}} \equiv \sum_{\mathbf k}\frac{|\mathbf k \rangle \langle \mathbf k|}{T - E(\mathbf k)}, \quad T \equiv E + i\epsilon
$$
Next, the Green function which connects the site $l$ of the lattice with itself (which is exactly the self energy) is
$$
\tag 1 \Sigma (T, l, l) = \sum_{\mathbf k}\frac{\langle \mathbf{l}|\mathbf k \rangle \langle \mathbf k| \mathbf{l}\rangle}{T - E(\mathbf k)}
$$
Lets talk about graphene in the nearest neigbours approximation. Its lattice is hexagonal (honeycomb), which can be represented by two interpenetrating triangular lattices with the strength of interaction given by $t$. Only the nearest cites (say, $A$ and $B$) of these lattices interact, so $\hat{H}$ lives in the space which it the direct product of spaces of two triangular lattices. Now this results in the fact that the hamiltonian can be given in the form of sum of two-dimensional matrix. Thus, for given cite the denominator of $(1)$ is
$$
\tag 2 \begin{pmatrix} T & -\mu t \\ -\mu^{*}t & T\end{pmatrix}
$$
Here $\mu$ defines the character of lattice, being
$$
\mu = e^{ik_{x}a} + e^{i\left(\frac{\sqrt{3}k_{y}a}{2} - \frac{k_{x}a}{2} \right)} + e^{-i\left(\frac{\sqrt{3}k_{y}a}{2} + \frac{k_{x}a}{2} \right)}
$$
Substituting $(2)$ into $(1)$, you can convert $(1)$ to the form
$$
\tag 3 \rho (E) = -\frac{1}{\pi}\text{Im}\left[ \lim_{\epsilon \to 0^{+}}\int \limits_{\text{1st Br. zone}}\frac{d^{2}\mathbf k}{(2 \pi)^{2}}\frac{T}{T^{2} - t^{2}|\mu|^{2}}\right] = -\frac{E}{8t^{2}\pi}\text{Im} \left[\lim_{\epsilon \to 0^{+}}\tilde{G}\left(\frac{T}{t}\right)\right],
$$
where
$$
\tilde{G}\left(\frac{T}{t}\right) \equiv \frac{1}{\pi^2}\int \limits_{-\pi}^{\pi}\frac{dxdy}{
\frac{\frac{T^{2}}{t^{2}} - 3}{2} - \cos(2y) - 2\cos (y)cos(3x)}
$$
Such quantity can be computed (there is no derivation of this result) for $\frac{\frac{T^{2}}{t^{2}} - 3}{2} > 3$, and
$$
\tag 4 \tilde{G}\left(\frac{T}{t}\right) = \frac{T}{t\pi}\frac{1}{\sqrt{\left(\frac{T}{t} - 1\right)^3\left(\frac{T}{t} + 3\right)}}K\left( \frac{4\sqrt{\frac{T}{t}}}{\sqrt{\left(\frac{T}{t} - 1\right)^3\left(\frac{T}{t} + 3\right)}}\right),
$$
where $K(x)$ is the elliptic integral of the first kind:
$$
K(x) = \int \limits_{0}^{\frac{\pi}{2}}\frac{dy}{\sqrt{1 - x^{2}\sin^{2}(y)}}
$$
The only thing which you have to do is to compute analytic continuation of $(4)$ and then to compute its imaginary part multiplied by four (which corresponds to the degeneracy of spins and two sites).
An edit
Here is the full derivation of the density of states in graphene.
