Improving speed of multidimensional Gaussian quadrature

I am trying to do a numerical integration using Gauss-Legendre quadrature over rectangular regions. The equations used in the code below is given as:

$$Q_{1}= \int_{c_{i}}^{d_{i}}\int_{a_{i}}^{b_{i}} I_{fn}(x_{o}, y_{o},x_{i},y_{i},x_{p},y_{p},z_{p},z_{o},z_{d},p_D,p_A) \ \mathrm{d} x_{i} \ \mathrm{d} y_{i} \\ \approx \frac{(d_i - c_i)}{2} \frac{(b_i - a_i)}{2} \sum_{u=1}^{n_{i}} \sum_{v=1}^{n_{i}}w_uw_v \ I_{fn}(x_{o}, y_{o},\frac{(a_{i} + b_{i})}{2} + \frac{(b_{i} - a_{i})}{2}x_{u},\frac{(c_{i} + d_{i})}{2} + \frac{(d_{i} - c_{i})}{2}y_{v},x_{p},y_{p},z_{p},z_{o},z_{d},p_D,p_A) \\ Q_{2}= \int_{c_{o}}^{d_{o}}\int_{a_{o}}^{b_{o}} Q_{1}(x_{o}, y_{o},x_{p},y_{p},z_{p},z_{o},z_{d},p_D,p_A) \ \mathrm{d} x_{o} \ \mathrm{d} y_{o} \\ \approx \frac{(d_o - c_o)}{2} \frac{(b_o - a_o)}{2} \sum_{u=1}^{n_{o}} \sum_{v=1}^{n_{o}}w_uw_v \ Q_1(\frac{(a_{o} + b_{o})}{2} + \frac{(b_{o} - a_{o})}{2}x_{u},\frac{(c_{o} + d_{o})}{2} + \frac{(d_{o} - c_{o})}{2}y_{v},x_{p},y_{p},z_{p},z_{o},z_{d},p_D,p_A)$$

My code segment is given here:

Needs["NumericalDifferentialEquationAnalysis"];
Remove[Result, GaussWeightsRin, GaussWeightsRout, FirstTerm, Ifn, Q1, Q2]

FirstTerm[xi_, yi_, xp_, yp_, zp_, Zo_, Zd_, pD_, pA_] :=
Piecewise[{{(Abs[Zo - zp Cos[pA] + Zd Cos[pA - pD] + xp Sin[pA] - xi Sin[pA + ArcTan[-Cos[pD], Sin[pD]]]]^3 (Zd + Zo Cos[pA - pD] - zp Cos[pD] + xp Sin[pD]) (Zd Cos[pA - pD] - xi Sin[pA + ArcTan[-Cos[pD], Sin[pD]]]))/(Sqrt[Abs[Cos[pD]]^2 +  Abs[Sin[pD]]^2] (Abs[yi - yp]^2 +  Abs[zp - Zo Cos[pA] - Zd Cos[pD] + xi Sin[pD]]^2 +  Abs[xp + xi Cos[pD] + Zo Sin[pA] + Zd Sin[pD]]^2)^(3/2) Abs[Zd Cos[pA - pD] - xi Sin[pA + ArcTan[-Cos[pD], Sin[pD]]]]^3 (Zo - zp Cos[pA] + Zd Cos[pA - pD] + xp Sin[pA] - xi Sin[pA + ArcTan[-Cos[pD], Sin[pD]]])),
Cos[pD] != 0 || Sin[pD] != 0}}, ComplexInfinity];

Ifn = Compile[{{xo, _Real, 0}, {yo, _Real, 0}, {xi, _Real, 0}, {yi, _Real, 0},
{xp, _Real, 0}, {yp, _Real, 0}, {zp, _Real, 0}, {Zo, _Real, 0},
{Zd, _Real, 0}, {pD, _Real, 0}, {pA, _Real, 0}},

FirstTerm[xi, yi, xp, yp, zp, Zo, Zd, pD, pA] 2 Exp[-(((xo - xi - 0)/5)^2 + ((yo - yi - 0)/5)^2)],

RuntimeAttributes -> {Listable}, CompilationTarget -> "C", Parallelization -> True, RuntimeOptions -> "Speed"];

With[{Ifn = Ifn},
Q1 = Compile[{{xo, _Real, 0}, {yo, _Real, 0}, {ai, _Real, 0}, {bi, _Real, 0},
{ci, _Real, 0}, {di, _Real, 0}, {ni, _Real, 0}, {xp, _Real, 0},
{yp, _Real, 0}, {zp, _Real, 0}, {Zo, _Real, 0}, {Zd, _Real, 0},
{pD, _Real, 0}, {pA, _Real, 0}},

(di - ci)/2*(bi - ai)/2*
Sum[Sum[GaussWeightsRin[[u, 2]] GaussWeightsRin[[v, 2]]
Ifn[xo, yo, (ai + bi)/ 2 + (bi - ai)/2 GaussWeightsRin[[u, 1]], (ci + di)/ 2 + (di - ci)/2 GaussWeightsRin[[v, 1]], xp, yp, zp, Zo, Zd, pD, pA], {u, 1, ni}], {v, 1, ni}],

RuntimeAttributes -> {Listable}, CompilationTarget -> "C", Parallelization -> True, RuntimeOptions -> "Speed"]];

With[{Q1 = Q1},
Q2 = Compile[{{ma, _Integer, 0}, {mb, _Integer, 0}, {no, _Integer, 0}, {ni, _Integer, 0},
{xp, _Real, 0}, {yp, _Real, 0}, {zp, _Real, 0}, {Zo, _Real, 0},
{Zd, _Real, 0}, {pD, _Real, 0}, {pA, _Real, 0}},

Module[{ao = ma*0.048046875 (1 - 1/(2 Abs[ma])) - 0.048/2,
bo = ma*0.048046875 (1 - 1/(2 Abs[ma])) + 0.048/2,
co = mb*0.048046875 (1 - 1/(2 Abs[mb])) - 0.048/2,
do = mb*0.048046875 (1 - 1/(2 Abs[mb])) + 0.048/2,
ai = ma*0.048046875 (1 - 1/(2 Abs[ma])) - 8,
bi = ma*0.048046875 (1 - 1/(2 Abs[ma])) + 8,
ci = mb*0.048046875 (1 - 1/(2 Abs[mb])) - 8,
di = mb*0.048046875 (1 - 1/(2 Abs[mb])) + 8 },

(do - co)/2*(bo - ao)/2*
Sum[Sum[GaussWeightsRout[[u, 2]] GaussWeightsRout[[v, 2]] Q1[(ao + bo)/2 + (bo - ao)/2 GaussWeightsRout[[u, 1]], (co + do)/2 + (do - co)/2 GaussWeightsRout[[v, 1]], ai, bi, ci, di, ni, xp, yp, zp, Zo, Zd, pD, pA], {u, 1, no}], {v, 1, no}]],

RuntimeAttributes -> {Listable}, CompilationTarget -> "C", Parallelization -> True, RuntimeOptions -> "Speed"]];

With[{Q2 = Q2},
Result = Compile[{{M, _Integer, 1}, {no, _Integer, 0}, {ni, _Integer, 0}, {xp, _Real, 0}, {yp, _Real, 0}, {zp, _Real, 0}, {Zo, _Real, 0}, {Zd, _Real, 0}, {pD, _Real, 0}, {pA, _Real, 0}},

Q2[CompileGetElement[M, 1], CompileGetElement[M, 2], no, ni, xp, yp, zp, Zo, Zd, pD, pA],

RuntimeAttributes -> {Listable}, CompilationTarget -> "C", Parallelization -> True, RuntimeOptions -> "Speed"]];


I want to run this code (Result[]) for a set of different points {ma,mb} and for that, I use Outer.

Num = 10;
CordsX = DeleteCases[Range[Num/2, -Num/2, -1], 0, Infinity];
CordsY = {1};
Outer[Result[{#2, #1}, 2, 15, 10, 10, 10, 100, 200, \[Pi], 0] &, CordsY, CordsX, 1] // AbsoluteTiming


This is very slow for my application as I need to run this over more points (i.e. for higher values of Num). Have I correctly used compile in my code? and how do I speed up my code?.

• (1) You might want to take in mathematica.stackexchange.com/questions/24595/…. (2) GaussianQuadratureWeights is undefined. Are you missing Needs@"NumericalDifferentialEquationAnalysis"? – Michael E2 Jun 12 '18 at 23:20
• – Michael E2 Jun 12 '18 at 23:24
• So I need to use pure functions and use "InlineCompiledFunctions" -> True wherever I want to use compiled functions?. (2) I forgot to write Needs@"NumericalDifferentialEquationAnalysis" in the code above. – dykes Jun 13 '18 at 8:54