parallel processing issue

Share icon

I want to extract and sum up "beams" from a 3D wave depending on a ROI or mask wave, which I do using the code below:

function GetSum(cube, mask)
    wave cube, mask
 
    Make/O/N=(DimSize(cube, 2)) total=0
        
    variable cx, cy 
    for(cx = 0; cx < DimSize(Cube, 0); cx += 1)         
            for(cy = 0; cy < DimSize(Cube, 1); cy += 1)
                    
                if (mask[cx][cy] ==1)           // mask[cx][cy] can be 1 or 0               
                    ImageTransform /Beam ={(cx), (cy)} getBeam cube
                    wave W_Beam
                    total+=W_Beam   
                endif       
            
            endfor
    endfor      
    
    KillWaves/Z W_Beam
end


For the size of my 3D wave this process takes about 2-3 seconds. To improve performance I copied and modified the parallel processing example from the help files to this:

Threadsafe function workerFunc(cube, total, mask, col)
    wave cube, total, mask
    variable col
    
    variable i
    for (i=0; i<DimSize(cube, 0);i+=1)
        
        if (mask[i][col] ==1)   
            ImageTransform /Beam ={(i), (col)} getBeam cube
            wave W_Beam
            total+=W_Beam   
        endif
    endfor
 
end
 
 
Function MT_GetSum(cube, mask)
    WAVE cube, mask
    
    Make/O/N=(DimSize(cube, 2)) total=0
    
    Variable ncol= DimSize(cube,1)
    Variable i,col,nthreads= ThreadProcessorCount
    Variable threadGroupID= ThreadGroupCreate(nthreads)
    
    for(col=0; col<ncol;)
        for(i=0; i<nthreads; i+=1)
        
            ThreadStart threadGroupID,i,workerFunc(cube,total, mask,col)
            
            col+=1
            if( col>=ncol )
                break
            endif
        endfor
        
        do
            Variable threadGroupStatus = ThreadGroupWait(threadGroupID,100)
        while( threadGroupStatus != 0 )
    endfor
    Variable dummy= ThreadGroupRelease(threadGroupID)
End



This gives a 60-70% speedup but it appears that wave "total" from the two methods does not contain identical values. The parallel version has slightly lower totals between 0-3% (0.25 avg) compared to the single thread version.
Any ideas what I'm doing wrong?

I don't have a lot of time but from a quick look at your post:

ChrLie wrote: 

Threadsafe function workerFunc(cube, total, mask, col)
    wave cube, total, mask
    variable col
    
    variable i
    for (i=0; i<DimSize(cube, 0);i+=1)
        
        if (mask[i][col] ==1)   
            ImageTransform /Beam ={(i), (col)} getBeam cube
            wave W_Beam
            total+=W_Beam   
        endif
    endfor
 
end



Think about what happens in your parallel calculation. Igor spawns a bunch of threads and these start crunching away, independently and regardless of what other threads may be doing. This also means that if different threads write to the same memory, each will modify this memory as it pleases and the result is unpredictable. Read the section called "Example atomic operation" here to understand what that means with respect to your problem.

Igor tries to prevent these problems from occurring by giving each thread its own, private data folder hierarchy, such that all waves created by the thread are accessible only by that thread. The exception, of course, are function arguments. In other words, every thread refers to the same "total" wave.

So when you do this:
total+=W_Beam  

The result is undefined because multiple threads can modify it simultaneously.

Your only option is to provide each thread with its own copy of total (so if there are 8 threads, create 8 copies of total, all initialized to zero), and to add all of these waves together into the final result after the threads have completed.


Also, the loop in which you call ThreadStart has issues, but I don't have time now so maybe someone else will chip in on that.

October 19, 2012 at 07:36 am - Permalink

Looking at this again, your ThreadStart loop does look OK. I misread it due to reading it in a hurry, sorry. However, my comment above still stands.

I played with this a bit:
function GetSum(cube, mask) // this is your original function
    wave cube, mask
 
    Make/O/N=(DimSize(cube, 2))/D total=0
 
    variable cx, cy 
    for(cx = 0; cx < DimSize(Cube, 0); cx += 1)         
            for(cy = 0; cy < DimSize(Cube, 1); cy += 1)
 
                if (mask[cx][cy] ==1)           // mask[cx][cy] can be 1 or 0               
                    ImageTransform /Beam ={(cx), (cy)} getBeam cube
                    wave W_Beam
                    total+=W_Beam[p]
                endif       
 
            endfor
    endfor      
 
    KillWaves/Z W_Beam
end
 
Function GetSum2(cube, mask)    // MatrixOP-based
    wave cube, mask
    
    Make/O/N=(DimSize(cube, 2))/D total=0
    
    MatrixOP /FREE Multiplied = cube * mask
    
    variable i
    for (i = 0; i < DimSize(Multiplied, 2); i+=1)
        ImageTransform /P=(i) sumPlane, Multiplied
        total[i] += V_value
    endfor
End
 
function GetSum3(cube, mask)    // MultiThreaded
    wave cube, mask
 
    Make/O/N=(DimSize(cube, 2))/D total=0
    
    MultiThread total = WorkerFunc3(cube, mask, p)
end
 
ThreadSafe Function WorkerFunc3(cube, mask, plane)
    wave cube, mask
    variable plane
    
    variable i, j, total = 0
    for (i = 0; i < DimSize(cube, 0); i+=1)
        for (j = 0; j < DimSize(cube, 1); j+=1)
            if (mask[i][j] == 1)
                total += cube[i][j][plane]
            endif
        endfor
    endfor
    
    return total
End
 
function GetSum4(cube, mask)    // MultiThreaded + MatrixOP
    wave cube, mask
 
    Make/O/N=(DimSize(cube, 2))/D total=0
    
    MultiThread total = WorkerFunc4(cube, mask, p)
end
 
ThreadSafe Function WorkerFunc4(cube, mask, plane)
    wave cube, mask
    variable plane
    
    MatrixOP /FREE W_Sum = sum(mask * cube[][][plane])
    return W_Sum[0]
End
 
Function DoTiming()
    Make /o/n=(128,128) /D M_Mask = round(enoise(0.5) + 0.5)
    Make /o/n=(128,128, 512) /D/O M_cube = enoise(1)
    
    variable timerRef, time1, time2, time3, time4
    
    timerRef = StartMSTimer
    GetSum(m_cube, m_mask)
    time1 = StopMSTimer(timerRef) / 1e6
    
    timerRef = StartMSTimer
    GetSum2(m_cube, m_mask)
    time2 = StopMSTimer(timerRef) / 1e6
    
    timerRef = StartMSTimer
    GetSum3(m_cube, m_mask)
    time3 = StopMSTimer(timerRef) / 1e6
    
    timerRef = StartMSTimer
    GetSum4(m_cube, m_mask)
    time4 = StopMSTimer(timerRef) / 1e6
    
    Printf "Original: %f, MatrixOP: %f s, MultiThread: %f s, MatrixOP + MultiThread: %f s\r", time1, time2, time3, time4
End


Result on my computer:
•DoTiming()
  Original: 0.343414 s, MatrixOP: 0.056957 s, MultiThread: 0.389801 s, MatrixOP + MultiThread: 0.008289 s


In my implementation, multithreading of the 'naive' Igor code actually leads to a slowdown. Using MatrixOP, however, is much faster than the naive code and multithreaded MatrixOP appears to be the performance winner.

Note the dirty trick that I'm using here – I make use of the fact that your mask contains only zero or one, and therefore the undesired values can be removed simply by multiplying the data with the mask (but of course you need to make a copy of the data to do that). Doing this multiplication is much faster than the explicit Igor code because it runs entirely in machine code, and computers are very good at multiplying lots of numbers.

ChrLie wrote:
The parallel version has slightly lower totals between 0-3% (0.25 avg) compared to the single thread version.

I think that this is probably due to a mistake in your code. But be wary of using single-precision waves when not necessary. Use doubles (Make /D) when practical.

I'd be interested to hear what the performance difference is on your dataset.

October 22, 2012 at 01:02 pm - Permalink

741 wrote:

Note the dirty trick that I'm using here – I make use of the fact that your mask contains only zero or one, and therefore the undesired values can be removed simply by multiplying the data with the mask (but of course you need to make a copy of the data to do that). Doing this multiplication is much faster than the explicit Igor code because it runs entirely in machine code, and computers are very good at multiplying lots of numbers.


GetSum4() is awesome!

Thinking of "planes" instead of "beams" clearly makes a difference. The result for my cube:

  Original: 2.831555, MatrixOP: 0.350114 s, MultiThread: 6.146490 s, MatrixOP + MultiThread: 0.097558 s


As for differences in totals in the initial post; I somehow had the feeling that all threads working on one wave (or memory) might be the problem but your solution is much more elegant anyway!

Thanks a lot!

October 23, 2012 at 12:42 am - Permalink