1
$\begingroup$

I am trying to implement the Thomas algorithm using OpenCL.

__kernel void thomas(__global float *a_d, __global float *b_d, __global float *c_d, __global float *d_d, __global float *x_d, 
                                     __local float *shared, int system_size, int num_systems, int iterations)
{
    int thid = get_local_id(0);
    int blid = get_group_id(0);

    int delta = 1;

    __local float* a = shared;
    __local float* b = &a[system_size+1];
    __local float* c = &b[system_size+1];
    __local float* d = &c[system_size+1];
    __local float* cs = &d[system_size+1];
    __local float* ds = &cs[system_size+1];

    a[thid] = a_d[thid + blid * system_size];
    b[thid] = b_d[thid + blid * system_size];
    c[thid] = c_d[thid + blid * system_size];
    d[thid] = d_d[thid + blid * system_size];

    float tmp;

    barrier(CLK_LOCAL_MEM_FENCE);

    cs[0] = c[0] / b[0];
    ds[0] = d[0] / b[0];


    // Forward "warp"
    barrier(CLK_LOCAL_MEM_FENCE);

    if(thid>0)
    {
        tmp = (b[thid] - a[thid] * c[thid - 1]);
        cs[thid] = c[thid] / tmp;
        ds[thid] = (d[thid] - d[thid - 1] * a[thid]) / tmp;
    }
    barrier(CLK_LOCAL_MEM_FENCE);
    if(thid < system_size)
    {
        d[thid] = ds[thid] - cs[thid]*d[thid + 1];
    }
    barrier(CLK_LOCAL_MEM_FENCE);

    x_d[thid + blid * system_size] = d[thid];
}

This is giving me an incorrect result. The host code is working perfectly fine. It has been used on several kernels that I have been experimenting with. I the problem is arising because the threads are not executed in any order (?).

Has anyone got some advice?

$\endgroup$
4
$\begingroup$

The Thomas algorithm as proposed in most textbooks is inherently serial and anything that just executes the steps concurrently will in general lead to wrong results. You will need a clever modification to parallelize the algorithm which usually comes at the price of at least a logarithmic factor (w.r.t. the number of threads). However, I doubt that you will be able to ever beat the Thomas algorithm if the problem size does not justify parallelization.

Related discussion: Parallel vs Serial Thomas Algorithm

$\endgroup$
  • $\begingroup$ I see. In this case would it be better to use something like PCR? I have successfully implemented PCR. $\endgroup$ – HMPARTICLE Sep 3 '18 at 9:04

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.