# Multithread a for loop

Fri, 09/28/2018 - 10:21 am

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