Thomas 算法内核 OpenCL

计算科学 线性代数 矩阵 开放式
2021-11-30 20:58:25

我正在尝试使用 OpenCL 实现 Thomas 算法。

__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];
}

这给了我一个不正确的结果。主机代码运行良好。它已在我一直在试验的几个内核上使用。我出现问题是因为线程没有按任何顺序执行(?)。

有人给点建议吗?

1个回答

大多数教科书中提出的 Thomas 算法本质上是串行的,任何只是同时执行这些步骤的东西通常都会导致错误的结果。您将需要一个巧妙的修改来并行化通常以至少一个对数因子(线程数)为代价的算法。但是,如果问题大小不能证明并行化是合理的,我怀疑您是否能够击败 Thomas 算法。

相关讨论:并行与串行 Thomas 算法