Speeding up double integration by nested Integrate1D

Sandbo
Posts: 39
Joined: 2017-05-10
Location: Canada

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


Igor's picture
Posts: 817
Joined: 2007-06-29
Location: United States

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.


Sandbo
Posts: 39
Joined: 2017-05-10
Location: Canada

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.


Igor's picture
Posts: 817
Joined: 2007-06-29
Location: United States

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.


Sandbo
Posts: 39
Joined: 2017-05-10
Location: Canada

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.


Back to top