Speeding up double integration by nested Integrate1D

Share icon

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.

June 11, 2018 at 10:13 am - Permalink

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.

June 11, 2018 at 11:17 am - Permalink

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.

June 12, 2018 at 01:01 pm - Permalink

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.

June 12, 2018 at 06:41 pm - Permalink