# Speeding up double integration by nested Integrate1D

Hi, I am currently running some integrations. As they are over the complex plane, I started by modifying the given example by Igor to perform a double integration. While it works, it is stuck at being single-threaded. May I know if there is a way to use multi-threading to speed it up? I tried putting multithread keywords in it but they didn't help. The below is the full function: (To give more details, it it reproducing eq(6), the Wigner function, in this paper: https://arxiv.org/abs/1011.6668)
Function/c do2DintForWF(xmin,xmax,ymin,ymax,pWave)
Variable xmin,xmax,ymin,ymax
wave/d pWave

Variable/G globalXmin=xmin
Variable/G globalXmax=xmax
Variable/G globalY

variable/d/c integralT
multithread integralT = Integrate1d(outerInt,ymin,ymax,1,0,pWave)       // Romberg integration

return  integralT
End

Function/c innerInt(pWave1,inX)
wave/d pWave1
Variable inX

//I chose inX to be lambda_Real, globalY to be lambda_Imaginary

NVAR globalY=globalY

//obtaining parameters from pWave
variable m,n
n=pWave1[2]
m=pWave1[3]

variable x,p
x = pWave1[0]
p = pWave1[1]

//Breaking down the equation into parts
variable/d/c integral1
variable/c/d inB=cmplx(1,0), inC=cmplx(1,0), inH, inI, inJ, inK
variable k,l
if (m==0)
inB = cmplx(1,0)
else
multithread inB = (-cmplx(inX,-globalY))^m
endif
if (n==0)
inC = cmplx(1,0)
else
multithread inC = (cmplx(inX,globalY))^n
endif

inH = exp(-0.5*cmplx(inX,globalY)*cmplx(inX,-globalY))  //tested
inI = exp(cmplx(X,P)*cmplx(inX,-globalY))       //tested
inJ = exp(cmplx(X,-P)*cmplx(inX,globalY))       //tested

integral1 = inB*inC*inH*inI/inJ

return integral1

End

Function/c outerInt(pWave2,inY)
wave/d pWave2
Variable inY

NVAR globalY=globalY
globalY=inY
NVAR globalXmin=globalXmin
NVAR globalXmax=globalXmax

variable/d/c integral2
multithread integral2 = Integrate1D(innerInt,globalXmin,globalXmax,1,0,pWave2)  // Romberg integration

return integral2
End
In the past you posted questions relating to IP6 and IP7. Here the version of Igor that you are using is critical. In IP6 Integrate1d was NOT thread safe. I recommend using IP7 or IP8 for this task. In the later versions of Igor there is automatic multithreading that might kick in at some point. If you want to experiment with the threshold for the multi-threading read the documentation for the MultiThreadingControl operation. Also note that there is an Integrate2D operation that may be useful for your application. Note that when Integrate1D is automatically multithreaded you may not gain much by calling MultiThread inside the integrand function. A.G.
Igor wrote:
In the past you posted questions relating to IP6 and IP7. Here the version of Igor that you are using is critical. In IP6 Integrate1d was NOT thread safe. I recommend using IP7 or IP8 for this task. In the later versions of Igor there is automatic multithreading that might kick in at some point. If you want to experiment with the threshold for the multi-threading read the documentation for the MultiThreadingControl operation. Also note that there is an Integrate2D operation that may be useful for your application. Note that when Integrate1D is automatically multithreaded you may not gain much by calling MultiThread inside the integrand function. A.G.
Thanks for the note, I was using Igor 7 above. I was looking at also Integrate2D but thought that might not work if the function return were complex, maybe I misread it? It said "real-valued user define function" For thread-safe issue, I am still exploring it. Not knowing good enough about this, I tried to naively flag them thread-safe and it crushed Igor. I believe the major point where things got slow down was that, inside the innerInt function I had the "array raised to power n or m" which were unnecessarily repeated. The problem I had was that, inX and globalY were variables that I cannot simply raise them outside of the function like I could with my function in another post.
Indeed, Integrate2D applies to real functions but you might be able to write an analytic form that breaks your task into evaluating two real 2D integrals. If you were able to crash Igor in any way you should report it to support@wavemetrics.com. Clearly working with threads could sometimes lead to unexpected results but if you found a way that you can crash the application we would appreciate a report. A.G. WaveMetrics, Inc.
Igor wrote:
Indeed, Integrate2D applies to real functions but you might be able to write an analytic form that breaks your task into evaluating two real 2D integrals. If you were able to crash Igor in any way you should report it to support@wavemetrics.com. Clearly working with threads could sometimes lead to unexpected results but if you found a way that you can crash the application we would appreciate a report. A.G. WaveMetrics, Inc.
Sure, I will report it next time; I thought that was a trivial thing as I didn't know what I was doing. I do manage to crush it often when dealing with XOP and stuff, I will send a feedback next time.