能否在Parallel.For的Local Finally段使用Kahan求和算法？

首先明确：并行场景下完全可以用卡汉求和减少数值误差，但你的实现存在核心逻辑错误，导致结果和串行版本有差异。

你的代码核心错误点 #

1. 共享全局补偿数组c是致命问题
   卡汉求和的补偿项c是和当前累加器的状态强绑定的，必须每个线程维护自己的局部补偿项，不能多个线程共用全局的c。你现在在lock里修改全局c，会导致不同线程的补偿计算互相干扰，完全破坏了卡汉求和的逻辑。

2. 局部累加阶段没有用卡汉求和
   你的并行循环里，局部数组q只是普通的+=累加，没有应用卡汉求和，局部累加已经产生了数值误差，最后合并时再用卡汉也无法挽回这部分误差。

3. 合并逻辑错误
   合并局部q到全局qTotal时，使用全局c计算补偿，这和卡汉求和的要求不符——合并过程也需要基于全局当前的补偿状态来正确累加局部结果。

修正方案 #

需要让每个线程维护自己的局部累加数组+局部补偿数组，局部累加时就用卡汉求和，最后合并时再用卡汉的方式把局部结果合并到全局。

修正后的代码示例 #

// 定义局部状态：每个线程维护自己的累加数组和补偿数组
class LocalState
{
    public double[] Q { get; set; }
    public double[] C { get; set; }
}

double[] qTotal = new double[order];
double[] globalC = new double[order]; // 全局补偿数组

Parallel.For(0,
    order,
    () => new LocalState // LocalInit：每个线程初始化自己的状态
    {
        Q = new double[order],
        C = new double[order]
    },
    (i, loopState, localState) => // Main Body：局部累加用卡汉求和
    {
        int indexFirst = I[i];

        // 对q[i]做卡汉累加
        double y = A[indexFirst] * d[J[indexFirst]] - localState.C[i];
        double t = localState.Q[i] + y;
        localState.C[i] = (t - localState.Q[i]) - y;
        localState.Q[i] = t;

        for (int j = indexFirst + 1; j < I[i + 1]; j++)
        {
            int col = J[j];
            double a = A[j];
            double val = a * d[col];

            // 累加q[i]
            y = val - localState.C[i];
            t = localState.Q[i] + y;
            localState.C[i] = (t - localState.Q[i]) - y;
            localState.Q[i] = t;

            val = a * d[i];
            // 累加q[col]
            y = val - localState.C[col];
            t = localState.Q[col] + y;
            localState.C[col] = (t - localState.Q[col]) - y;
            localState.Q[col] = t;
        }
        return localState;
    },
    localState => // Local Finally：合并局部状态到全局，同样用卡汉求和
    {
        lock (qTotal)
        {
            for (int i = 0; i < order; i++)
            {
                double y = localState.Q[i] - globalC[i];
                double t = qTotal[i] + y;
                globalC[i] = (t - qTotal[i]) - y;
                qTotal[i] = t;
            }
        }
    });

额外说明 #

并行版本的结果和串行版本完全一致是很难的，因为浮点数累加的顺序会影响结果（即使都是卡汉求和）。但修正后的代码会把数值误差控制在和串行卡汉求和相当的范围，差异会远小于你之前的实现。

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