Multithread a for loop
Hi,
Is there a general way to multithread a for-loop if they are embarrassingly parallel?
I asked about this earlier but I think I didn't understand the problem sufficiently to continue. Now that I have better ideas about what I am doing, it would be great if I can have your inputs.
Problem:
I am doing a double (complex) integral of a function to create a color map. The inputs are essentially x and y coordinates and a 5-by-5 matrix. For each pair of (x,y) combination, a 2D integration over a complex variable "lambda" is performed using the function Integrate1D, which is equation (6) of this paper:
https://pdfs.semanticscholar.org/a5fe/278566175bb7852aa107d5d213c0e4030…
The difficulty was that I found no easy way to multithread the Integrate1D function, as Integrate1D doesn't seem to allow it after all. Is there another way to enable the multithreading on a higher level? The outermost for-loop which loop through x and y can be done independently from each other as pair.
Appreciated if you can give me some ideas.
I can post the code here but they are 4 nested functions so it's not easy to display them, let me know if they will be useful.
The Integrate1D operation is listed as automatically multithreaded, and also threadsafe. Automatic threading may depend on your integrand function being threadsafe- you need to qualify the function declaration with the Threadsafe keyword.
But if you are doing a number of separate integrations for each cell of a matrix, you may get a better boost by running integrations in threads. If you can write the integration as a function returning a value to assign to the matrix wave, you may be able to do it with the Multithread keyword.
To learn more about threading in Igor, read our help starting with DisplayHelpTopic "ThreadSafe Functions and Multitasking"
September 28, 2018 at 02:14 pm - Permalink
In reply to The Integrate1D operation is… by johnweeks
Thanks for the ideas, it clears a few problem I have. I didn't know Integrate1D is automatically multithreaded.
Naively, I tried to just flag them threadsafe but it crashed Igor. I then found that I have been doing things that Threadsafe disallowed. I will go through the help file in more details.
Update, going through the code, I think one problem I have (from the example as well) is that I need to update the globalY coordinate in the outerloop, which is updated by the outer Integrate1D. Here is the code:
do2DintForWF is the first function which perfoms a 2D integration for one coordinate.
It does the first Integrate1d and call the outerInt.
In outerInt, the globalY is being updated by the inY produced when outerInt is being called. This globalY is then used when doing the second Integrate1D with innerInt. This prevents me from flagging the outerInt function as Threadsafe (or it will crash).
How should I solve this?
Variable xmin,xmax,ymin,ymax
wave/d pWave
Variable/G globalXmin=xmin
Variable/G globalXmax=xmax
Variable/G globalY
variable/d/c integralT
integralT = Integrate1d(outerInt,ymin,ymax,2,0,pWave) // Romberg integration
return integralT
End
Function/c outerInt(pWave2,inY)
wave/d pWave2
Variable inY
NVAR globalY//=globalY
globalY=inY//<--prevents outerInt from becoming Threadsafe
NVAR globalXmin//=globalXmin
NVAR globalXmax//=globalXmax
variable/d/c integral2
integral2 = Integrate1D(innerInt,globalXmin,globalXmax,2,0,pWave2) // Romberg integration
return integral2
End
Threadsafe Function/c innerInt(pWave1,inX)
wave/d pWave1
Variable inX
//I chose inX to be lambda_R, globalY to be lambda_I
NVAR globalY//=globalY
//make pWave more reable
variable m,n
n=pWave1[2]
m=pWave1[3]
// nfac=factorial(n)
// mfac=factorial(m)
variable x,p
x = pWave1[0]
p = pWave1[1]
// variable/d AMatReal, AMatImag
// AMatReal = pWave1[4]
// AMatImag = pWave1[5]
variable/d/c integral1
//integral = inX
//break down the integral1
variable/c/d inA, inB=cmplx(1,0), inC=cmplx(1,0), inH, inI, inJ, inK
//inA = cmplx(AMatReal,AMatImag)/(pi^2*nfac*mfac)
//use if statement to separate out the m = 0 case
variable k,l
if (m==0)
inB = cmplx(1,0)
else
inB = (-cmplx(inX,-globalY))^m
// for (k=0;k<m;k=k+1)
// inB = (-cmplx(inX,-globalY))*inB
// endfor
endif
if (n==0)
inC = cmplx(1,0)
else
inC = (cmplx(inX,globalY))^n
// for (l=0;l<n;l=l+1)
// inC = (cmplx(inX,globalY))*inC
// endfor
endif
inH = exp(-0.5*cmplx(inX,globalY)*cmplx(inX,-globalY)) //tested
// inH = exp((-0.5*(inx^2+globalY^2))+(cmplx(X,P)*cmplx(inX,-globalY))-(cmplx(X,-P)*cmplx(inX,globalY)))
inI = exp(cmplx(X,P)*cmplx(inX,-globalY)) //tested
inJ = exp(cmplx(X,-P)*cmplx(inX,globalY)) //tested
integral1 = inB*inC*inH*inI/inJ
//integral1 = -cmplx(inX,-globalY)*cmplx(inX,globalY)*inH*inI/inJ
return integral1
//return cmplx(AMatReal,AMatImag)/(pi^2*nfac*mfac)*(cmplx(-inX,globalY))^m*(cmplx(inX,globalY))^n*exp(-0.5*(inx^2+globalY^2))*exp(cmplx(X,P)*cmplx(inX,-globalY))*exp(cmplx(X,-P)*cmplx(inX,globalY))
End
September 28, 2018 at 09:14 pm - Permalink
In reply to johnweeks wrote: The… by Sandbo
Hello Sandbo,
Any reason for not using Integrate2D?
AG
October 1, 2018 at 08:39 am - Permalink
In reply to Hello Sandbo, Any reason… by Igor
Maybe it didn't work with complex number?
I think you gave me this suggestion last time, but I haven't really tried to make it work that way yet.
If I want to use it with complex number, how should I start?
October 1, 2018 at 09:26 am - Permalink
In order to deal with the complex number you would have to express your integrand in terms of real and imaginary parts. Not so easy unless you are considering a limited range of m,n. Also, before trying any of these it would help to know what are alpha and alpha*.
October 1, 2018 at 01:03 pm - Permalink
In reply to In order to deal with the… by Igor
If you wouldn't mind taking a look, I am trying to implement equ (6) in this paper:
https://arxiv.org/pdf/1011.6668.pdf
alpha = x+i*p is a complex number as the input argument to equation (6) (sorry I forgot to include a definition for the pwave, the parameter wave). Alpha represent the coordinate (x,p) on a complex plane.
Then, equ (6) integrates for a given coordinate (as x and p, with alpha=x+i*p) for a certain xrange and yrange.
My bad to use the x again, let's called it x' here. Then xrange (x') and yrange (y') represent the integration variable of the complex number lambda (lambda = x'+i*y').
In short, my Igor integration will take the parameter alpha=x+i*p which is a constant in the integration, then integrate over a complex variable lambda=x'+i*y' with equation (6).
Let me know if there is anything else I can clarify, thank you.
October 1, 2018 at 03:22 pm - Permalink
1. Please start by changing x and p in innerInt() to different variable names. This may save you from unintended consequences.
2. It may be trivial but I think it is important to re-write equation (6) so that all the factors that are independent of the integration are out of the integral. You do not want even a single multiplication inside any of the integration functions if it can be computed outside. You are also potentially setting twice both inB=cmplx(1,0), inC=cmplx(1,0).
3. Before actually attempting the integration, it may be worth your time to check in some tables if the integration can't be done in closed form for some values of m, n and alpha. Even if you can find a closed form solution for one point, it can be used to check your integration and IMO worth the trouble.
October 2, 2018 at 12:51 pm - Permalink