Interpolant defined on irregularly spaced grid called on 3D matrix

I am trying to define a 9D interpolant on an irregularly spaced grid. In particular, I am defining the interpolant as follows:

itp = interpolate((x1, x2, x3, x4, x5, x6, x7, x8, x9), y, Gridded(Linear()));
where x1, x2, x3, x7, x8 and x9 have 3 points each, while x4, x5 and x6 have 2 points each. This creates a 3 * 3 * 3 * 2 * 2 * 2 * 3 * 3 * 3 interpolant itp.

The problem is that I need to call itp on 5832 * 120 * 7 values for each of the 9 inputs, which is taking a lot of time and memory when I use for loops (over each dimension) to call the interpolant. In Matlab, I am able to pass the 3D array for each of the 9 inputs, and it returns a 5832 by 120 by 7 matrix in 1.4 seconds. What is a quick and easy way to do the same in Julia. Using for loops for each dimension took approximately 21 second in Julia.