I would like to evaluate $$ \sum_{q_1 = -\infty}^{\infty} \cdots \sum_{q_N = -\infty}^{\infty} e^{-\sum_{j}\sum_{k} q_{k} A_{kj} q_{j}} $$ with $A$ a real $N\times N$ symmetric matrix.
I know how to compute this when $q$ is continuous (the sum is an integral), and I know how to compute this when $A$ is a scalar (a $1\times 1$, this leads to the Jacobi theta), but If I try to diagonalize $A$, I end up with a transformation for the $q$s that I don't know how to write down in terms of a sum. The transformed $q$s, $q' = O^{T}q$, are linear combinations of the $q$s, and I don't know what the analogous Jacobian-like object would be for summation (in place of integration). Thanks.
This cannot be expressed in terms of elementary (or Jacobi theta) functions: in fact the series $$\Theta\left(\mathbf z | \Omega\right)=\sum_{\mathbf q\in\mathbb Z^N}e^{\pi i \mathbf q\cdot \Omega\cdot \mathbf q+2\pi i \mathbf q \cdot \mathbf z}$$ is a multidimensional generalization of the Jacobi theta function called Riemann theta function. Your case corresponds to setting $\mathbf z=\mathbf 0$, $\Omega=\frac{iA}{\pi}$, i.e. to Riemann theta constants.
For the numerical evaluation, you may have a look at this paper.