如何优化CSR格式矩阵与向量乘法的并行实现？

问题概述 #

你需要并行化一段CSR格式对称矩阵与向量的乘法代码，当前通过手动拆分任务到多个局部数组再合并的方案过于繁琐，寻求更优实现思路。

原串行实现代码：

/// <summary>
/// q = A * d, with A stored in CSR format.
/// </summary>
public static double[] Multiplication(int order, double[] d, List<double> A, List<int> J, int[] I)
{ 
    double[] q = new double[order];

    for (int i = 0; i < order; i++)
    {                    
       int indexFirst = I[i];

       q[i] += A[indexFirst] * d[J[indexFirst]];

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

          q[i] += a * d[col];

          q[col] += a * d[i];
       }
   }

   return q;
}

当前手动拆分的并行实现：

Parallel.For(0, cpuCount, delegate(int i)
{
    MultiplySection(startIndices[i], endIndices[i], d, A, J, I, qSection[i]);
});

double[] q = new double[order];

for (int j = 0; j < cpuCount; j++)
{
    for (int i = startIndices[j]; i < order; i++)
    {
        q[i] += qSection[j][i];
    }
}

private static void MultiplySection(int startIndex, int endIndex, double[] d, List<double> A, List<int> J, int[] I, double[] q)
{
    for (int i = startIndex; i < endIndex; i++)
    {
        int indexFirst = I[i];

        q[i] += A[indexFirst]*d[J[indexFirst]];

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

            q[i] += a*d[col];

            q[col] += a*d[i];
        }
    }
}

更优实现方案 #

1. 利用Parallel.For的本地存储特性简化实现 #

Parallel.For支持线程本地存储，可以让每个线程先在自己的局部数组中完成累加，最后再一次性合并到全局数组，避免手动拆分任务和管理多个局部数组的繁琐操作，同时减少线程间的竞争。

实现代码：

public static double[] ParallelMultiplication(int order, double[] d, List<double> A, List<int> J, int[] I)
{
    double[] q = new double[order];

    Parallel.For(
        // 遍历范围：矩阵的所有行
        0, 
        order,
        // 线程本地初始化：每个线程创建独立的局部q数组
        () => new double[order],
        // 循环体：在局部数组上执行与串行逻辑一致的累加操作
        (row, state, localQ) =>
        {
            int rowStart = I[row];
            // 处理对角线元素（或行首元素）
            localQ[row] += A[rowStart] * d[J[rowStart]];

            // 处理行内其余非零元素
            for (int idx = rowStart + 1; idx < I[row + 1]; idx++)
            {
                int col = J[idx];
                double val = A[idx];
                localQ[row] += val * d[col];
                localQ[col] += val * d[row];
            }
            return localQ;
        },
        // 合并逻辑：将每个线程的局部数组结果累加到全局q
        localQ =>
        {
            lock (q) // 加锁避免合并时的线程竞争
            {
                for (int i = 0; i < order; i++)
                {
                    q[i] += localQ[i];
                }
            }
        });

    return q;
}

2. 优化：利用矩阵对称性减少竞争（可选） #

从代码逻辑可以看出，你处理的是对称矩阵（仅存储一半元素，通过q[col] += val*d[row]完成对称部分的计算）。可以调整遍历逻辑，只处理row <= col的元素对，避免重复操作：

• 遍历行时，仅处理当前行中col >= row的元素
• 对col > row的元素，仅在处理row行时更新q[row]和q[col]，不再在col行重复处理

这种方式能减少一半的计算量，同时彻底避免线程间对同一q元素的竞争，不过需要确保CSR存储的是矩阵的上三角/下三角部分。

示例代码（假设CSR存储下三角，即col <= row）：

public static double[] SymmetricParallelMultiplication(int order, double[] d, List<double> A, List<int> J, int[] I)
{
    double[] q = new double[order];

    Parallel.For(0, order, row =>
    {
        int rowStart = I[row];
        for (int idx = rowStart; idx < I[row + 1]; idx++)
        {
            int col = J[idx];
            double val = A[idx];
            q[row] += val * d[col];
            if (col != row)
            {
                // 仅在col < row时更新对称位置，避免重复计算
                q[col] += val * d[row];
            }
        }
    });

    return q;
}

注：此方案无需本地存储和锁，因为每个(row, col)对仅被处理一次，线程间不会竞争同一q元素的更新操作。

方案优势 #

• 简化代码结构：无需手动拆分任务、管理局部数组索引，Parallel.For自动处理线程分配和任务调度
• 降低竞争开销：线程本地存储方案将竞争延迟到最后合并阶段，对称矩阵优化方案则彻底消除竞争
• 可读性强：逻辑与串行代码高度一致，便于维护和调试

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