# Computation performance involving multiple numeric integrals

I have the following code that has been running for more than three days without an answer, unfortunately. Now, my institution provides us with a computer with Intel Xeon CPU E-52650 v3 @ 20 CPU, 2.3GHz and 32 GB RAM. I have tried the same code on this computer, and it is taking the same time and still running without getting any answer.

I heard that we can use Needs["CUDALink"] but it gave us False. I don't know now how I can benefit from the multi-core computer to speed up our computations or even how to rewrite the code in different way to help in the process. Is there any package we have to download to benefit from the supercomputer? Also, I saw online some info regarding parallel computing but I really don't understand it. Can any one help me speed up this computation?

Here is the code I wrote to evaluate my desired result, ae:

psi1[x_] :=1/24 (-4 (-2 + x)^3 Sign[-2 + x] + (-3 + 2 x)^3 Sign[-(3/2) + x] - 8 Sign[-1 + x] + 24 x Sign[-1 + x] - 24 x^2 Sign[-1 + x] +
8 x^3 Sign[-1 + x] + 3 Sign[-(1/2) + x] -
18 x Sign[-(1/2) + x] + 36 x^2 Sign[-(1/2) + x] -
24 x^3 Sign[-(1/2) + x] + 3 Sign[1/2 + x] + 18 x Sign[1/2 + x] +
36 x^2 Sign[1/2 + x] + 24 x^3 Sign[1/2 + x] - 8 Sign[1 + x] -
24 x Sign[1 + x] - 24 x^2 Sign[1 + x] - 8 x^3 Sign[1 + x] -
27 Sign[3/2 + x] - 54 x Sign[3/2 + x] - 36 x^2 Sign[3/2 + x] -
8 x^3 Sign[3/2 + x] + 32 Sign[2 + x] + 48 x Sign[2 + x] +
24 x^2 Sign[2 + x] + 4 x^3 Sign[2 + x]);

psi2[x_] := 1/(4 Sqrt[
6]) (Sign[1/2 - x] + 2 (-2 + x)^3 Sign[-2 + x] +
27 Sign[-(3/2) + x] - 54 x Sign[-(3/2) + x] +
36 x^2 Sign[-(3/2) + x] - 8 x^3 Sign[-(3/2) + x] -
8 Sign[-1 + x] + 24 x Sign[-1 + x] - 24 x^2 Sign[-1 + x] +
8 x^3 Sign[-1 + x] + 6 x Sign[-(1/2) + x] -
12 x^2 Sign[-(1/2) + x] + 8 x^3 Sign[-(1/2) + x] -
20 x^3 Sign[x] + Sign[1/2 + x] + 6 x Sign[1/2 + x] +
12 x^2 Sign[1/2 + x] + 8 x^3 Sign[1/2 + x] + 8 Sign[1 + x] +
24 x Sign[1 + x] + 24 x^2 Sign[1 + x] + 8 x^3 Sign[1 + x] -
27 Sign[3/2 + x] - 54 x Sign[3/2 + x] - 36 x^2 Sign[3/2 + x] -
8 x^3 Sign[3/2 + x] + 16 Sign[2 + x] + 24 x Sign[2 + x] +
12 x^2 Sign[2 + x] + 2 x^3 Sign[2 + x]);

psi3[x_] :=1/24 (-4 (-2 + x)^3 Sign[-2 + x] +
3 (-3 + 2 x)^3 Sign[-(3/2) + x] + 56 Sign[-1 + x] -
168 x Sign[-1 + x] + 168 x^2 Sign[-1 + x] -
56 x^3 Sign[-1 + x] - 7 Sign[-(1/2) + x] +
42 x Sign[-(1/2) + x] - 84 x^2 Sign[-(1/2) + x] +
56 x^3 Sign[-(1/2) + x] - 7 Sign[1/2 + x] - 42 x Sign[1/2 + x] -
84 x^2 Sign[1/2 + x] - 56 x^3 Sign[1/2 + x] + 56 Sign[1 + x] +
168 x Sign[1 + x] + 168 x^2 Sign[1 + x] + 56 x^3 Sign[1 + x] -
81 Sign[3/2 + x] - 162 x Sign[3/2 + x] - 108 x^2 Sign[3/2 + x] -
24 x^3 Sign[3/2 + x] + 32 Sign[2 + x] + 48 x Sign[2 + x] +
24 x^2 Sign[2 + x] + 4 x^3 Sign[2 + x]);

psi4[x_] :=1/12 (-(-2 + x)^3 Sign[-2 + x] + (-3 + 2 x)^3 Sign[-(3/2) + x] +
28 Sign[-1 + x] - 84 x Sign[-1 + x] + 84 x^2 Sign[-1 + x] -
28 x^3 Sign[-1 + x] - 7 Sign[-(1/2) + x] +
42 x Sign[-(1/2) + x] - 84 x^2 Sign[-(1/2) + x] +
56 x^3 Sign[-(1/2) + x] - 70 x^3 Sign[x] + 7 Sign[1/2 + x] +
42 x Sign[1/2 + x] + 84 x^2 Sign[1/2 + x] +
56 x^3 Sign[1/2 + x] - 28 Sign[1 + x] - 84 x Sign[1 + x] -
84 x^2 Sign[1 + x] - 28 x^3 Sign[1 + x] + 27 Sign[3/2 + x] +
54 x Sign[3/2 + x] + 36 x^2 Sign[3/2 + x] +
8 x^3 Sign[3/2 + x] - 8 Sign[2 + x] - 12 x Sign[2 + x] -
6 x^2 Sign[2 + x] - x^3 Sign[2 + x]);

psi1jk[x_, j_, k_] :=Piecewise[{{(Sqrt[2])^j psi1[2^j x - k],
0 <= j}, {2^j psi1[2^j (x - k)], j < 0}}];
psi2jk[x_, j_, k_] :=Piecewise[{{(Sqrt[2])^j psi2[2^j x - k],
0 <= j}, {2^j psi2[2^j (x - k)], j < 0}}];

psi3jk[x_, j_, k_] :=Piecewise[{{(Sqrt[2])^j psi3[2^j x - k],
0 <= j}, {2^j psi3[2^j (x - k)], j < 0}}];

psi4jk[x_, j_, k_] :=Piecewise[{{(Sqrt[2])^j psi4[2^j x - k],
0 <= j}, {2^j psi4[2^j (x - k)], j < 0}}];

PSI[j_, k_, l_, s_] :=NIntegrate[
psi1jk[x, j,
k]*(psi1jk[x, l, s] + psi2jk[x, l, s] + psi3jk[x, l, s] +
psi4jk[x, l, s]), {x, -1, 1}] -
NIntegrate[(x t + (x t)^2)*
psi1jk[t, j,
k]*(psi1jk[x, l, s] + psi2jk[x, l, s] + psi3jk[x, l, s] +
psi4jk[x, l, s]), {x, -1, 1}, {t, -1, 1}] +
NIntegrate[
psi2jk[x, j,
k]*(psi1jk[x, l, s] + psi2jk[x, l, s] + psi3jk[x, l, s] +
psi4jk[x, l, s]), {x, -1, 1}] -
NIntegrate[(x t + (x t)^2)*
psi2jk[t, j,
k]*(psi1jk[x, l, s] + psi2jk[x, l, s] + psi3jk[x, l, s] +
psi4jk[x, l, s]), {x, -1, 1}, {t, -1, 1}] +
NIntegrate[
psi3jk[x, j,
k]*(psi1jk[x, l, s] + psi2jk[x, l, s] + psi3jk[x, l, s] +
psi4jk[x, l, s]), {x, -1, 1}] -
NIntegrate[(x t + (x t)^2)*
psi3jk[t, j,
k]*(psi1jk[x, l, s] + psi2jk[x, l, s] + psi3jk[x, l, s] +
psi4jk[x, l, s]), {x, -1, 1}, {t, -1, 1}] +
NIntegrate[
psi4jk[x, j,
k]*(psi1jk[x, l, s] + psi2jk[x, l, s] + psi3jk[x, l, s] +
psi4jk[x, l, s]), {x, -1, 1}] -
NIntegrate[(x t + (x t)^2)*
psi4jk[t, j,
k]*(psi1jk[x, l, s] + psi2jk[x, l, s] + psi3jk[x, l, s] +
psi4jk[x, l, s]), {x, -1, 1}, {t, -1, 1}];

ce = ArrayReshape[Table[PSI[j, k, l, s], {j, -4., 4.}, {k, -32., 17.}, {l, -4.,4.}, {s, -32., 17.}], {450, 450}];

coef = Flatten[Table[NIntegrate[(psi1jk[x, l, s] + psi2jk[x, l, s] +
psi3jk[x, l, s] + psi4jk[x, l, s]), {x, -1, 1}], {l, -4.,4.}, {s, -32., 17.}]];

ae = PseudoInverse[ce].coef


• Maybe you'd like to say some words about what you are doing there? Are all those Signs really necessary? Why do all the Signs pop up? Are they necessary? They might be a major obstruction for using more efficient integration schemes. Does the matrix ce in the end happen to contain many zeroes? Then most of the integrals are superfluous. – Henrik Schumacher May 29 '18 at 19:32
• And forget about Cuda for the moment. This code won't run on a GPU without major refactorization. – Henrik Schumacher May 29 '18 at 19:33
• @HenrikSchumacher I am working on some simulation that used wavelets generated by B-splines. Sign functions there after I did the Fourier transform to get the form of B-splines in time domain. We can't remove all of them as it is essential to get the domain of each spline. What do you mean for using efficient integration schemes? – Mutaz May 29 '18 at 19:57
• Btw.: A 20 core machine is not a "supercomputer" nowerdays. And your code employs no parallelization, so it runs only on one core at the moment. I guess you can get rid of almost all integrations but I don't see any meaning in the integrals. This looks a bit as if wavelets were involved... – Henrik Schumacher May 29 '18 at 19:57
• @HenrikSchumacher The matrix ce yes has many zeros and may have multiple rows with a zero values, but I cant ignore them as I have rows after them that are not zero. If we skip them we will get bad approximation. So I don't know if there is a way to skip them and consider those nonzero without affecting evaluating the vector ae. – Mutaz May 29 '18 at 19:57

Start with a simpler problem. I you change the iterators below (after With[ ) gradually you will get an idea how the complexity of your problem behaves. If there are no surprises with NIntegrate your 405000 calls to NIntegrate should be doable on your multicore machine during a night or so.

    psi1[x_] := (1/24)*(-4*(-2 + x)^3*Sign[-2 + x] + (-3 + 2*x)^3*Sign[-(3/2) + x] - 8*Sign[-1 + x] +
24*x*Sign[-1 + x] - 24*x^2*Sign[-1 + x] + 8*x^3*Sign[-1 + x] + 3*Sign[-(1/2) + x] -
18*x*Sign[-(1/2) + x] + 36*x^2*Sign[-(1/2) + x] - 24*x^3*Sign[-(1/2) + x] + 3*Sign[1/2 + x] +
18*x*Sign[1/2 + x] + 36*x^2*Sign[1/2 + x] + 24*x^3*Sign[1/2 + x] - 8*Sign[1 + x] - 24*x*Sign[1 + x] -
24*x^2*Sign[1 + x] - 8*x^3*Sign[1 + x] - 27*Sign[3/2 + x] - 54*x*Sign[3/2 + x] - 36*x^2*Sign[3/2 + x] -
8*x^3*Sign[3/2 + x] + 32*Sign[2 + x] + 48*x*Sign[2 + x] + 24*x^2*Sign[2 + x] + 4*x^3*Sign[2 + x]);
psi2[x_] := (1*(Sign[1/2 - x] + 2*(-2 + x)^3*Sign[-2 + x] + 27*Sign[-(3/2) + x] - 54*x*Sign[-(3/2) + x] +
36*x^2*Sign[-(3/2) + x] - 8*x^3*Sign[-(3/2) + x] - 8*Sign[-1 + x] + 24*x*Sign[-1 + x] -
24*x^2*Sign[-1 + x] + 8*x^3*Sign[-1 + x] + 6*x*Sign[-(1/2) + x] - 12*x^2*Sign[-(1/2) + x] +
8*x^3*Sign[-(1/2) + x] - 20*x^3*Sign[x] + Sign[1/2 + x] + 6*x*Sign[1/2 + x] + 12*x^2*Sign[1/2 + x] +
8*x^3*Sign[1/2 + x] + 8*Sign[1 + x] + 24*x*Sign[1 + x] + 24*x^2*Sign[1 + x] + 8*x^3*Sign[1 + x] -
27*Sign[3/2 + x] - 54*x*Sign[3/2 + x] - 36*x^2*Sign[3/2 + x] - 8*x^3*Sign[3/2 + x] + 16*Sign[2 + x] +
24*x*Sign[2 + x] + 12*x^2*Sign[2 + x] + 2*x^3*Sign[2 + x]))/(4*Sqrt[6]);
psi3[x_] := (1*(-4*(-2 + x)^3*Sign[-2 + x] + 3*(-3 + 2*x)^3*Sign[-(3/2) + x] + 56*Sign[-1 + x] -
168*x*Sign[-1 + x] + 168*x^2*Sign[-1 + x] - 56*x^3*Sign[-1 + x] - 7*Sign[-(1/2) + x] +
42*x*Sign[-(1/2) + x] - 84*x^2*Sign[-(1/2) + x] + 56*x^3*Sign[-(1/2) + x] - 7*Sign[1/2 + x] -
42*x*Sign[1/2 + x] - 84*x^2*Sign[1/2 + x] - 56*x^3*Sign[1/2 + x] + 56*Sign[1 + x] + 168*x*Sign[1 + x] +
168*x^2*Sign[1 + x] + 56*x^3*Sign[1 + x] - 81*Sign[3/2 + x] - 162*x*Sign[3/2 + x] -
108*x^2*Sign[3/2 + x] - 24*x^3*Sign[3/2 + x] + 32*Sign[2 + x] + 48*x*Sign[2 + x] + 24*x^2*Sign[2 + x] +
4*x^3*Sign[2 + x]))/24;
psi4[x_] := (1*(-((-2 + x)^3*Sign[-2 + x]) + (-3 + 2*x)^3*Sign[-(3/2) + x] + 28*Sign[-1 + x] -
84*x*Sign[-1 + x] + 84*x^2*Sign[-1 + x] - 28*x^3*Sign[-1 + x] - 7*Sign[-(1/2) + x] +
42*x*Sign[-(1/2) + x] - 84*x^2*Sign[-(1/2) + x] + 56*x^3*Sign[-(1/2) + x] - 70*x^3*Sign[x] +
7*Sign[1/2 + x] + 42*x*Sign[1/2 + x] + 84*x^2*Sign[1/2 + x] + 56*x^3*Sign[1/2 + x] - 28*Sign[1 + x] -
84*x*Sign[1 + x] - 84*x^2*Sign[1 + x] - 28*x^3*Sign[1 + x] + 27*Sign[3/2 + x] + 54*x*Sign[3/2 + x] +
36*x^2*Sign[3/2 + x] + 8*x^3*Sign[3/2 + x] - 8*Sign[2 + x] - 12*x*Sign[2 + x] - 6*x^2*Sign[2 + x] -
x^3*Sign[2 + x]))/12;
psi1jk[x_, j_, k_] := Piecewise[{{Sqrt[2]^j*psi1[2^j*x - k], 0 <= j}, {2^j*psi1[2^j*(x - k)], j < 0}}];
psi2jk[x_, j_, k_] := Piecewise[{{Sqrt[2]^j*psi2[2^j*x - k], 0 <= j}, {2^j*psi2[2^j*(x - k)], j < 0}}];
psi3jk[x_, j_, k_] := Piecewise[{{Sqrt[2]^j*psi3[2^j*x - k], 0 <= j}, {2^j*psi3[2^j*(x - k)], j < 0}}];
psi4jk[x_, j_, k_] := Piecewise[{{Sqrt[2]^j*psi4[2^j*x - k], 0 <= j}, {2^j*psi4[2^j*(x - k)], j < 0}}];
PSI[j_, k_, l_, s_] :=NIntegrate[(psi1jk[x, j, k] + psi2jk[x, j, k] + psi3jk[x, j, k] + psi4jk[x, j, k])*
(psi1jk[x, l, s] + psi2jk[x, l, s] + psi3jk[x, l, s] + psi4jk[x, l, s]), {x, -1, 1}]      -
NIntegrate[(x*t + (x*t)^2)*(psi1jk[t, j, k] + psi2jk[t, j, k] + psi3jk[t, j, k] + psi4jk[t, j, k])*
(psi1jk[x, l, s] + psi2jk[x, l, s] + psi3jk[x, l, s] + psi4jk[x, l, s]), {x, -1, 1}, {t, -1, 1}]

With[{
jj = {j, -1, 1}
,
kk = {k, -1, 1}
,
ll = {l, -1, 1}
,
ss = {s, -1, 1}
},
tot = Length[Range @@ Rest[jj]]*Length[Range @@ Rest[kk]]*
Length[Range @@ Rest[ll]]*Length[Range @@ Rest[ss]];
If[Length[Kernels[]] === 0, LaunchKernels[]];
klist = ConstantArray[0, Length[Kernels[]]]; SetSharedVariable[klist];
SetAttributes[parProgress, HoldFirst];
parProgress[ke_, n_] :=
Row[{ProgressIndicator[Dynamic[Total[ke]], {0, n}], " ",
Dynamic[Round[100*(Total[ke]/n)]], " %"}];
Print[parProgress[klist, tot]];
(* Returning the total wall-clock time needed: *)
AbsoluteTiming[
ce = (ArrayReshape[#1, {Sqrt[tot], Sqrt[tot]}] & )[
ParallelTable[klist[[\$KernelID]]++; PSI[j, k, l, s],
jj, kk, ll, ss]
];
]]

• BTW: parProgress is a really nice utility function to monitor parallel progress and I really think WRI should build something like this into the next version. – Rolf Mertig May 29 '18 at 21:19
• Thank you for your email. I run your code and I got Out[1113]= {5.84482, Null}. What is that mean? How this will benefit me to get the vector ae? BTW I was nit able to get ride of these surprises came from NIntegrate – Mutaz May 29 '18 at 21:27
• I put a comment into the code before AbsoluteTiming. I am afraid you have to study my code a bit by yourself and try to understand it. And then go on by yourself. – Rolf Mertig May 29 '18 at 21:40

If you simplify the integrals to

PSI[j_, k_, l_, s_] :=
NIntegrate[(psi1jk[x, j, k] + psi2jk[x, j, k] + psi3jk[x, j, k] +
psi4jk[x, j, k])*(psi1jk[x, l, s] + psi2jk[x, l, s] +
psi3jk[x, l, s] + psi4jk[x, l, s]), {x,-1,1}] -
NIntegrate[(x t + (x t)^2)*(psi1jk[t, j, k] + psi2jk[t, j, k] +
psi3jk[t, j, k] + psi4jk[t, j, k])*(psi1jk[x, l, s] +
psi2jk[x, l, s] + psi3jk[x, l, s] + psi4jk[x, l, s]), {x,-1,1}, {t,-1,1}]


you see that they are symmetric under exchange of (j,k) with (l,s). So if you define

PSI[j_, k_, l_, s_] := PSI[j,k,l,s] = PSI[l,s,j,k] =
NIntegrate[(psi1jk[x, j, k] + psi2jk[x, j, k] + psi3jk[x, j, k] +
psi4jk[x, j, k])*(psi1jk[x, l, s] + psi2jk[x, l, s] +
psi3jk[x, l, s] + psi4jk[x, l, s]), {x,-1,1}] -
NIntegrate[(x t + (x t)^2)*(psi1jk[t, j, k] + psi2jk[t, j, k] +
psi3jk[t, j, k] + psi4jk[t, j, k])*(psi1jk[x, l, s] +
psi2jk[x, l, s] + psi3jk[x, l, s] + psi4jk[x, l, s]), {x,-1,1}, {t,-1,1}]


you save at least a factor of 2, and I'd think the code runs within a few hours on my laptop.

Also, you may gain some speed if you simplify the psi functions like this:

psi1[x_] = Piecewise[{{(2 + x)^3/3, -2 < x < -3/2},
{(5 - 2*x*(3 + 2*x*(3 + x)))/12, -3/2 < x < -1},
{-1/4 - (x*(5 + 2*x*(3 + x)))/2, -1 < x < -1/2},
{x*(-1 + x^2), -1/2 < x < 0 || 0 < x < 1/2},
{1/4 - (x*(5 + 2*(-3 + x)*x))/2, 1/2 < x < 1},
{-5/12 - (x*(3 + 2*(-3 + x)*x))/6, 1 < x < 3/2},
{(-2 + x)^3/3, 3/2 < x < 2}}, 0]
psi2[x_] = Piecewise[{{(2 + x)^3/Sqrt[6], -2 < x < -3/2},
{(-11 - 6*x*(5 + x*(4 + x)))/(2*Sqrt[6]), -3/2 < x < -1},
{(-3 - 6*x + 2*x^3)/(2*Sqrt[6]), -1 < x < -1/2},
{((1 + x)*(-1 + x + 5*x^2))/Sqrt[6], -1/2 < x < 0},
{(-1 + (6 - 5*x)*x^2)/Sqrt[6], 0 < x < 1/2},
{-(3 - 6*x + 2*x^3)/(2*Sqrt[6]), 1/2 < x < 1},
{(-11 + 6*x*(5 + (-4 + x)*x))/(2*Sqrt[6]), 1 < x < 3/2},
{-((-2 + x)^3/Sqrt[6]), 3/2 < x < 2}}, 0]
psi3[x_] = Piecewise[{{(2 + x)^3/3, -2 < x < -3/2},
{-49/12 - (x*(57 + 2*x*(21 + 5*x)))/6, -3/2 < x < -1},
{7/12 + (9*x)/2 + 7*x^2 + 3*x^3, -1 < x < -1/2},
{x - (5*x^3)/3, -1/2 < x < 0 || 0 < x < 1/2},
{-7/12 + (9*x)/2 - 7*x^2 + 3*x^3, 1/2 < x < 1},
{49/12 - (x*(57 + 2*x*(-21 + 5*x)))/6, 1 < x < 3/2},
{(-2 + x)^3/3, 3/2 < x < 2}}, 0]
psi4[x_] = Piecewise[{{-(2 + x)^3/6, -2 < x < -3/2},
{((1 + x)*(19 + x*(23 + 7*x)))/6, -3/2 < x < -1},
{-((1 + x)*(3 + x*(11 + 7*x)))/2, -1 < x < -1/2},
{-1/3 + (5*x^2*(6 + 7*x))/6, -1/2 < x < 0},
{-1/3 + (5*(6 - 7*x)*x^2)/6, 0 < x < 1/2},
{((-1 + x)*(3 + x*(-11 + 7*x)))/2, 1/2 < x < 1},
{-((-1 + x)*(19 + x*(-23 + 7*x)))/6, 1 < x < 3/2},
{(-2 + x)^3/6, 3/2 < x < 2}}, 0]

• Thank you for the note. Could please tell me how did you get ride of the sign function? What command you have used? – Mutaz May 30 '18 at 7:14
• following Akku14's answer: psi1a[x_] = FullSimplify[PiecewiseExpand[psi1[x], Assumptions -> x \[Element] Reals]]` – Roman May 30 '18 at 10:43
• Dear friend, do you know where I can find Mr. Akku14 comments? I suddenly couldnt find it. I think he suggested a nice comments like all of you but it is nit there? Please can you help me how I can get his comments. – Mutaz May 31 '18 at 22:36
• Maybe s/he deleted it? – Roman Jun 2 '18 at 8:17
• I hope not, but it seems you alright 😌 – Mutaz Jun 2 '18 at 9:06