Crout分解并行循环结果错误问题求助

问题根源：上三角循环的迭代间数据依赖 #

在上三角的j循环（for (j = 0; j < i; j++)）中，存在迭代间的真依赖（流依赖）：

• 对于j1 < j2，j2的计算会读取matx[j1 * size + i]（当k = j1时），而这个值是由j1迭代写入的。
• 串行执行时，j从小到大依次处理，确保j1先完成写入，j2再读取正确的更新后值。
• 并行化后，j的迭代可能乱序执行，j2可能在j1完成写入前就读取matx[j1 * size + i]，导致使用原始未更新的值计算，最终结果偏差。

下三角循环可以安全并行的原因：
每个j的计算仅依赖k < i的matx[k * size + i]（这些值在当前i循环的上三角处理前已确定），以及matx[j * size + k]（属于已处理的列，数据稳定）；且不同j之间仅写入各自的matx[j * size + i]，无交叉依赖，因此并行不会引发数据错误。

修正方案 #

方案1：仅并行下三角（最直接可靠） #

保留上三角的串行执行，只对无依赖的下三角部分并行，这是保证正确性的最优选择：

for (i = 0; i < size; i++)
{
    // Upper Triangle - 保持串行，避免迭代间依赖问题
    for (j = 0; j < i; j++)
    {
        q = matx[j * size + i];
        for (k = 0; k < j; k++)
        {
            q -= matx[j * size + k] * matx[k * size + i];
        }
        matx[j * size + i] = q;
    }

    // Lower Triangle - 安全并行
#pragma omp parallel for private(j,k,q)
    for (j = i; j < size; j++)
    {
        q = matx[j * size + i];
        for (k = 0; k < i; k++)
        {
            q -= matx[j * size + k] * matx[k * size + i];
        }
        matx[j * size + i] = q;
    }
}

方案2：用OpenMP任务依赖实现上三角并行（复杂度较高） #

如果追求最大化并行度，可使用OpenMP的任务依赖机制强制j迭代按顺序执行，但任务调度会带来额外开销，仅适合大规模矩阵：

for (i = 0; i < size; i++)
{
    // Upper Triangle - 任务依赖保证j的执行顺序
#pragma omp parallel
#pragma omp single
    {
        for (j = 0; j < i; j++)
        {
            // 声明任务依赖：当前任务需等待所有k<j的matx[k*size+i]写入完成，再写入matx[j*size+i]
#pragma omp task depend(in: matx[0:j*size+i]) depend(out: matx[j*size+i]) private(k,q)
            {
                q = matx[j * size + i];
                for (k = 0; k < j; k++)
                {
                    q -= matx[j * size + k] * matx[k * size + i];
                }
                matx[j * size + i] = q;
            }
        }
#pragma omp taskwait // 等待所有上三角任务完成后再处理下三角
    }

    // Lower Triangle - 并行
#pragma omp parallel for private(j,k,q)
    for (j = i; j < size; j++)
    {
        q = matx[j * size + i];
        for (k = 0; k < i; k++)
        {
            q -= matx[j * size + k] * matx[k * size + i];
        }
        matx[j * size + i] = q;
    }
}

注意：该方案需要编译器支持OpenMP 3.0及以上版本，且需根据实际矩阵规模测试性能收益。

内容的提问来源于stack exchange，提问作者MikeBushell