Thomas Algorithm Kernel OpenCL

Question

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?

Answer 1

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