Find the value of $\int_{-\infty}^\infty \int_{-\infty}^\infty e^{-(x^2+xy+y^2)} \, dx\,dy$ 
Given that $\int_{-\infty}^\infty e^{-x^2} \, dx=\sqrt{\pi}$. Find the value of $$\int_{-\infty}^\infty\int_{-\infty}^\infty e^{-(x^2+xy+y^2)} \, dx\,dy$$

I don't understand how I find this double integral by using the given data. Please help.

For fun, I want to point out that much more can be said. The spectral theorem for real symmetric matrices tells us that real symmetric matrices are orthogonally diagonalizable. Thus, if $A$ is some symmetric matrix, then there exists an orthogonal matrix $U$ and diagonal matrix $D$ such that $A$ can be written as $A=U^{-1} DU$. Hence $x^TAx=x^TU^TDUx=(Ux)^TD(Ux)$. The Jacobian det of the transformation $x\mapsto Ux$ is simply $1$, so we have a change of variables $u=Ux$:
$$\begin{align} \int_{{\bf R}^n}\exp\left(-x^TAx\right)\,dV & =\int_{{\bf R}^n}\exp\left(-(\lambda_1u_1^2+\cdots+\lambda_nu_n^2)\right)\,dV \\[6pt]
& =\prod_{i=1}^n\int_{-\infty}^{+\infty}\exp\left(-\lambda_i u_i^2\right)du_i \\[6pt]
& =\prod_{i=1}^n\left[\frac{1}{\sqrt{\lambda_i}}\int_{-\infty}^{+\infty}e^{-u^2} \, du\right] \\[6pt]
& =\sqrt{\frac{\pi^n}{\det A}}.\end{align}$$
Note that we don't even have to compute $U$ or $D$. The above calculation is the generalized version of the "completing the squares" approach when $n=2$ (which is invoked elsewhere in this thread).
This formula is in fact the basis for the Feynman path integral formulation of the functional determinant from quantum field theory; since technically the integral diverges we need to compare them instead of looking at individual ones outright. Wikipedia has some more details.
