Can NIntegrate be used with the Levin method in several dimensions?

I've got some data in the form of an interpolating function. It's a function of three variables, $\rho(x,y,z)$. I'd basically like to integrate this with some phase over a cube of known size, like $$\iiint_{cube}\rho(x,y,z) e^{i(x+y+2z)}\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z$$ I've been trying to do this with NIntegrate, but I'm having a lot of difficulty, probably because of the exponential part's oscillation. Many of the integration routines don't converge. Note that I don't know for certain that $\rho$ changes slowly in comparison to the exponential bit. I know that Mathematica has an algorithm for dealing with oscillating integrands, called the Levin Rule, but I haven't been able to figure out how to work that in more than one dimension numerically. I was wondering if someone could either give insight on whether or not this can be done or point me to some references on multidimensional Levin integration using Mathematica. I've read the help pages on NIntegrate integration rules, and it hasn't born fruit.

To be more concrete, here's the most functional code I've got so far. It always says it fails to converge, but it DOES spit out a number that's within 10% of the true answer. I'd like to get something that's more predictably accurate. Speed is somewhat secondary, but still important. I've uploaded rho.mx here.

rho = Import["rho.mx"];
Clear[Sg]
Sg = Compile[{},
NIntegrate[rho[x, y, z]*Exp[-I (x + y + 2 z)], {x, 0, 2}, {y,
0, 2}, {z, 0, 2}, Method -> "QuasiMonteCarlo"]];

• You can just use Method -> "LevinRule" for single or multiple integrals. (It might be using it automatically anyway.) Here's an overview for the strategies, rules and their uses. – Michael E2 Oct 22 '14 at 4:11